Charles Gueunet
email: charles.gueunet@kitware.com

Pierre Fortin
email: pierre.fortin@lip6.fr

Julien Jomier
email: julien.jomier@kitware.com

Julien Tierny
email: julien.tierny@lip6.fr

Contour Forests: Fast Multi-threaded Augmented Contour Trees

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

As scientific data-sets become more intricate and larger in size, advanced data analysis algorithms are needed for their efficient visualization and interactive exploration. For scalar field visualization, topological analysis techniques have shown to be practical solutions in various contexts by enabling the concise and complete capture of the structure of the input data into high-level topological abstractions such as contour trees [START_REF] Carr | Computing contour trees in all dimensions[END_REF], Reeb graphs [START_REF] Pascucci | Robust on-line computation of Reeb graphs: simplicity and speed[END_REF][START_REF] Biasotti | Reeb graphs for shape analysis and applications[END_REF][START_REF] Tierny | Loop surgery for volumetric meshes: Reeb graphs reduced to contour trees[END_REF], or Morse-Smale complexes [START_REF] Gyulassy | A practical approach to Morse-Smale complex computation: scalabity and generality[END_REF][START_REF] Floriani | Morse complexes for shape segmentation and homological analysis: discrete models and algorithms[END_REF]. Such topological abstractions are fundamental data-structures that enable the development of advanced data analysis, exploration and visualization techniques, including for instance: small seed set extraction for fast isosurface traversal [START_REF] Van Kreveld | Contour trees and small seed sets for isosurface traversal[END_REF][START_REF] Carr | Simplifying flexible isosurfaces using local geometric measures[END_REF], feature tracking [START_REF] Sohn | Time varying contour topology[END_REF], data-summarization [START_REF] Pascucci | Multi-resolution computation and presentation of contour trees[END_REF], transfer function design for volume rendering [START_REF] Weber | Topology-controlled volume rendering[END_REF], similarity estimation [START_REF] Thomas | Multiscale symmetry detection in scalar fields by clustering contours[END_REF], or application-driven segmentation and analysis tasks [START_REF] Laney | Understanding the structure of the turbulent mixing layer in hydrodynamic instabilities[END_REF][START_REF] Gyulassy | Topologically Clean Distance Fields[END_REF][START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF][START_REF] Gyulassy | Stability of dissipation elements: A case study in combustion[END_REF][START_REF] Gyulassy | Interstitial and interlayer ion diffusion geometry extraction in graphitic nanosphere battery materials[END_REF].

However, with the ongoing development of computational resources on the one hand and of acquisition devices on the other, the resolution of the geometrical domains on which scalar fields are defined is continuously increasing. This resolution increase yields several technical challenges for topological data analysis, including that of computation time efficiency. In particular, to enable truly interactive exploration sessions, highly efficient algorithms are required both for the computation of topological abstractions as well as for their interactive manipulation (i.e. topological simplification for instance). A natural direction towards the improvement of the time efficiency of topological data analysis is parallelism, as all commodity hardware now embeds processors with multiple cores.

However, most topological analysis algorithms are originally intrinsically sequential as they often require a global view on the data.

Regarding the contour tree -a fundamental topology-based data structure in scalar field visualization -several algorithms have been proposed for its parallel computation [START_REF] Pascucci | Parallel computation of the topology of level sets[END_REF][START_REF] Maadasamy | A hybrid parallel algorithm for computing and tracking level set topology[END_REF][START_REF] Acharya | A parallel and memory efficient algorithm for constructing the contour tree[END_REF]. However, these algorithms only compute non-augmented contour trees [START_REF] Carr | Computing contour trees in all dimensions[END_REF], trees that only represent the connectivity evolution of the level-sets, and not the corresponding data-segmentation (i.e. the arcs are not augmented with regular vertices). While such non-augmented trees enable some of the traditional visualization applications of the contour tree, they do not enable them all. For instance, they do not readily support seed set extraction for fast isosurface extraction or topology based data segmentation. In both examples, additional postprocessing will be needed to support these applications. Moreover, fully augmenting in a post-process non-augmented trees is a non trivial task, for which no linear-time algorithm has ever been documented to the best of our knowledge.

This paper addresses this problem by presenting a fast, shared memory multi-threaded algorithm for the computation of augmented contour trees on tetrahedral meshes. Such a tree augmentation makes our output data-structures generic application-wise and enables the full extent of contour tree based applications, including data segmentation. Extensive experiments demonstrate the efficiency of our approach in comparison to a reference sequential implementation (libtourtre [START_REF] Dillard | libtourtre: A contour tree library[END_REF]), despite the strong memory-bound aspect of the graph traversal tasks [START_REF] Asanovic | The landscape of parallel computing research: a view from berkeley[END_REF] involved in augmented contour tree computation, and despite other limitations that we detail. We illustrate the utility of our approach with specific use cases for the interactive exploration of hierarchies of topology-based data segmentations that were enabled by our algorithm. We also provide a lightweight VTK-based C++ reference implementation of our approach for reproduction purposes.

Related work

The notion of contour tree has first been introduced in Computer Science by Boyell and Ruston [START_REF] Boyell | Hybrid techniques for real-time radar simulation[END_REF]. Algorithms for their efficient computation have first been investigated for 2D domains [START_REF] Van Kreveld | Contour trees and small seed sets for isosurface traversal[END_REF], then for 3D domains [START_REF] Tarasov | Construction of contour trees in 3d in o(n log n) steps[END_REF] and last for domains of arbitrary dimension [START_REF] Carr | Computing contour trees in all dimensions[END_REF], with an algorithm that is simple to implement, efficient in practice and with optimal time complexity. In particular, this algorithm allows for the computation of both augmented and non-augmented contour trees. An open source reference implementation (libtourtre [START_REF] Dillard | libtourtre: A contour tree library[END_REF]) of this algorithm is provided by Scott Dillard. Chiang et al. [START_REF] Chiang | Simple and optimal outputsensitive construction of contour trees using monotone paths[END_REF] later presented an output-sensitive algorithm for the computation of non-augmented contour trees based on monotone paths, where the arcs of the intermediate trees (called join and split trees, see Sec. 2) were evaluated by considering monotone paths connecting the critical points of the input scalar field. Applications of the contour tree in data analysis and visualization include small seed set extraction for isosurface traversal [START_REF] Van Kreveld | Contour trees and small seed sets for isosurface traversal[END_REF], topological simplification of isosurfaces [START_REF] Carr | Simplifying flexible isosurfaces using local geometric measures[END_REF], feature tracking [START_REF] Sohn | Time varying contour topology[END_REF], transfer function design for volume rendering [START_REF] Weber | Topology-controlled volume rendering[END_REF], similarity estimation [START_REF] Thomas | Multiscale symmetry detection in scalar fields by clustering contours[END_REF], or applicationdriven segmentation and analysis tasks [START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF]. Note that all of the applications mentioned above require the augmented contour tree as they rely on the identification of the sets of regular vertices mapping to each arc of the contour tree, to find isosurface traversal seeds, to compute geometrical measures on the volume regions corresponding to each arc of the contour tree, to measure overlaps between such regions, to compare geometrical measures on such regions or to simply extract these regions for visualization purpose.

Pascucci and Cole-McLaughlin introduced the first parallel algorithm for contour tree computation [START_REF] Pascucci | Parallel computation of the topology of level sets[END_REF]. However, this algorithm does not compute contour trees that are augmented with regular vertices 1 . Moreover, the algorithm proceeds with a hierarchical domain decomposition. This means that fewer and fewer threads will be active along the execution, eventually finishing with a nonnegligible sequential pass. Also, while this paper describes a parallel algorithm for the computation of intermediate trees (the join and split trees, see Sec. 2), the authors still use a sequential pass to combine these two trees into the final output contour tree.

Maadasamy et al. [START_REF] Maadasamy | A hybrid parallel algorithm for computing and tracking level set topology[END_REF] introduced a multi-threaded variant of the output-sensitive algorithm by Chiang et al. [START_REF] Chiang | Simple and optimal outputsensitive construction of contour trees using monotone paths[END_REF], which results in good scaling performances on tetrahedral meshes. However, we note that, in practice, the sequential version of this algorithm is up to three times slower than the reference implementation (libtourtre [START_REF] Dillard | libtourtre: A contour tree library[END_REF]) of the algorithm by Carr et al. [START_REF] Carr | Computing contour trees in all dimensions[END_REF] (see Tab. 1 in [START_REF] Maadasamy | A hybrid parallel algorithm for computing and tracking level set topology[END_REF]). This only yields eventually speedups between 1.6 and 2.8 with regard to libtourtre [START_REF] Dillard | libtourtre: A contour tree library[END_REF] on a 8-core CPU [START_REF] Maadasamy | A hybrid parallel algorithm for computing and tracking level set topology[END_REF]. We suspect that these moderate speedups over libtourtre are due to the lack of efficiency of the sequential algorithm based on monotone paths by Chiang et al. [START_REF] Chiang | Simple and optimal outputsensitive construction of contour trees using monotone paths[END_REF] in comparison to that of Carr et al. [START_REF] Carr | Computing contour trees in all dimensions[END_REF]. Indeed, from our experience, although the extraction of the critical points of the field is a local operation [START_REF] Banchoff | Critical points and curvature for embedded polyhedral surfaces[END_REF], we found in practice that its overall computation time is often larger than that of the contour tree itself. Moreover, this algorithm triggers monotone path computations for each saddle point [START_REF] Chiang | Simple and optimal outputsensitive construction of contour trees using monotone paths[END_REF], even if it does not yield branching in the contour tree (which induces unnecessary computations). Note also that the approach by Maadasamy et al. [START_REF] Maadasamy | A hybrid parallel algorithm for computing and tracking level set topology[END_REF] only parallelizes the construction of the join and split trees and not their combination into the final contour tree, which is done sequentially. Finally, since it connects critical points through monotone paths, this algorithm does not visit all the vertices of the input mesh. Thus it cannot provide a contour tree-based data segmentation and does not produce an augmented contour tree. Acharya and Natarajan [START_REF] Acharya | A parallel and memory efficient algorithm for constructing the contour tree[END_REF] specialized and improved this approach for the special case of regular grids. In this paper, we focus however on tetrahedral meshes because of the genericity of this mesh representation (any mesh can be triangulated).

Morozov and Weber [START_REF] Morozov | Distributed merge trees[END_REF][START_REF] Morozov | Distributed contour trees[END_REF] and Landge et al. [START_REF] Landge | In-situ feature extraction of large scale combustion simulations using segmented merge trees[END_REF] presented three approaches for contour tree-based visualization in a distributed environment, with minimal inter-node communications. However, these approaches focus more on the reduction of the communication between the processes than on the efficient computation on a single shared memory node as we do here with the target of an efficient interactive exploration in mind.

Contributions

This paper presents the following new contributions:

1. a fast, shared memory multi-threaded algorithm for the computation of augmented contour trees on tetrahedral meshes;

2. a lightweight VTK-based C++ reference implementation of our approach for reproduction purposes.

PRELIMINARIES

This section briefly describes our formal setting and presents an overview of our approach. An introduction to Topological Data Analysis can be found in [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF].

Background

The input to our algorithm is a piecewise linear (PL) scalar field f : M → R defined on a simply-connected PL d-manifold M . Without loss of generality, we will assume that d = 3 (tetrahedral meshes) in most of our discussion, although our algorithm, as it extends Carr's, supports manifolds of arbitrary dimension. The scalar field f is provided on the vertices of M and is linearly interpolated on the simplices of higher dimension. We will additionally require that the restriction of f to the vertices of M is injective (which can be easily enforced with a mechanism inspired by simulation of simplicity [START_REF] Edelsbrunner | Simulation of simplicity: a technique to cope with degenerate cases in geometric algorithms[END_REF]). A level-set is defined as the pre-image of an isovalue i ∈ R onto 1(a)). Each connected component of a level-set is called a contour. In Fig. 1(b), each contour of the level-set of Fig. 1(a) is shown with a distinct color. Let f -1 ( f (p)) p be the contour that contains the point p. The Reeb graph [START_REF] Reeb | Sur les points singuliers d'une forme de Pfaff complètement intégrable ou d'une fonction numérique[END_REF] is a 1-dimensional simplicial complex (Fig. 1(b)) that is defined as the quotient space R( f ) = M / ∼ by the equivalence relation p 1 ∼ p 2 : [START_REF] Tierny | Loop surgery for volumetric meshes: Reeb graphs reduced to contour trees[END_REF] (i.e. a continuous map that retracts each contour to a point such that the restriction of such a map to its image is the identity) and where ψ : R( f ) → R is a continuous function that maps points of R( f ) to their f values. Note that the pre-image by φ of R( f ) induces a complete partition of M . In particular, the pre-image φ -1 (Σ 1 ) of a 1-simplex Σ 1 ∈ R( f ) is guaranteed by construction to be connected. This latter property is at the basis of the usage of the Reeb graph in visualization as a data segmentation tool (Fig. 1(b)) for feature extraction or isosurface traversal acceleration. In practice, φ -1 is represented explicitly by maintaining, for each 1-simplex Σ ∈ R( f ) (i.e. for each arc), the list of regular vertices of M that retracts to Σ. Moreover, since the Reeb graph is a simplicial complex, persistent homology concepts [START_REF] Edelsbrunner | Topological persistence and simplification[END_REF] can be readily applied to it by considering a filtration based on ψ. Intuitively, this progressively simplifies R( f ), by iteratively removing its shortest arcs, as described in further details in [START_REF] Pascucci | Robust on-line computation of Reeb graphs: simplicity and speed[END_REF]. This yields hierarchies of Reeb graphs that are accompanied by hierarchies of data segmentations, that the user can interactively explore in practice (see Fig. 1(c)). As discussed by Cole-McLaughlin et al. [START_REF] Cole-Mclaughlin | Loops in Reeb graphs of 2-manifolds[END_REF], the construction of the Reeb graph can lead to the removal of 1-cycles, but not to the creation of new ones. This means that the Reeb graphs of PL scalar fields defined on simply-connected domains are loop-free. Such a Reeb graph is called a contour tree and we will note it C ( f ). Contour trees can be computed efficiently by considering the combinations of two trees called join and split trees [START_REF] Carr | Computing contour trees in all dimensions[END_REF]. Given an isovalue

M through f : f -1 (i) = {p ∈ M | f (p) = i} (Fig.
f (p 1 ) = f (p 2 ) p 2 ∈ f -1 ( f (p 1 )) p 1 Note that f can be decomposed into f = ψ • φ where φ : M → R( f ) is a contour retraction
i ∈ R, the sub-level set L -(i) is defined as the pre-image of the open interval (-∞, i] onto M through f : L -(i) = {p ∈ M | f (p) ≤ i}.
The join tree, noted J ( f ), is a 1-dimensional simplicial complex obtained by contracting each connected component of a sub-level set to a point (similarly to the Reeb graph). The 0-simplices of J ( f ) that are attached to only one 1-simplex are either the root of this tree (corresponding to the global maximum of f ) or its leaves (corresponding to the local minima of f ). All other 0-simplicies (which induce branching in the tree) correspond to saddles of f that join distinct connected components of L -(i). The split tree S ( f ) of f is defined symmetrically, by considering the sur-level sets L + (i):

L + (i) = {p ∈ M | f (p) ≥ i}.

Overview

Our approach is based on a range-driven partitioning strategy, as illustrated in Fig. 2. First, given n t threads, the image of the domain f (M ) is divided into n t /2 contiguous, non-overlapping intervals I i that contain (nearly) the same amount of vertices of M (Sec. 3.1):

f (M ) = I 0 ∪ I 1 ∪ • • • ∪ I n t -1 (1) |σ 0 | i ≈ |σ 0 | j ∀i = j
where |σ 0 | i refers to the number of vertices of M mapping to I i . Next, two threads are assigned to each partition P i , P i being the pre-image of the corresponding interval, P i = f -1 (I i ) (Sec. 3.1). The two threads then compute the augmented contour tree of the restriction of the function to its partition (Sec. 3.2), with a variant of the algorithm by Carr et al. [START_REF] Carr | Computing contour trees in all dimensions[END_REF]: one thread builds the join tree, and the other the split tree 2 . This yields a forest of contour trees:

{C ( f ) 0 , C ( f ) 1 , . . . C ( f ) n t -1 }.
Finally, the output contour tree is retrieved by connecting the trees of the forest along common connected components of partition boundaries (Sec. 3.3).

Despite its simplicity, our range-driven approach exhibits many advantages. In particular, our strategy enables the computation of augmented contour trees, since it extends Carr's algorithm [START_REF] Carr | Computing contour trees in all dimensions[END_REF]. Second, since the input mesh is split into partitions of roughly equal size (in terms of vertices), the work load should be well balanced between the threads. Third, since it is range-driven, our approach allows for a full computation of the local contour tree within each partition (join and split tree computations, plus their combination in the contour tree) while previous approaches systematically delayed the combination to a post-process pass (implemented in serial). Finally , since it is range-based, our approach allows for a simple stitching of the local trees of the forest into the output contour tree, while previous approaches needed to run a special procedure on the common boundary of merged partitions [START_REF] Pascucci | Parallel computation of the topology of level sets[END_REF][START_REF] Maadasamy | A hybrid parallel algorithm for computing and tracking level set topology[END_REF].

ALGORITHM

This section details the algorithms for each step of our approach.

Domain partitioning

The first step of our approach consists in sorting the vertices of M by increasing function value, which can be efficiently done in parallel [START_REF] Tsigas | A simple, fast parallel implementation of quicksort and its performance evaluation on sun enterprise 10000[END_REF][START_REF] Singler | The Multi-Core Standard Template Library[END_REF][START_REF] Gnu | C++ Standard Library, Parallel Mode[END_REF]. This step can be done in

O(|σ 0 | × log(|σ 0 |))
where |σ 0 | is the number of vertices in M . Nex, the sorted list of vertices is split into n t /2 contiguous sets P i of roughly equal size, whose images correspond to the intervals I i described in Eq. 1. Next, each vertex set P i is extended into a set P i with the following procedure. Let f i -and f i + be the two extremities of the interval

I i : I i = ( f i -, f i + ).
The level-sets for the isovalues f i -and f i + are called interface level-sets. Let (σ 1 ) i -and (σ 1 ) i + be the the set of edges of M whose image contains f i -and f i + respectively. The vertex set P i is extended into P i by adding the vertices of (σ 1 ) i - and (σ 1 ) i + (red circles, Fig. 3). We call such vertices boundary vertices. Note that with this approach, two adjacent partitions P i and P j will overlap, precisely along the simplices crossed by f i -or f i + (triangles with red edges, Fig. 3).

This strategy guarantees that each connected component of an interface level-set is captured by the overlaps in between the partitions (triangles with red edges, Fig. 3). Therefore, all possible contours living in the interval I i are completely captured by P i . This guarantees that the restriction of the local contour tree C ( f ) i (computed on P i ) to the interval I i (in green in Fig. 3(a) and blue in Fig. 3(b)) is equal the restriction of the output contour tree C ( f ) to I i . This property will be of paramount importance to guarantee an efficient stitching of the contour forest into the output contour tree (Sec. 3.3).

In practice, this expansion procedure is performed efficiently by visiting in parallel all the edges (σ 1 ) of M and tracking the vertices of the edges crossing f i -and f i + for a given interval I i , in O(|σ 1 |) steps. In particular, each of the n t threads maintains its own list of boundary vertices, which are merged globally (and sequentially in practice since this merge implies minor computation times).

Local computations

The contour tree C ( f ) i of each of the n t /2 partitions P i is computed by two distinct threads. Note that in practice, the partitions P i are not copied, but represented implicitly. In particular, the list of vertices of the initial partition P i is represented by an interval in the global sorted list of vertices. The boundary vertices added in the expansion procedure described in Sec. 3.1 (red circles, Fig. 3) are represented by two sorted lists of vertices B i -and B i + , representing the boundary vertices below f i -and above f i + respectively.

Given a partition P i , its augmented contour tree is computed with a variant of the algorithm by Carr et al. [START_REF] Carr | Computing contour trees in all dimensions[END_REF], for which we briefly sketch the main steps here for completeness. This algorithm constructs the join tree J ( f ) (and symmetrically the split tree S ( f )) by visiting the vertices of the domain in increasing function value and by keeping track of the connectivity evolution of the sub-level sets L -(i) with a Union-Find data-structure [START_REF] Cormen | Introduction to Algorithms[END_REF]. Given an input sorted list of vertices, this steps requires

O(|σ i | × α(|σ i |))
steps, where |σ i | is the number of simplices in the partition P i and where α is an extremely slow-growing function (inverse of the Ackermann function). Next, the join tree J ( f ) and the split tree S ( f ) are combined into the output augmented contour tree C ( f ) by iteratively removing arcs connected to leaves either in J ( f ) or S ( f ) and adding them to C ( f ) [START_REF] Carr | Computing contour trees in all dimensions[END_REF]. This final steps takes a linear time with the number of arcs in C ( f ) (which is of the same order of magnitude than |σ i | in practice).

Our approach to the local computation of the augmented contour tree C ( f ) i for each partition P i only requires a slight modification to this algorithm. In particular, when constructing the join tree J ( f ) i , our algorithm first visits the boundary vertices B i -(if any) in increasing order. Next, it visits the vertices of P i by traversing the global sorted vertex list within the interval prescribed by the domain partitioning step (Sec. 3.1). Finally, it completes the traversal by considering the vertices of B i + (if any) in increasing order. For each of these three traversals, the join tree construction algorithm by Carr et al. [START_REF] Carr | Computing contour trees in all dimensions[END_REF] is applied as-is by the corresponding thread. The split tree S ( f ) i is constructed with a symmetrical pass by the other thread: B i + , then P i and finally B i -. Once the join and split trees are constructed, they are combined into the augmented contour tree C ( f ) i with the original algorithm [START_REF] Carr | Computing contour trees in all dimensions[END_REF] by one of the two threads. This combination does not require parallelization within each partition since its computation time is not significant, and since parallelization already applies among partitions.

Contour forest stitching

Once the n t threads have finished the computation of each local augmented contour tree C ( f ) i , the resulting forest is stitched into the final augmented contour tree C ( f ) with the following procedure.

During the local computation of the contour tree C ( f ) i (Sec. 3.2), each 1-simplex (each arc) that crosses an interface level set is added to a list of crossing arcs, noted X i . This corresponds in Fig. 3 to the arcs (a, j) in Fig. 3(a) and (hn) in Fig. 3(b) Then, the stitching procedure consists in visiting sequentially the list of crossing arcs X i for each local contour tree C ( f ) i . Given such an arc a i , its regular vertex v exactly above f i + (or below f i -) is identified through dichotomic search (vertex j in Fig. 3(a)). Since augmented contour trees store the destination of each vertex into the tree, it is possible to retrieve in constant time the homologous arc a j from the adjacent tree C ( f ) j which contains v. This corresponds to the arc (hn) in Fig. 3(b). Finally, a i is updated to form the union of the arcs a i and a j . This operation includes the modification of the higher extremity of a i (to use a j 's instead) as well as the concatenation of the two sorted lists of regular vertices (see Fig. 3(c)). Note that the vertex v can belong to multiple partitions. In such a case, a i will be updated iteratively to form the union of multiple arcs (a i , a j , a k , etc.), by successively applying this pairwise stitching in increasing order of function value (i.e. a i and a j will first be stitched, then the result of this stitching will be stitched with a k and so on). As discussed in Sec. 4, this final stitching procedure is extremely fast in practice (hence performed sequantially), since only a small portions of the arcs are visited (only the crossing arcs) and since the merging operation is simple (it simply consists in stitching pairs of arcs across interface level-sets). Note that the simplicity of this stitching procedure is due to our domain partitioning strategy, which guarantees that the restriction of a local tree C ( f ) i to the interval I i is equal to the restriction of C ( f ) to I i (Sec. 3.1). Note that the zipping procedure employed in the streaming Reeb graph computation algorithm [START_REF] Pascucci | Robust on-line computation of Reeb graphs: simplicity and speed[END_REF] could also be employed for the stitching of the local trees. However, this procedure admits a quadratic time complexity in the number of nodes of C ( f ), which is prohibitive in our approach, where only sub-quadratic routines have been used.

EXPERIMENTAL RESULTS

In this section, we present practical results obtained with a VTKbased C++ implementation of our algorithm (provided as additional material). Experiments were performed on a desktop computer with an Intel Xeon CPU E5-2630 v3 (2.4 GHz, 8 cores) with 64GB of RAM. All parallel tests are run with n t = 8 threads for n p = 4 partitions.

Detailed performance results

Table 1 first presents detailed performance results for various datasets. Plasma, Bucky and SF Earthquake are standard data-sets available at the AIM@SHAPE repository [START_REF]AIM@SHAPE Shape Repository[END_REF]: these are however too small to fully exploit our 8-core CPU. In fact, with our experimental setting, even the sort step is slower in parallel than in sequential for these data-sets (results not shown). This explains the low parallel speedups for such tests.

We thus focus in the following on larger tetrahedral meshes (upper part of Table 1) which have been obtained by triangulating regular grids. For the sake of comparison, these have systematically been upsampled to 256 3 vertices. One can first see that the additionnal Overlap step required to build the B i -and B i + lists (see Sec. 3.2) leads to low overheads. Moreover, the stitching step is efficiently performed in sequential which results in small run-times. As far as parallel speedups are concerned, the Elevation data-set is a synthetic and very simple one that shows good speedups, with a parallel efficiency of 5.38/8 = 67%. Moving to more complex data-sets (i.e. resulting in larger contour trees), one can see that we obtain good or average speedups (parallel efficiencies ranging between 55% and 40%), except for the Foot data-set which shows a limited speedup. We will detail these limitations and their causes in the next sub-section. The scalability of our approach is evaluated with Fig. 4, which presents the evolution of the speedup obtained by our algorithm as a function of the number of threads. The slope of these curves shows that the scalability of our approach, similarly to the speedups discussed above, is data-set dependent as well. In particular, our algorithm seems less scalable for the datasets which result in complex output trees (see the third column of Table 1 which shows the number of arcs per tree). Also, the smallest slope (nearly constant) is also observed with the Foot data-set, as further discussed in the next sub-section. Next, we compare our parallel implementation with the reference sequential libtourtre implementation in Table 2. Our speedups with respect to libtourtre range between 2.6 and 3.8 (except for Foot), which compares favorably (even if the data-sets differ) with the 1.6-2.8 speedups of Maadasamy et al. [START_REF] Maadasamy | A hybrid parallel algorithm for computing and tracking level set topology[END_REF] on tetrahedral meshes (see Sec. 1.1). Note additionally that in contrast to Maadasamy et al. [START_REF] Maadasamy | A hybrid parallel algorithm for computing and tracking level set topology[END_REF], our approach computes the augmented contour tree, which constitutes a more generic and versatile version of the contour tree.

We also added in Table 2 performance results of a naive parallel implementation of libtourtre [START_REF] Dillard | libtourtre: A contour tree library[END_REF] which uses a parallel sort and two OpenMP threads to build independently the join and split trees. Our parallel algorithm outperforms this naive implementation in all our test cases with exception of the foot, which further stresses the efficiency of our approach.

Limitations

In this section, we detail the three factors that limit our parallel speedups in practice.

As presented in Table 3, we can see that the actual number of vertices per partition is always greater than the ideal one (obtained by dividing the total number of vertices by the number of partitions). This is due to the boundary vertices that have to be added to each partition, which imply redundant computations that directly impact the parallel speedups. In particular, there can be important variations in the size of the interface level sets (in terms of crossed edges, in red in Fig. 3) within a single data-set, as shown in Fig. 5. Therefore the size of the overlaps between the partitions (expressed as the number of vertices in the lists B i -and B i + , see Sec. 3.2) can also vary. This induces redundant computations of varying importance within a single data-set. One can also see larger imbalance in the number of vertices per partition for more complex data-sets. This adds load imbalance to the parallel computations which further decreases the speedups.

This load imbalance is worsened by the fact that, depending on its impact on the join and split tree constructions, each vertex of M does not require the same processing time in practice. This effect is shown in Table 4 which shows varying computation speeds among the different partitions of a given data-set. In particular, the more complex is the contour tree, the larger is the gap among the computation speeds. One could choose to use several partitions per thread, with dynamic load balancing, in order to minimize such load imbalance, but this would introduce even more redundant computations (because of the boundary vertices). That is why we choose to use n t /2 partitions for n t threads: this indeed minimizes the redundant computations, while fully exploiting the complete independency between the join and split tree computations.

These two factors (redundant computations and load imbalance) jointly explain our lower speedups with more complex data-sets (especially for Foot).

Finally, a third factor also limits our parallel efficiencies for any data-set. The contour tree computation requires on average very few operations with respect to the number of memory accesses. Its operational intensity [START_REF] Williams | Roofline: An insightful visual performance model for multicore architectures[END_REF] is therefore low which makes such an application memory-bound, like most graph traversal algorithms [START_REF] Asanovic | The landscape of parallel computing research: a view from berkeley[END_REF]. As shown in Table 4, giving the parallel computations of the join and split trees to a single thread leads to higher computations speeds. Hence speedups linear in the number of cores cannot be obtained for such memory-bound applications: the memory bandwidth of the processor can not cope with the memory requirement of the 8 threads at a time. This also justifies our choice not to rely on the 2-way SMT (Simultaneous MultiThreading) capability of our CPU, and to use only 1 thread (instead of 2) per physical CPU core.

APPLICATION: DATA SEGMENTATION

The purpose of this work is to improve interactivity in contour-tree based data analysis. Our algorithm enables the computation of augmented contour trees within interactive run-times (up to a dozen seconds, Table 1), even for large and complex data-sets. Our aug- mented output data-structure is readily usable for data segmentation tasks, especially when features of interest coincide with connected components of level-sets. Combined with persistent homology based topological simplification (see Sec. 2.1), this enables the interactive exploration of hierarchies of data segmentations. Such a segmentation capability is an important feature for quantitative analysis [START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF]. Fig. 6 shows an example of such a contour tree based data segmentation on the EthaneDiol data-set (electron density). In particular, the contour tree is shown with a 3D embedding (colored cylinders). In this example, the pre-images of the arcs attached to maxima (which correspond to atom centers in the electron density field) reveal the influence regions for each atom (Fig. 6(a)) or for each atom group if topological simplification is employed (Fig. 6(b)). Such a segmentation enables the quantitative analysis of these features (for instance volume measurement). Saddles of the electron density (red spheres, Fig. 6) are located in configurations at the boundary between multiple atom influence zones. These correspond to covalent bonds. Last, the minima (blue spheres, Fig. 6) correspond to noise induced by the boundary effects. This noise is easily removed with topological simplification (Fig. 6

(b)).

A second use case is shown with the Enzo data-set (Fig. 7) which represents matter density in universe expansion simulations. Regions of high density correspond to the galaxies forming the core structure of the data (often called cosmic web). These are surrounded by large zones, called voids, of low matter density. Originally the augmented contour tree is composed of hundreds of thousands of arcs (Fig. 7(a)), which challenges its interpretation. Here, the interactive topological simplification of our contour tree based data segmentations progressively reveals the core structures of the cosmic web (Fig. 7(b)). In particular, regions of high matter density correspond to pre-images of arcs attached to maxima (opaque surfaces) while voids correspond to pre-images of arcs attached to minima (transparent surfaces).

In both use cases, such data segmentations were provided by the augmented contour tree, for which our approach enables a computation within interactive times (at most a dozen seconds, Table 1). In comparison, the topological simplification takes typically less than a second (1.25 for the simplification of 436,726 persistence pairs on the Foot data-set) while the update of the segmentation is instantaneous.

CONCLUSION

In this paper, we have presented a fast, shared memory multithreaded algorithm, based on a range-driven partitioning strategy, for the complete parallel computation of augmented contour trees on tetrahedral meshes. We have also provided a lightweight VTKbased C++ reference implementation of our approach based on OpenMP multi-threading.

In practice, especially for complex and noisy contour trees arising from acquired data-sets, the parallel speedups of our implementation are limited by redundant computations and load imbalance among the partitions, as well as by the memory-bound nature of the contour tree computation, which we have demonstrated experimentally (Table 4). Nevertheless, on simulated data-sets we managed to obtain good speedups with respect to the reference sequential libtourtre [START_REF] Dillard | libtourtre: A contour tree library[END_REF] implementation: these speedups exceed the ones of Maadasamy et al. [START_REF] Maadasamy | A hybrid parallel algorithm for computing and tracking level set topology[END_REF] (although the data-sets differ). Also, our experiments demonstrate the superiority of our approach over a naive parallelization of libtourtre. Moreover, our parallel approach allows for the augmented contour tree computation which enables data-segmentation.

In future work, we plan to investigate processing of larger datasets that do not fit in memory by designing efficient out-of-core algorithms. Also, our range-driven partitioning strategy seems particularly conducive to a parallelization of the persistence-driven simplification of the trees. We will explore this direction in the future to further improve the interactivity of the user-guided exploration of persistence thresholds for interactive data segmentation. Moreover, we believe our range-driven partitioning scheme could be further improved, in particular by trying to minimize the number of boundary vertices, to reduce computation redundancy. In that perspective, the contour spectrum [START_REF] Bajaj | The contour spectrum[END_REF] could be a good candidate for the selection of optimal cutting isovalues, which would cross only few edges. However, the research of the optimal trade-off between the balancing of the partition size and the computation redundancy remains an open question. Last, we would like to investigate in the future the extension of our approach to distributed systems. However, such use cases are particularly appealing in simulation contexts where the data is already distributed at the time of the simulation. This means that the data partitioning scheme is likely to be imposed by the simulation code, which would challenge our approach, based on a range-driven partitioning.

Figure 1 :

 1 Figure 1: Topology driven hierarchical data segmentation. (a) Input scalar field f (color gradient), level-set (light green) and critical points (blue: minimum, white: saddle, green: maximum). (b) Contour tree of f and its corresponding segmentation (arcs and their pre-images by φ are shown with the same color). (c) Contour tree of f and its corresponding segmentation, simplified according to persistence.

Figure 2 :

 2 Figure 2: Algorithm overview on the height function f of a volume M with two threads. (a) Input scalar field f (color gradient) with its critical points (blue: min, white: saddle, green: max). The domain is split into two partitions P i and P j of roughly equal size corresponding to the pre-images of contiguous intervals I i and I j of f (M ). The interface level-set between such two partitions is shown in red. (b) The augmented contour trees C ( f ) i (top) and C ( f ) j (bottom) are constructed in parallel for each partition. These local trees can be easily and efficiently stitched together to form the output augmented contour tree (right).

Figure 3 :

 3 Figure 3: Domain partitioning for a 2D toy example (height function). (a) Partition P i (green) with its overlap simplices (red) and its augmented contour tree C ( f ) i . (b) Partition P j (blue) with its overlap simplices (red) and its augmented contour tree C ( f ) j . The common region between the two partitions is made of all the triangles containing red edges and being crossed by the interface level-set (black dashes). (c) The output, stitched, augmented contour tree C ( f ).

Figure 4 :

 4 Figure 4: Speedups obtained by our algorithm as a function of the number of threads (one curve per data set).

Figure 5 :

 5 Figure 5: Size difference between two interface level-sets on the Foot data-set. (a) Interface crossing 1,507,357 edges (blue). (b) Interface crossing 606,276 edges (green). This difference can lead to load imbalance between our partitions.

Figure 6 :

 6 Figure 6: Augmented contour-trees for the EthaneDiol data-set (electron density) without (a) and (b) with topological simplification. An isosurface of electron density is shown in transparent grey in both cases. Pre-images through φ of arcs attached to maxima of electron density (i.e. atoms) are shown in color. This segmentation extracts the regions of influence for each atom (a) or for the main atom groups in the molecule (b).

Figure 7 :

 7 Figure 7: Augmented contour-trees for the Enzo data-set without (left) and with (right) topological simplification. The interactive topological simplification of our contour tree based data-segmentation progressively reveals the core structures of the cosmic web (region of high matter density, opaque surfaces), surrounded by voids of low density (transparent surfaces).

Table 1 :

 1 Running time of the different steps of the algorithm (in seconds). |M | denotes the number of vertices in the data-set, and |C ( f )| A the number of arcs in the output contour tree. Overall corresponds to the complete application, including memory allocations, etc.

	Data-set	|M | |C ( f )| A Sequential Sort Overlap Local trees Stitching Overall Speedup
	Elevation	82,906,875	1	29.18 0.91	0.18	4.18	0.14	5.42	5.38
	EthaneDiol	82,906,875	29	33.09 0.67	0.33	6.64	0.14	7.81	4.37
	Combustion	82,906,875	3649	28.04 0.61	0.34	6.19	0.15	7.31	3.83
	Boat	82,906,875	3235	29.94 0.69	0.41	6.17	0.14	7.44	4.02
	Jet	82,906,875	4171	26.82 0.65	0.36	6.03	0.15	7.21	3.72
	Enzo	82,906,875	282800	39.63 0.74	1.50	9.48	0.66	12.40	3.20
	Foot	82,906,875	844463	18.09 0.49	0.99	7.12	1.10	9.72	1.86
	Plasma	1,310,720	2851	0.18 0.01	0.01	0.06	0.01	0.09	2
	Bucky	1,250,235	4377	0.11 0.01	0.01	0.05	0.01	0.08	1.38
	SF Earthquake	2,067,739	11887	0.19 0.01	0.02	0.09	0.02	0.13	1.46

Table 2 :

 2 Overall running time comparison (in seconds) between the sequential libTourtre implementation (sTourtre), a naive parallel implementation of libTourtre (pTourtre) and our approach.

	Data-set	sTourtre pTourtre Speedup wrt.	Ours	Speedup wrt.
				sTourtre		sTourtre pTourtre
	Elevation	20.63	10.07	2.04	5.42	3.81	2.64
	EthaneDiol	23.47	13.96	1.68	7.81	3.00	1.79
	Combustion	21.26	12.39	1.72	7.31	2.91	1.70
	Boat	23.26	12.52	1.85	7.44	3.13	1.68
	Jet	20.60	11.50	1.79	7.21	2.86	1.60
	Enzo	32.51	18.07	1.80 12.40	2.62	1.46
	Foot	13.52	8.40	1.60	9.72	1.39	0.86
	Plasma	0.08	0.08	1.00	0.09	0.89	0.89
	Bucky	0.07	0.06	1.16	0.08	0.88	0.75
	SF Earthquake	0.12	0.10	1.20	0.13	0.92	0.77

Table 3 :

 3 Partition sizes (in vertices).

	Data-set	ideal	min	max
	Elevation	4,194,304 4,259,840 4,325,376
	EthaneDiol	4,194,304 4,362,086 4,616,938
	Combustion 4,194,304 4,353,986 4,635,078
	Boat	4,194,304 4,418,409 4,791,092
	Jet	4,194,304 4,358,176 4,701,586
	Enzo	4,194,304 5,234,144 6,474,322
	Foot	4,194,304 4,499,572 6,044,708

Table 4 :

 4 Computation speeds (in vertices/second) for join or split tree computations with our parallel implementation, and with one single thread to perform all computations required by our parallel approach.

	Data-set	Parallel: min	max 1 thread: min	max
	Elevation	1,623,500 1,768,720	2,151,250 2,292,390
	EthaneDiol	963,962 1,108,170	1,470,960 1,804,410
	Combustion	1,029,050 1,190,160	1,688,080 2,006,210
	Boat	1,055,410 1,237,030	1,463,880 1,985,720
	Jet	1,065,720 1,256,730	1,754,240 2,094,010
	Enzo	860,937	933,616	1,166,540 1,366,070
	Foot	1,120,560 4,031,030	1,220,80 5,195,250

Pascucci and Cole-McLaughlin use the term "augmented" with a different meaning than in the current paper. Here, it refers to the augmentation of the arcs with regular vertices (as introduced by Carr et al.[START_REF] Carr | Computing contour trees in all dimensions[END_REF]) while they refer to an evaluation of the Betti numbers of the level sets.

Note that n t threads could have been used (one thread per partition), by building these trees sequentially. However, our experiments showed that it was less efficient than using two threads per partition.

ACKNOWLEDGEMENTS

This work is supported by the French National Association for Research and Technology (ANRT), in the framework of the LIP6 -Kitware SAS CIFRE partnership reference 2015/1039.


Matthieu Martin

Alexandre Barzyk

Eric Bertin

Philippe Peyla

Salima Rafai

Photofocusing: Light and flow of phototactic microswimmer suspension

We explore in this paper the phenomenon of photofocusing: a coupling between flow vorticity and biased swimming of microalgae toward a light source that produces a focusing of the microswimmer suspension. We combine experiments that investigate the stationary state of this phenomenon as well as the transition regime with analytical and numerical modeling. We show that the experimentally observed scalings on the width of the focalized region and the establishment length as a function of the flow velocity are well described by a simple theoretical model.

Microorganisms such as phytoplankton have been shown very often to present some patchiness in their spatial distribution [1,2]. These spatial inhomogeneities occur either at large scales as in the ocean or at smaller scales in lakes or ponds for instance. Although this patchiness might in some case be due to inhomogeneities of the ecosystem -gradients in temperature or nutrients, gravity, etc. -the coupling between flow and plankton's ability to swim is suspected to play a major role in the spatial localization of plankton. Indeed, patchiness has been shown to be more important in the case of swimming as compared to non-swimming plankton [3].

More generally, swimming plankton represents one among many systems that form what is called active matter [4], i.e., systems composed of units able to self-propel autonomously. These can be molecular motors, bacteria, algae, fish to only mention some living systems. Active matter attracts a lot of attention from scientists as it represents a new "state of (non-equilibrium) matter" to be investigated. In particular, an interesting feature to understand is how the coupling between flow and motility occurs [5,6]. This question has been addressed for instance in the case of gyrotactic algae in a flow [7,8]. Because of the eccentric location of their chloroplast, these algae are subject to a torque induced by gravity and consequently tend to swim upward. Gyrotactic algae have been shown to concentrate in downwelling regions of turbulent flows [8]; this is a direct consequence of the migration of gyrotactic algae towards low vorticity regions as shown by Kessler [7].

Here, we investigate an analogous phenomenon that appears when flow vorticity is coupled to phototaxisa biased swimming of algae toward a light source. We propose a detailed quantitative description of the experimental dynamics and stationary regimes of what we call photofocusing [9] coupled to an analytical and numerical treatment of the problem in a spirit similar to [10,11].

A. Experimental part

The green microalga Chlamydomonas Reinhardtii is a biflagellated photosynthetic cell of 10 micrometer diameter [12]. Cells are grown under a 14h/10h light/dark cycle at 22 o C. Cells are harvested in the middle of the exponential growth phase. This microalga propels itself in a break-stroke type swimming using its two front flagella [13]. This way of swimming has been shown to be well characterized by a persistent random walk [13,14] in absence of tropism. The typical time of persistence is of the order of a few seconds. When subjected to a light stimulus (green wavelength, i.e., around 510 nm), cells tend to swim toward the light source and perform a quasi ballistic motion.

Let us first quantify the useful characteristics of the phototaxism mechanism of the cells. Whereas in absence of light bias, microalgae perform a persistent random walk, they adopt a ballistic motion toward the direction of a light source. Their response time is related to the mean reorientation time in the random walk of the cells; indeed, these reorientations allow the cells to scan the space for the presence of a light source. This time can be quantified by measuring the time correlation function of direction over a whole population of cells as a function of the time [14] and is found to be of the order of 3s. Once microswimmers swim toward the light, trajectories still show a dispersion in their orientation. Figure 1 shows the probability distribution of angles ψ(θ -π) in absence of light source and when a light is switched on in the direction θ = π in a fluid at rest.

The cells phototaxism is then coupled to a Poiseuille flow. A squared section channel is made of PDMS (1 × 1 × 50 mm 3 ). A white led is used as a light source placed upstream (fig. 2). A syringe pump imposes a flow with flow rates Q ranging from 0.5 to 3 mm 3 .s -1 . Observation is made with a macroscope (Leica APOZ16) coupled to a CMOS camera. A low magnification objective is used that provides a field of view of size 1.5 mm ×20 mm; this insures the observation of a stationary state (contrary to our previous study [9]). The flowing channel is enclosed in an occulting box with two red filtered windows for visualization. This is to avoid any parasite light that would trigger phototaxis.

In the stationary state, images of the focused cells (fig. 3) are taken and analyzed to deduce the concentration profile along y direction. For this, Beer-Lambert law is used to convert grey levels of intensity into concentration values. Volume fractions in the experiments are below 0.4%. The cell density distribution along y direction is then measured as a function of flow rate and in presence of light upstream. An example of distribution is shown in fig. 4 and compared to the analytical and numerical distributions (the models are described in the following sections). We measure the exponential length L x over which a stationary density profile is reached. This length is plotted as a function of the flow velocity in fig. model presented in the following.

B. Analytical modeling

Dynamics of individual particles

We consider a two-dimensional Poiseuille flow, seeded with point-like microswimmers which tend to reorient their swimming direction towards the light source, situated up-stream, at randomly chosen times. We denote as x the direction of the flow, y the direction transverse to the flow, and 2L the width of the channel (-L ≤ y ≤ L).

The velocity profile of the flow is given by

u(y) = U 0 1 - y 2 L 2 (1)
and the strain rate is:

γ = du dy = - 2U 0 L 2 y. (2)
Each swimmer moves at a constant speed v 0 in its own swimming direction, characterized by an angle θ (fig. 2) which evolves in time due to two different contributions. First, due to the vorticity of the flow, the swimmer rotates at an angular velocity

ω = - η γ 2 = ηU 0 L 2 y (3)
where η is a partial drive coefficient (0 < η ≤ 1) resulting from the fact that the microswimmers partially resist against the vorticity of the flow [15]. Secondly, each swimmer reorients, at randomly chosen times, its direction of motion close to the direction of the light source (fig. 1). More precisely, a new angle θ is chosen from a probability distribution ψ(θ -π), where ψ(θ ) is a distribution centered around θ = 0, and where θ = π corresponds to the exact direction of the light source (fig. 1).

Reorientation times occur randomly with rate α, meaning that the time duration τ between two successive reorientations is exponentially distributed, p(τ ) = α e -ατ .

Statistical description

We assume that the concentration of swimmers is small enough so that hydrodynamic interactions between swimmers can be neglected. The key quantity to describe the statistics of swimmers is the probability f (r, θ, t) for a swimmer to be at position r = (x, y), with a direction θ, at time t. The evolution equation for f reads

∂ t f + u(y)e x +v 0 e(θ) •∇f + ηU 0 y L 2 ∂ θ f = -αf +αρ ψ(θ-π) (4) where ρ is the concentration of swimmers, ρ(r, t) = 2π 0 f (r, θ, t) dθ (5)
and e x is the unit vector along the x-direction (direction of the flow), while e(θ) is the unit vector in the direction θ, with e(θ = 0) = e x . Eq. ( 4) contains three different contributions to the evolution of the probability f (r, θ, t).

The second term on the l.h.s. of Eq. ( 4) corresponds to the advection of particles under the combined effect of self-propulsion and Poiseuille flow. The third term on the l.h.s. describes the rotation of the swimmers induced by the shear flow. Finally, the r.h.s. of Eq. ( 4) describes the reorientation dynamics, through which the current angle θ is instantaneously changed into a new angle drawn from the distribution ψ(θ -π). It is convenient to define dimensionless variables by taking the channel half-width L as the unit of length, and the inverse of the reorientation rate α as the unit of time:

x = x L , ỹ = y L , t = αt , f = L 2 f . (6)
In the following, we drop the tildes to lighten notations. The dimensionless evolution equation for f is then given by

∂ t f + b η (1-y 2 )e x +β e(θ) •∇f +by ∂ θ f = -f +ρ ψ(θ-π) (7) with -1 ≤ y ≤ 1 and β = v 0 αL , b = ηU 0 αL . (8)
In the experiment, one has v 0 ≈ 10 -1 mm.s -1 , α ≈ 0.33 s -1 , 2L = 1 mm and a shear rate of a few s -1 . With these values, one finds β ≈ 0.6 and the parameter b, which measures the ratio between the shear rate and the reorientation rate, is of the order of 5.

Considering the coefficient b/η of the x-derivative in Eq. ( 7) suggests that the density profile relaxes, from the entrance of the channel, to its x-independent shape over a typical length scale L x ≈ (b/η)L = U 0 /α. Figure 5 shows the experimentally determined rescaled focusing length in the x-direction and the numerical values as a function of the rescaled flow speed. Both are well described by a linear law with a prefactor of 0.5 which is consistent with the expected value of v 0 /αL = 0.6.

Analytical determination of the concentration profile

We now wish to determine an analytical approximation of the stationary concentration profile far from the entrance of the channel, where the profile becomes invariant along the direction of the channel. Starting from Eq. ( 4), we expand the distribution f (r, θ) in angular Fourier modes, yielding a hierarchy of coupled equations. Using a simple closure relation, this hierarchy can be truncated at second (nematic) order, allowing for a simple Gaussian solution of the density profile given by ρ(y) = ρ max e -y 2 /2σ 2 (9)

with σ = (1 -ψ 2 ) 2ηψ 1 v 0 U 0 (10)
The coefficients ψ k are the Fourier modes of the distribution ψ(θ) (fig. 1). Technical aspects of the calculations are reported in Appendix A. The expression (10) of the width σ confirms the scaling in v 0 /U 0 observed in experimental data.

C. Numerical simulations

We consider a swimmer moving with a velocity v 0 = v 0 e(θ) in a Poiseuille flow along the x-direction. Its velocity v in the lab frame is given by

     v x = v 0 cos θ + U 0 1 -(y/L) 2 v y = v 0 sin θ , ( 11
)
where θ is, as previously, the angle between the swimmer and the x-axis (if 0 < |θ| < π/2, the swimmer is oriented downstream, while if π/2 < |θ| < π, the swimmer is oriented upstream). The swimmer is constantly rotated by the vorticity γ (see Eq. 2) of the external flow:

θ(t) = t 0 γ(y(t )) dt . ( 12
)
The light being situated upstream, the angle is set to π at a frequency α, where α -1 is the persistence time which corresponds to a ballistic motion of a cell in a fluid at rest. The parameter η corresponds to the fact that a cell resists the flow rotation as observed experimentally (i.e. η < 1) [15]. Here, we find η = 0.25 as a single adjustable parameter when fitting both d ∞ and L x . In our simulations, we calculate the trajectories of 10 000 swimmers by integrating Eq.( 11), r(t) =

D. Discussion and conclusion

In this paper, we show that the photofocusing experimentally observed can be well described both by numerical and analytical models which neglect hydrodynamics interactions (a reasonable assumption for such dilute suspensions). The experiments described in this paper allowed us to investigate both the transition regime of the photofocusing phenomenon as well as the stationary distributions of cells in the channel. An analytical continuous model is shown to satisfactorily describe the experimental data and provides us with useful scaling laws. Moreover, numerical simulations confirm these findings. Interestingly, the consistency between experimental data and modeling backs up a hypothesis that we previously emitted concerning the resistance of the cells to vorticity and that is a key ingredient to understand the peculiar rheological behaviour of Chlamydomonas Reinhardtii suspensions [15].

FIG. 1 .

 1 FIG. 1. Probability distribution of angles ψ(θ -π) in absence of light source (filled symbols) and when a light is switched on in the direction θ = π (open symbols); the cells are swimming in a fluid at rest. A truncated Lorentzian of width 1.2 is found to adjust the experimental points.

  FIG. 2. Experimental setup and coordinates system.

  FIG. 3. Examples of typical experiments. In absence of light (a), microswimmers are dispersed uniformly along y direction; in presence of light upstream (b), the cells are focalized in the central region of the channel. Note that here only zoomed pictured are shown for the sake of clarity but the experiments have explored a much wider field of view: 1.5 mm x 20 mm.

FIG. 5 .

 5 FIG. 5. Normalized focusing length scale Lx as a function of the normalized flow speed U0/v0. Experimental data (dots) and numerical data (squares). Data are obtained by varying U0. The fit of the experimental data represents a linear scaling as suggested in the text Lx ≈ U0/α.

t

  0 v(t )dt with different initial conditions: -L < y < L and x = 0. Then, we calculate the distribution of swimmers accross the channel for different x-values ρ(x, y) and evaluate the width of this distribution: d(x)/L = x)/L, we evaluate the length L x (fig.5) on which the focusing reaches its stationnary values d ∞ (fig.6).

ACKNOWLEDGMENTS

The authors thank A. Chardac for measuring the probability distribution of angles in the presence of light. This work has been partially supported by the LabEx Tec 21 (Investissements d'Avenir -grant agreement n • ANR-11-LABX-0030)

Appendix A: Analytical derivation of the density profile

We provide in this Appendix the technical aspects of the analytical derivation of the velocity profile given in Eq. ( 9). Let us define, for integer k, the angular Fourier mode as

Note that f -k = f * k , where f * k denotes the complex conjugate of f k , and that f 0 = ρ. Expanding Eq. ( 4) in Fourier modes then leads to

with ψ k the corresponding Fourier mode of the distribution ψ(θ),

(note that the symmetry ψ(-θ) = ψ(θ) has been taken into account). To shorten notations, we have introduced in Eq. (A2) the complex differential operators

For k = 0, Eq. (A2) reduces to the continuity equation

Note that v 0 f 1 is the exact analog, using the natural mapping between complex numbers and two-dimensional vectors, of the usual mass flux ρv, where v is the local collective velocity of swimmers. Hence Eq. ( A5) is equivalent to a standard continuity equation of the form

with the mass current ρv given by

where Re(z) is the real part of the complex number z, and Im(z) its imaginary part. We now search for the stationary density profile ρ(y), assumed to be invariant along the x-direction, that is the direction of the flow. Similarly, f 1 depends only on y. Under these assumptions, we find that in the stationary state, f 1 satisties ∂ y Imf 1 (y) = 0. Hence Imf 1 (y) is a constant, independent of y. Given that Imf 1 (y) is proportional to the flux of swimmers in the y-direction [see Eq. (A7)], Imf 1 has to be zero at the walls of the channel, implying Imf 1 = 0 for all y. The density profile is thus obtained by solving Imf 1 = 0, which requires to be able to express f 1 as a function of ρ and its derivatives. This can be done be taking into account the evolution equation of the mode f 1 , namely Eq. ( A2) for k = 1, which reads

From this equation, one can express f 1 as a function of ∂ y ρ and ∂ y f 2 ,

In order to close the equation, we need to express f 2 as a function of ρ and f 1 (and possibly their derivatives). Writing Eq. (A2) for k = 2, one finds

This equation also involves the higher order mode f 3 . To obtain a simple expression for f 2 , we thus need to make an approximation in order to close the set of equations. We first note that in the absence of flow (U 0 = 0) and selfpropulsion (v 0 = 0), the distribution f (θ) simply relaxes to f (θ) = ρψ(θ), yielding for the Fourier modes fk = ρψ k . The flow and self-propulsion can then be considered as driving mechanisms that perturb this distribution f (θ).

A simple closure relation is then to neglect the effect of these perturbations on f 2 (while, of course, keeping them on f 1 ). We thus simply assume in the following that f 2 = ρψ 2 . Under this approximation, the equation Imf 1 = 0 can be rewritten as, using Eq. (A9),

After integration, we obtain for the density profile

(where ρ max is a constant), as given in Eq. ( 9).