Guilhem Dupuy
email: guilhem.dupuy@univ-pau.fr

Bruno Jobard
email: bruno.jobard@univ-pau.fr

Sébastien Guillon
email: sebastien.guillon@total.com

Noomane Keskes

Dimitri Komatitsch
email: dimitri.komatitsch@univ-pau.fr

Bruno Liuppa

Jobard

INRIA-Magique3D Sebastion Liuppa

Guillon

Noomame Keskes
email: noomame.keskes@total.com

§ Total

¶ Migp

Isosurface extraction and interpretation on very large datasets in geophysics

Keywords: CR Categories: I.3.5 [Computer graphics]: Computational Geometry and Object Modeling F.2.2 [Analysis of algorithms and pronlem complexity]: Nonnumerical Algorithm and problems-Geometrical problems and computations I.4.10 [Image processing and computer vision]: Image Representation-Volumetric J.2 [Computer Applications]: Physical sciences and engineering D.1.3 [Programming techniques]: Concurrent Programming-Distributed programming Geometric computation and processing, Collaborative and distributed design, Geo-scientific applications

Introduction

Surface reconstruction for shape modeling is widely used in a large variety of domains (medicine, geophysics, ...). The marching cubes algorithm, introduced by Lorensen and al. [START_REF] Lorensen | Marching cubes: high resolution 3d surface construction algorithm[END_REF], is the most classical algorithm used for isosurface extraction. Due to the increasing size of processed datasets and extracted surfaces, a lot of improvements of the marching cubes have been proposed. Those improvements concern for instance ambiguities treatment [START_REF] Gelder | Topological considerations in isosurface generation[END_REF], reduction of the num-ber of traversed cells [Wilhelms andGelder 1992, Livnat et al. 1996], load balancing in parallel approaches [START_REF] Niguet | A load-balanced parallel implementation of the marching-cubes algorithm[END_REF]Nicod 1995, Mackerras 1992] and reduction of the number of generated triangles [START_REF] Shekhar | Octree-based decimation of marching cubes surfaces[END_REF].

In order to cope with the increasing size of datasets, isosurfaces might have to be simplified for reducing the number of generated triangles both for memory storage or more likely for visualization purposes. In a brute-force approach one would extract the full mesh and simplify it in a second pass. The problem with this method lies in the generation of a first heavy mesh which may not fit into the main memory before further simplification. In [START_REF] Attali | Extraction and simplification of iso-surfaces in tandem[END_REF]], Attali et al contribute in lowering memory requirements by introducing a tandem algorithm which combines isosurface extraction and simplification stages in one pass. Their solution reduces drastically the amount of vertices and triangles stored in memory during extraction allowing larger datasets to be processed.

The same way Attali et al addressed the memory problem raised by larger datasets, we propose in this paper a parallel, bloc-wise extended version of the tandem algorithm to accelerate the computation of simplified isosurfaces.

With our approach, the dataset is split in blocs and sent on computing nodes. On each node a local isosurface is extracted and semisimplified with a slightly modified tandem algorithm. Then nodes might receive an adjacent semi-simplified isosurface that will be merged with the local one. This merge operation ends with a simplification stage with relaxed edge constraints at their common interface that remove seams between them. The algorithm finishes when all local semi-isosurfaces have been merged and simplified. Our splitting/merging strategy ensures adjacent nodes to be processed together maintaining the memory budget as low as possible during the overall isosurface extraction.

We furthermore introduce an early components dumping mechanism that frees the memory as soon as independent surface components have been extracted, simplified and stored on disk along with meta-data for later high level exploration. This strategy has proved to be very useful on particular datasets such as in geophysics where the extracted features are numerous but small compared to the global size of the volume. The meta-data stored along these disconnected components can be used for filtering purpose during their visualization, for instance discarding those with too small volumes or keeping particularly aligned ones. This dumping mechanism can also be beneficial for noisy datasets or when a pertinent isovalue is not yet well determined and leads to many small disconnected components.

In the first part of this paper we describe the tandem algorithm. In the second part we propose a parallel extension of this algorithm with dumping of completed objects. The third part presents some computational experiments on our method.

The "Tandem Algorithm"

The main idea of Attali et al's algorithm is to alternate the extraction of a layer (see figure 1) and the simplification of the current overall extracted surface in order to reduce the amount of occupied memory. Indeed, a global simplification after a complete extraction The tandem algorithm is then a simple loop which iterates over the cross-sections applying two operations:

• extract which adds a layer (vertices, edges and triangles obtained with the marching cube algorithm between two subsequent cross-sections) to the current triangulation.

• simplify which simplifies the current triangulation. The last added layer is not simplified so that the next layer can be appended.

The simplification stage consists in applying the classical edge collapse algorithm [START_REF] Garland | Surface simplification using quadric error metrics[END_REF]. Each edge collapse operation has a cost which measures the numerical error it would introduce in the triangulation. Edge candidates for edge collapse are kept in a priority queue Q ordered by their cost. To evaluate this cost Attali et al have revisited the quadratic metric of [START_REF] Garland | Surface simplification using quadric error metrics[END_REF]] and proposed a quadratic error which is a weighted sum of a shape measure criterion and a mesh isotropy criterion.

The shape measure is similar to the one used in the original edge contraction algorithm. It measures the deviation introduced by collapse operation between the new vertices and the original surface.

The shape measure of a point x is then defined by

hc(x) = 1 Wc X t∈Uc wtd 2 (x, Pt) = 1 Wc x T Hcx
where Uc is the patch defined by all the neighboring triangles of point c. Pt is the plane spanned by the triangle t, wt its area and Wc = P t∈Uc wt. Hc is a positive definite matrix. Edge contraction ab → c leads to Hc = Ha + H b and Wc = Wa + W b .

The mesh isotropy criterion is introduced in order to counteract the creation of long and skinny triangles. Considering S ab the set of triangles containing the vertices a, b or both in the current triangulation, the mesh isotropy criterion for the point c is defined as the square distance of the point to the patch :

gc(x) = X t∈S ab wt(||x -t|| 2 + avg(t)) = x T Gcx
where t is the barycenter of the triangle t and avg(t) = 1 12 (||p|| 2 + ||q|| 2 + ||r|| 2 ) with p, q, r the vectors from t to the vertices of t.

The term gc is normalize by W = 3area(S ab )W 1/2 c /E0 in order to balance its influence with hc in the global cost defined by :

εα(c) = r c T [(1 -α) Hc Wc + α Gc W ]c
α is called the isotropy parameter, and represents a compromise between the shape measure criterion and the anisotropy measure criterion. In practice α is set to 0.4 for a good compromise between the two criteria.

For more details about those mathematical formulations, please refers to the original paper [START_REF] Attali | Extraction and simplification of iso-surfaces in tandem[END_REF].

During edge contraction, the resulting vertex position c is obtained by minimizing the local error function εα and the final error value εα(c) is used to order the priority queue Q. In any case, a candidate cannot be accepted if its shape measure, ε0(c), exceeds a positive constant error threshold E0. The simplify function consists then in emptying the queue Q by applying consecutive edge collapses.

Just after the extraction stage, the simplification stage has to cope with the heavy edge constraints on the last extracted layer. To counteract the artifacts theses blocked edges would introduce, an innovating way of scheduling the edge collapse was proposed, called time lag. The main idea is to delay edge collapses near the advancing front.

The time lag is based on the rank of a vertex, equals to its coordinate along the extraction direction (for example z if the cross-sections are taken along the z-axis). The rank of front is the maximum rank that has been extracted. Considering, height(u) = rank(u) and rad(u) = 1 for new vertices introduced by extraction, the contraction of an edge ab → c leads to :

height(c) = (height(a) + height(b))/2 rad(c) = (||a -b|| + rad(a) + rad(b))/2 reach(c) = height(c) + rad(c)
The contraction of an edge is prevented a long as its reach value is greater or equal to the rank of the advancing front. As detailed in [START_REF] Attali | Extraction and simplification of iso-surfaces in tandem[END_REF], if a and b belong to the last extracted layer, height(c) = rank(front) and since rad(c) ¿ 0, the contraction of ab would be prevented. Similarly, if a vertex of ab lies in the front plane, reach(c) would be greater than the rank of the advancing front.

The blocked edges are kept in a priority queue W, ordered by reach value. The function delay adds all edges of the last layer k in W. An edge is moved from W to Q if its reach value is lesser than the rank value of the advancing front. The function activate, described in algorithm 1, moves edges from W to Q. The figure 2 illustrates the effect of the time lag near the advancing front.

Procedure activate(

k : integer) While (reach(top(W)) ¡ k) do add top(W) in Q; pop(W); done End
Algorithm 1: At the k-th layer, the activate function fills up the edge collapse candidates queue Q with previously delayed edges of W. The tandem algorithm can then be written as in algorithm 2. It takes E0 as a parameter, the maximum shape error allowed.

Extended Tandem algorithm

Tandem algorithm was written to work on large datasets but experimental results (test on a noisy dataset of size 1626x7028x2000) showed us the limit of the scalability of this algorithm. On this cube, each extraction on the advancing front generate 1 250 000 triangles (see figure 3). In this case, queues' updates and simplification steps are too memory consuming. We propose then an extension of the tandem algorithm to increase its scalability. This extension consists in dumping parts of the extracted surface and dispatching the extraction process on sub parts of the dataset. The dumping requires to detect the end of a components extraction, but as we use the time lag technique, a component could be entirely extracted but its simplification could have been prevented. So we need to detect the end of its extraction as well as the end of its simplification.

To formalize the finalization of a component, we define an active component γ as the set of the vertices defining its shape. We consider Γ the set of all the active components, Γ = {γi}. We define Vw the set of all the vertices defining edges in the queue W (all the edges prevented by the time lag technique). A component is then finalized if it has no more edges to collapse (γ has no more edges in W), this proposition could be written as :

γ is f inalized ⇔ γ ∩ Vw = ∅
We define the function save(γ) which migrate a component to disk and the function clear(γ) which delete γ in current triangulation.

Dumping function could be written as in the algorithm 3. 5). One exploration scenario could be based on volume filtering like in [START_REF] Pivot | Complex volumic seismic interpretation by new geobodies extraction strategy[END_REF] where Pivot et al propose a workflow for complex volume seismic interpretation. They propose to extract isosurfaces on seismic attributes and to delete automatically inconsistent small "bubbles" (disconnected components with very small volumes regarding to other components). Another useful way to separate components is an analysis based on depositional direction (azimuth direction).

Parallel Algorithm

A feature of the Marching cubes is its granularity, facets of a grid cell (composed by 8 values of the data set) are computed independently of other grid's cells. So we propose to split datasets in sub blocks, to extract simplified surfaces in these sub blocks using the extended tandem algorithm and to merge then in a final surface.

Splitting We propose a hierarchical scheme for surface extraction. As in [START_REF] Muller | Adaptive generation of surfaces in volume data[END_REF], dataset is split recursively in two parts perpendicularly to its longest direction as long as sub blocks size is greater than a given threshold (see figure 6). Isosurfaces extraction are done on leafs of the generated tree and recursively merge to recompose the global surface.

Figure 6: Binary partitioning

Extraction in leaf nodes are done using the tandem algorithm combined with the dumping approach. In order to avoid refinement dependencies between adjacent blocks and prevent the artifacts due to bias in shape during simplification, we extend the time lag notion. As for original time lag the radius of a point c (result of the collapse ab → c) is given by rad(c) = (||a -b|| + rad(a) + rad(b))/2 with rad(x) = 1 if x is an original vertex of the isosurface. The contraction of edge ab is prevented as long as the sphere with center on the middle of ab and radius rad(c) is not totally include in the block. We show on the figure 7 how the introduction of the time lag influence the refinement progression near borders of the blocks.

Merging In our splitting strategy (binary partitioning), each node p is split in two parts p1 and p2, called children of p written cp. Reciprocally is the father of p1 and p2 written respectively fp1 and fp2. p1 is the brother of p2 written bp2 = p1 and reciprocally bp1 = p2.

Intuitively, the optimal strategy for block merging should be to recursively merge brothers together. In practice components are not equally distributed in sub blocks : extraction time of two brothers could be drastically different, and merging them implies to wait for the slower one. To counteract this problem we allow to merge a block with its brother or with one of its brother's children if they have a common border.

We thus need to define a criterion determining if two blocks have a common border. We consider that a block p could be split by three different planes (according to its longest direction). We define the type of a node as the position of the node after the split of its father. If a node is split along x p1 and p2 are typed respectively left and right. Similarly, cut along y types children as up and down, and cut along z as front and bottom. We define tree of root p (Tp), the set of nodes containing p and its recursive children. In this tree, ancestry of a node n (An, p) is defined as the set of nodes containing n and its recursive fathers up to p :

An, p = ( n if n = p n ∪ A fn, p else
Considering a block split in p1 and p2, then a node p belonging to Tp1 has a common border with p2 if

∀ p ′ ∈ Ap, p1, type(p ′ ) = type(p2)
Then if p fulfills this requirement it could be merge with p2.

Merging two sub blocks requires first to identify common points extracted on each block. Due to numerical inaccuracy in vertices extraction, we cannot identify common points with a simple coordinates comparison. Identification is then done by using a property of the Marching cubes algorithm : as a cell edge could have only one or zero points generated by the algorithm, 2 points belonging to common cell edges are identical and then merged. Once this merging step done, one has to simplify pasted area with the same error criteria as for the surface extraction E0. Like during the extraction step, the contraction of an edge ab is prevented as long as the sphere with center on the middle of ab and radius rad(c) is not totally included in merged blocks. We show on the figure 7 the merge of surfaces and the effect of the time lag on the quality of result surfaces. Finally the last step of the merge consists in dumping finished components on the disk.

Master/Slaves Implementation

In our application the target parallel machine is a load balancing cluster managed by Platform LSF-HPC. The common use of this parallel machine is to decompose the process in independant tasks, the task scheduler (hear LSF-HPC) dispatch tasks on the distant grid. When a task is finished, processor used for this task is liberated. The task scheduler could realloc it for the process if resources are sufficient for other running processes.

We propose an architecture composed of workers and a worker manager. Workers are distant processes which have to extract, merge and dump surface. The worker manager is a process, distant or not, which assign tasks to workers. This architecture, and formalism, is adapted to our target parallel machine (number of allocated processor is time-varying) but could used with others parallel machines.

Workers could be written as an infinite loop which request tasks from the worker manager. Tasks are composed of a task type and a worker id (each worker as an unique id use to send messages). We distinguish four task types :

• EXTRACT : worker has to extract surface from a sub block of dataset. The extraction is made with the extended tandem algorithm. Sub block borders, excepted the ones which are common to dataset borders, would prevent the edge collapse at their proximity using time lag. The extracted surface is appended to its current surface.

• SEND : worker a worker id, referring to a target worker. Worker has to send its current mesh to the target worker would merge it with its current one. After this sent, current surface of the worker is emptied.

• MERGE : worker has to merge its current surface with a mesh sent by another worker. After the merge operation, the new current surface is simplified and finished components are dumped.

• FINISHED : worker has to stop. The time life of a worker is necessary equal to global process duration.

The algorithm 5 describes the pseudo code of the worker. Worker manager a loop which iterate as long as the global surface has not been totally extracted and recomposed. Worker manager has a global view of the process and can dynamically dispatch tasks to workers, which ask for tasks when they begin or when they have finished their last given task. In order to define the worker manager mechanism (algorithm 6), we define the three following operations : has block, has neighbor and has surfaces.

Function has block determines if there's some cubes which have not been extracted. If some cubes have not been extracted the worker manager send it to a worker. If current worker have a triangulation which could be merge, un-extracted block candidates must have a common border with blocks of the worker. If current worker have no triangulation, one of the un-extracted block is chosen.

Function has neighbor determines if a querying worker has a neighbor in the current working workers. there's some neighboring workers, querying worker has to send its surface to one of them. If no neighboring worker is found, the worker manager look for a sub block to extract.

Function has surface determines if a worker has some waiting surfaces sent by other workers. If there's some surfaces, the query worker has to merge them with it current surface. This operation has priority due to memory consuming minimization policy.

Procedure

Computational Experiments

In order to analyse performances of our algorithm, we focus on the memory consuming during the extraction and the quality of the generated mesh.

Memory consuming A critical point in Marching cubes is the amount of memory needed for the generated triangles. The figure 8 illustrates the comparison between brute-force approach, tandem algorithm and tandem algorithm with dumping. Triangles are counted right after the process of a layer k. We see that tandem reduce drastically the number of triangles generated comparing to the classical Approach. Our dumping strategy enhances the tandem algorithm by clearly reducing number of triangles. Mesh quality Many applications need "well-shaped" triangles instead of long and skinny ones. A well-known evaluation of the aspect ratio of a triangle t = abc is ρ(t) = p λ2/λ1 where λ1 ≥ λ2 are the greatest eigenvalues of the inertia matrix of t. This matrix is defined by :

Mt = 1 3 [(a -t)(a -t) T + (b -t)(b -t) T + (c -t)(c -t) T ]
with t the centroid of the triangle t. ρ(t) = 1 if the triangle t is equilateral (λ1 = λ2) and ρ(t) = 0 if it is degenerate (λ2 = 0). As a global measure of the triangulation K (composed by n triangles), Attali et al propose a formulation of the anisotropy :

anisotropy(K) = 1 - 1 n X t ρ(t)
with 0 ≤ anisotropy(K) ≤ 1 and anisotropy(K) = 0 if all are equilateral.

The figure 9 shows anisotropy variation in function of the isotropy parameter α for the three approaches : sequential tandem algorithm, parallel algorithm with and without time lag. We clearly see effect of the time lag on the quality of the generated meshes. Parallel extraction generates a mesh with an equivalent quality as the sequential approach. The increasing of anisotropy (α greater than 0.9) is due to the fact that the shape criterion is no more took into account in the choice of candidate for edge collapse. Then the collapse operation is only based on triangle shape quality, and proposed points far from the original surface : most of those candidate are rejected by the error threshold E0. In practice, cost function must always take care of the shape measure and then prevent these degenerated cases : a value of α equals to 0.4 is a good compromise.

Application to geophysics Seismic cubes are very common data used on hydrocarbon exploration and analysis of those data is very helpful to specify depositional environment. A common tool for seismic analysis is based on texture analysis and seismic attributes computation. A seismic attributes highlight specific features of the image. For example figure 10 shows an example of attribute which detect all the chaotic and high amplitude area of the seismic As those facies could be related to the presence of hydrocarbon, making an inventory of all these area is very helpful for geophysicists. In order to make this inventory, we apply the following workflow : attribute computation, surface extraction and disconnected components sorting. As the attribute highlight the interesting part, then by adjusting a threshold we extract all the 3D shapes with our proposed approach (see result on figure 10). This workflow is applied to 1626x7028x2000 dataset (40 Go). The first interest of our approach is its parallelism : according to the number of available processor, we reduce the extraction time (the extraction have been performed in 5 minutes with processors instead of more than 3 hours for the tandem). The figure 11 shows the all extracted surfaces, and can see clearly a lot of small extracted objects : these objects could be due to noise on the attribute or to very small objects. In order to make a ranking of all these surfaces, we can easily sort them according to geometric criteria. Figure 12 illustrates this sorting by elimination of small objects. Only 328 components are now selected instead of 18 111. This workflow is very helpful to locate all these area and to analyse their repartition and relation.

Conclusion

Due to the increasing size of the datasets, the literature on isosurfaces extraction has focused on approaches which propose extract simplified surfaces or to dump surfaces as a soup of triangles. But none of these approaches propose parallel extraction with simplification and dumping of disconnected components. Our component-based dumping approach induces a file organization allowing interactive exploration of the volume. This organization could be applied on all surface reconstruction algorithms which generate many disconnected components. Our parallel processing based on dataset splitting and time lag extension could benefit to simplification methods which split their surfaces into patches.

Figure 1 :

 1 Figure 1: The k-th layer (in red) is the set of vertices, edges and triangles extracted with the marching cube algorithm between the cross-sections k-1 and k.

Figure 2 :

 2 Figure 2: Two partial triangulations constructed (left) without and (right) with the time lag : edges' length get progressively longer as they get further to the front layer.

Figure 3 :

 3 Figure 3: Number of triangles generated for each layer of the (1626x7028x2000) dataset. The portion with no triangles extracted is due to a hole in the dataset.

3. 1

 1 Components DumpingClassical out of core approaches migrate mesh to disk during extraction. Generated files are then soups of triangles ordered by direction extraction. Isosurfaces extraction, especially in geophysics, generate many disconnected components. Figure4illustrate an usual repartition of geologic events : finished components are shown in light gray and components which are still growing or which are not completely simplified are shown in dark gray. This distribution of disconnected components allows us to dump completed surfaces as soon as their simplification is finished.

Figure 4 :

 4 Figure 4: Distribution of extracted components. Finished components are colored in light grey and growing ones are in dark grey. Black arrow indicates the direction used for extraction.

  Dumping function Dumping function is then introduced in the tandem algorithm (algorithm 4). Procedure tandem( E0 : float) For k from 1 to number of cross-sections -1 do extract(Tandem algorithm with dumping An intrinsic property of our component-based dumping is that our generated file is ordered by component. We can then access easily to a sub set of components by reading sub parts of generated file. To optimize access to components in this file (called raw data file) we generate an index file (called index component file). To enhance the exploration of extracted surface, meta-data (oriented bounding box, volume, number of vertices and facets, ...) are computed on each component stocked in the index file (see figure

Figure 5 :

 5 Figure 5: File organization. On the left, the index file referring to raw data file (on the right).

Figure 7 :

 7 Figure7: On the left, figures illustrate the effect of the time lag on refinement blocks' boundaries before simplification. On the right, surfaces have simplified. On the top row, just edges which have a least a vertex on a border have been preserved. On the bottom row, edges have been prevented using the time lag method. This figure shows that applying the time lag method on the borders of blocks leads to better quality triangulations.

Figure 8 :

 8 Figure 8: Evolving size of the triangulated surface for the three algorithms : classical algorithm (extraction of the entire surface and simplification of it), tandem algorithm and tandem algorithm with dumping of finished components. Computations are done on a cube of size 401x1051x2001.

Figure 9 :

 9 Figure 9: Anisotropy of the mesh variation in function of the isotropy parameter α.

Figure 10 :

 10 Figure 10: Seismic attribute. (a) Cross sectional view of a seismic dataset. (b) Cross sectional view of an attribute computed on the seismic dataset. (c) Superposition of the attribute (color palette is saturated at a given threshold) over the seismic. (d) Global view of the extracted surfaces on the seismic attribute.

Figure 11 :

 11 Figure 11: Example of surface generated on a dataset of size 1626x7028x2000. This surface contains 18 111 disconnected components.

Figure 12 :

 12 Figure 12: Components extracted on dataset of size 1626x7028x2000 (see figure 11) have been filtered according to their volumes. This surface contains 328 disconnected components.

end If done End Algorithm 6: Worker manager

	worker manager()
	While ( ¬finished) do
	id ← get task request id();
	If (has surface(id)) then
	send task(id, MERGE);
	else
	If (has neighbor(id) ) then
	send task(id, SEND, targetId);
	else
	If (has block() ) then
	send task(id, EXTRACT);
	else
	send task(id, FINISHED);
	end If
	end If

Acknowledgments

This work was supported by Total company.


Romain Pacanowski

Xavier Granier

Christophe Schlick
email: schlick@labri.fr

Poulin Pierre

Pierre Poulin
email: poulin@iro.umontreal.ca

Sketch and Paint-based Interfacefor Highlight Modeling

Keywords: Descriptors (according to ACM CCS): I, 3, 3 [Computer Graphics]: Line and Curve Generation

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Motivation

3D modeling is generally decomposed in three components: geometry, animation, and appearance (i.e., lighting behavior). Over the years, and with an increasing interest since the publication of Teddy [START_REF] Igarashi | Teddy: a sketching interface for 3d freeform design[END_REF], simple sketch-based interfaces have been proposed for geometry modeling and more recently, for animation (e.g., [START_REF] Thorne | doodles: an interface for sketching character motion[END_REF]). This research in intuitive modeling has only recently emerged for the design of lighting behavior (e.g., [START_REF] Okabe | Illumination brush: Interactive design of all-frequency lighting[END_REF][START_REF] Todo | Locally controllable stylized shading[END_REF]), but is still in its early stages. Over all possible lighting effects, some of the most noticeable are highlights: they play an important role in the final appearance of an object, providing the user with convincing materials. Unfortunately, a highlight is a complex 4D function that depends on both light and view directions. Therefore, in a creative and artistic context, edition and creation of plausible highlights are still challenging tasks. Currently, most systems rely on the selection over pre-defined shading models (such as Phong shading) and the modification of their intrinsic parameters: the user's choice and freedom is thus limited. Recent approaches extend the modeling freedom by using a painting approach: a user creates and edits highlights by painting (cf. Figure1a) the expected result on a sphere [START_REF] Pattanaik | Creating Physically Correct Bidirectional Reflectance Distribution Functions[END_REF]. This approach is still limited by the underlying analytical models.

In this paper, we present an intuitive and flexible system for highlight design. Our system allows the user to specify c The Eurographics Association 2008. interactively through sketching, painting, and manipulation of vectorial data, the highlight's shape and its color gradient. With our new modeling tools, the user can create a wide variety of appearances. Thanks to the simple underlying model, our approach can provide real-time feedback and interactive rendering. These results are due to our two main contributions:

• we provide new sketching and painting tools associated with vectorial gradient edition in order to define easily different highlight characteristics, such as their shape and color variations. • we present a new highlight model based on a curve as a global shape representation for a highlight, and a texture for its color and refined shape description.

After reviewing previous work with a focus on interactive shading edition, we present the main principles of our approach and its corresponding modeling tools. We also introduce the underlying highlight models. We then illustrate their user controls and the resulting appearance, and discuss the limitations of the current implementation.

Previous Work

In computer graphics, one of the common ways to represent surface appearance and thus highlights, is to use a Bidirectional Reflectance Distribution Function (BRDF). However, its specification relies on a non-intuitive choice by an artist of a pre-defined model (e.g., [START_REF] Phong | Illumination for computer generated pictures[END_REF][START_REF] Blinn | Models of light reflection for computer synthesized pictures[END_REF][START_REF] Ward | Measuring and modeling anisotropic reflection[END_REF][START_REF] Schlick | An Inexpensive BRDF Model for Physically-Based Rendering[END_REF]), and an adjustment of its parameters that can have a nonuniform perceptual behavior. For an improved selection of the expected appearance, Ngan et al. [START_REF] Durand | MATUSIK Imagedriven Navigation of Analytical BRDF Models[END_REF] introduced a perceptually uniform navigation in a space defined by the different models and their parameters. On the edition side and once a given model has been selected, Ben-Artzi et al. [START_REF] Ben-Artzi A | Real-time BRDF editing in complex lighting[END_REF] proposed to project it into a given basis and to factorize it in a set of 1D curves that can be edited directly, even under complex illumination. However, even if these solutions provide users with improved selection and edition tools of predefined highlights, they do not allow the creation of more freely designed ones.

To increase modeling freedom, Colbert et al. [START_REF] Pattanaik | Creating Physically Correct Bidirectional Reflectance Distribution Functions[END_REF] develop Poulin and Fournier's work [START_REF] Poulin | FOURNIER Painting Surface Characteristics[END_REF] by proposing a painting interface. In their tool, called BRDF-Shop, different "painting" operations are introduced for the design of highlights. The resulting highlight is approximated by a sum of Ward lobes [START_REF] Ward | Measuring and modeling anisotropic reflection[END_REF], but other lobe models can be used for this fitting process [START_REF] Durand | Experimental Analysis of BRDF Models[END_REF]. Their approach does not guarantee that the painted highlights correspond to a reasonable fit because an arbitrary large number of lobes can be required. An increasing number of lobes will thus result in reduced interactivity. In constrat, Edwards et by Neumann et al. [START_REF] Neumann | Reflectance Models with Fast Importance Sampling[END_REF], where the shape of a lobe is defined on the tangent plane of the surface. Unfortunately, these indirect controls can be non-intuitive.

By removing constraints on realism, more freedom is provided to users in highlight design in a non-photorealistic context. One of the first solutions, based also on a painting metaphor, is the Lit Sphere [START_REF] Martin | The Lit Sphere: A Model for Capturing NPR Shading from Art[END_REF] by Sloan et al.. In their approach, the painted appearance on a sphere, for given viewpoint and light direction, is used as a texture projected on 3D surface, taking into account the similarity of the configuration between the viewpoint and the normal. Therefore their approach is only possible for a fixed lighting direction. Recently, Okabe et al. [START_REF] Okabe | Illumination brush: Interactive design of all-frequency lighting[END_REF] directly painted the lighting on a 3D surface. From this input, a 3D environment map is constructed and used as light source, but no control on highlights can be provided. Anjyo et al. [START_REF] Wemler | Tweakable light and shade for cartoon animation[END_REF] tweaked the shape of a highlight in cartoon-like (cf. Figure 1b) shading. However they do not really edit the BRDF: they locally move, scale, split, and merge the light sources in order to obtain the requested shape which is the only editable parameter. Smoothness, glossiness, and color variations are not taken into account. Similarly, Todo et al. [START_REF] Todo | Locally controllable stylized shading[END_REF] remove or add some highlights in toon shading, using offset texture defined for light key-directions. Unfortunately, this solution is limited to stylized shading. Of all the previous solutions, the most closely related to sketching is the work of Pellacini and Lawrence [START_REF] Pellacini | AppWand: Editing Measured Materials using Appearance-Driven Optimiza[END_REF]: They use strokes to select different areas of a BTF to interactively modify materials. In our approach, we also want to provide users with a fully interactive edition, but for the design of highlights. For user friendliness, our solution is based on sketching and painting metaphors. For efficiency and interactivity reasons, we directly manipulate the highlight shape and colors. Overview. As illustrated in Figure 2, our approach is based on the direct edition of highlight features, displayed on a plane oriented perpendicularly to the light mirror direction (R). On this plane, the artist changes the highlight shape and color-gradient with 2D editing tools. Consequently, the highlight features are defined for a given light (L), similar to the solution of Todo et al. [START_REF] Todo | Locally controllable stylized shading[END_REF]. The set of highlight features associated with a light key-direction is called a lighting configuration. A typical workflow in our system is firstly to define the highlight characteristics with our 2D tools (first shape, then color and gradient, as shown in Figure 4), and secondly to observe the highlight behavior for different light or view directions. By default, our system replicates the same highlight features for all lighting configurations. However, the user can edit them for a defined light direction. For undefined ones, our system smoothly interpolates every highlight features to ensure a coherent highlight behavior. Furthermore, at any time our system provides realtime feedback to user's editing actions.

Sketch-based Interface for Highlight

To provide intuitive tools for highlight modeling, we rely on three different interaction approaches adapted for each modeling action. For simplicity of use, each interaction is associated with one specific screen area (cf. Figure 3). As illustrated in Figure 4, user can sketch the shape, paint the highlight color or edit its gradient with vectorial tools.

Sketching has been recognized as an efficient tool to define a global shape. Therefore, we rely on this approach to modify the highlight shape. Through sketched strokes in a specialized area (cf. Figure 3(c)), the user selects part of or the whole highlight shape (cf. (Figure 4(a)). Once selected, a new target curve is drawn. A fitting process then approximates the sketch (cf. Figure 4(b)) using the underlying representation of curves (defined in Section 4). This sketch-based definition of the shape also controls the highlight shininess: the smaller the shape is, the shinier the highlight is. To in- crease the shininess, users have simply to scale down the curve.

To define color variations and gradients, painting is the most well-suited approach. The user thus uses a set of brush and gradient tools to refine the highlight shape and to control its color behavior. Brushes allow the user to edit the highlight color while the filters let her adjust the highlight intensity. This is also done in specialized areas, one for the gradient manipulation (cf. Figure 3 For an interactive modeling system, real-time feedback is crucial. Thanks to our representation based on a unique lobe defined with a simple curve and a texture (see Section 4); each modification is directly viewed on a selected 3D model (cf. Figure 3(d)). Since the highlight depends both on the light and on the view directions we provide also the user a visual representation of the current light key-direction (cf. Figure 3(e)). Furthermore, the user can select a part of the 3D object to retrieve local parameters, such as the most important light configuration or the highlight color and shape. Although this is not a 3D painting interface, it greatly facilitates the user's work when it comes to modifying precisely parts of the highlight.

Highlight Representation

For each lighting configuration, a distance field defined by a curve represents the highlight shape and a texture represents the highlight colors and intensity.

Highlight Shape. We define the highlight shape as a curve, the size of which controls the highlight shininess. The key here is to use this curve as the outline for the whole highlight. Outside of this curve, the highlight intensity is null, whereas inside the intensity is modulated by a color texture. In order to fill the shape, our representation should c The Eurographics Association 2008. allow a simple and intuitive access to the color texture. Furthermore, we also need a representation that can be easily interpolated between light key-directions. This is required for the reconstruction of highlight shape for unspecified light directions.

Therefore, we use the spline-based polar curve proposed by Crespin et al. [START_REF] Crespin | Implicit sweep objects[END_REF]. Each control point for the curve is located by an angle θ is defined by its radius ρ, and its left and right tangents. These parameters (radius and tangents) are easily interpolated. Furthermore, this curve representation enables the definition of a distance field d(ρ, θ) that starts from the center of the Sketching Plane (cf. Figure 5(b)). This distance field is then used to fill the shape with a color texture.

Texture Definition. The color texture represents the highcolor and intensity inside the shape. Therefore, black texels indicate areas where the highlight intensity is null. This allows a refinement of the curve-defined shape. Furthermore, by controlling the intensity gradient, we control the highlight glossiness.

Two parameterizations (polar or Cartesian) are used for texture lookup. We found it easier to work with polar coordinates when the texture is used to define a simple color gradient, and Cartesian coordinates for more generic textures, such as those created using painting tools. When using the polar parameterization (cf. Figure 4), the horizontal axis of the texture color component represents the angular variation whereas the vertical axis represents the radial one.

Since the color texture is enclosed by the star-shape polar curve one, we can also use it to define more complex shapes. To further ease this process, we separate, as it is done in common painting software, the color texture into two multiplicative layers. One layer stores the highlight color and the other stores the highlight intensity as gray level (cf. Figures 9 and10).

Parameterization. Finally, we need to define the parameters ρ and θ as functions of the view and light directions:

ρ = V • R cos θ = V • X sin θ = V •Y (1)
where • denotes the dot product, R the light mirror direction, and V the view direction. The directions X and Y are defined as follows:

X = T × R Y = R × X (2)
where × denotes the cross product, and T the local geometry tangent (cf. Figure 5)(a)). With polar coordinates, to access the color texture, the system uses the following (u, v) texture

u = ρ 1 -d(ρ, θ) v = θ 2π . ( 3
)
To summarize, for a given light direction, the system smoothly interpolates a distance field d(ρ, θ) and a color texture that are both used and accessed, at rendering time, according to the current view direction.

Implementation. We only use the CPU to perform the least-square fitting of the polar curve. To achieve interactive feedback, we rely on current graphics hardware capabilities using OpenGL Shading Language.

The shape of the highlight is stored as a floating 3D texture. A slice of the 3D texture contains parameters (radius and tangents) of one curve's control point for each lighting configuration. The number of slices is the number of control points of the polar curve. Therefore, a slice parameterizes a spherical triangle representing all possible lighting configurations. We pack every color texture (one for each lighting configuration) into a single texture array (ARB_texture_array extension).

practice we use between 20 and 60 control points with 3 different lighting configurations, thus requiring less than 40 KB to store the highlight features (shape + color). Therefore, our representation is a very low-consumption memory solution and one may use it to store many different shaders on GPU resident memory.

Use-case Scenarios and Results

We perform all the presented experiments on an Intel Core 2 Duo T7500 CPU workstation with a GeForce 8600 M GT. For a rendering resolution of 1280 × 1024 and an object complexity of 35 000 triangles, we report a frame rate of 60 frames per second. On the interaction side, each action leads to an update of the interface in less than 100 ms. Through all presented results, keep in mind that the final highlight appearance results from the combination of diffuse and highlight colors. Figure 7 shows a scenario where the user sketches a trefoil-like shape with different sizes and vectorial gradients for each light key-direction (cf. Figure 7(d)). For the normal incidence direction, the color texture is a green gradient; at tangent (resp. bitangent) incidence, the color texture is a blue (resp. red) gradient. As demonstrated in Figures 7(b) and(c), when the light shifts toward the tangent (resp. bitangent), the highlight becomes more blue (resp. red). The yellow color in Figure 7(c) comes when interpolating between the red and green colors stored for two lighting key-directions. This explains why the upper part of the sphere, where light configurations for normal and bitangent key-directions have equal importance, is yellower while the lower part of it, where the configuration for bitangent direction is more important, appears redder. The user can sketch more complex shapes, such as the one presented in Figure 6. Contrary to the previous example, painting tools were used to create a more complex color texture. For both previous examples, our system uses polar coordinates for texture look-up.

With a similar approach, the user can create complex appearances, close to realistic ones, such as the dispersion-like effect of Figure 8. This is simply done by creating two different color textures (using Cartesian coordinates) and two different shapes (cf. Figure 8(d)). For the lighting configuration at normal incidence, the texture is a white gradient, where for grazing configurations, it is a rectangular rainbow gradient. Furthermore, a circular shape is used for the normal incident light key-direction and a rectangular shape for the others.

The two layers of the color texture can be used to refine the star shape created with the sketched curve (cf. Figures 9 and10). We model the highlights of Figure 9 with a smiley bitmap texture (cf. upper right corner of Figure 10). As explained earlier, points of the sketching plane falling outside the bitmap texture are not shaded (i.e., the highlight intensity is null). Since the color texture is enclosed in the polar curve shape, it can be deformed to achieve new effects, like the unhappy smiley of Figure 9(b). Finally, Figure 10 shows an example with more complex color texture. The color layer is filled according to the lighting configuration, whereas the intensity layer, which contains the shape, remains constant (cf. Figure 10(d)). Discussion. Our system and its highlight representation are only a first step and some limitations remain. First, since the painted and sketched highlights are projected from the sketching plane on the surface, which is not necessarily planar, some shape deformation can occur. However, our system allows changes and adjustments of the shape real time, such that the user may compensate for this deformation. This was done, for example, in Figures 9(a) and (b) to compensate the vase curvature. Even though our curve-based representation ensures a smooth interpolation between different shapes, some interpolation artefacts may occur when using complex color textures (as in Figure 10). This limitation is related to the partially (un)solved general problem on how to morph smoothly between two arbitrarily different images. Third, remember that our system models a highlight expressed as a function of the view and light directions. However, these directions are parameterized according to local Therefore, a complex surface with many normal variations exhibits many highlights but necessarily on a large area. This explains why highlights in Figures 8 and9 on the upper edge of the vase are almost indistinguishable.

Although it might still be complex to design a highlight from 2D sketches and paintings rather directly paint it the object, we think that this drawback is compensated by our interactive feedback and our visual hints. Finally, even though we show examples with only three different lighting configurations, one may use more lighting configurations (within hardware limits). However, we found, through empirical testing, that increasing the number of lighting configurations increases the highlight design complexity for the artist.

Conclusion and Future Work

We have presented a new system based on sketching and painting metaphors for highlight design. A highlight is represented with a polar curve for its shape and a texture for its color gradient and for precise shape enhancement. Our simple and compact representation, suitable for hardware rendering, allows real-time edition of highlights, which can be various and reach some realistic aspects. Our tools and representation also provide more editing and creative freedom than previous approaches. We believe that our solution is more intuitive for highlight control in shading design.

Our upcoming work follows two streams. First we would like to extend it as a full 3D painting and sketching interface where the user can draw directly the appearance for simple and complex illuminations. We believe that this would be possible with the same parameterization, but a new representation. Second, we want to develop new tools and a new highlight model with spatially varying properties to help the user control highlights and highly detailed geometry.

Figure 1 :

 1 Figure 1: Three different styles of highlight obtained with different systems. (a) BRDF-Shop Physic based highlight, (b) Anjyo [AWB06] cartoon based highlight, and (c) example of highlight obtained with our system.

Figure 2 :

 2 Figure 2: The Sketching Plane. The user sketches the highlight shape and paints the highlight color on a Sketching Plane oriented perpendicularly to the light mirror direction (R).

Figure 3 :

 3 Figure 3: User Interface Main View.

Figure 4 :

 4 Figure 4: (a) The external curve represents the user's selection sketch.(b) The cyan curve represents the user's sketch whereas the blue one represents the fitting result. (c) A vectorial color gradient defined by the user. (d) Sketching Plane view that displays the resulting highlight when applying (b) and (c).

  (b) and 4(c)), one for a precise edition on the final results (cf. Figure 3(a) and 4(d)). The combination of the shape and the gradient texture is explained in Section 4. If the embedded tools are not sufficient, an experimented artist can load an image created with any other image manipulation software.

Figure 5 :

 5 Figure 5: (a) The X and Y directions constitute the Sketching Plane local frame. They are computed according to the light mirror direction R and the geometry local tangent T . To compute the highlight's intensity, we project the view direction V on the Sketching Plane. (b) The distance field d(ρ, θ) by the polar curve parameterizes the Sketching Plane.

Figure 6 :

 6 Figure 6: A more complex highlight shape.

Figure 7 :Figure 8 :Figure 9 :Figure 10 :

 78910 Figure 7: Light key-directions illustration. The highlight color and size is set according to the light direction. When the light direction is collinear with the geometry normal (a), the highlight exhibits a green color, whereas when the light shifts toward the tangent (b) (resp. bitangent (c)), the highlight becomes more blue (resp. red). For configurations (b) and (c), the highlight size is larger than for (a) to increase the color-shift effect.

c The Eurographics Association 2008.


Christoph Weiss
email: weiss@theorie.physik.uni-oldenburg.de

Niklas Teichmann

Differences between mean-field dynamics and N-particle quantum dynamics as a signature of entanglement

Keywords: numbers: 03.75.Gg, 05.45.Mt, 74.50.+r Keywords: Chaos, mesoscopic entanglement, Bose-Einstein condensation

A Bose-Einstein condensate in a tilted double-well potential under the influence of time-periodic potential differences is investigated in the regime where the mean-field (Gross-Pitaevskii) dynamics become chaotic. For some parameters near stable regions, even averaging over several condensate oscillations does not remove the differences between mean-field and N -particle results. While introducing decoherence via piecewise deterministic processes reduces those differences, they are due to the emergence of mesoscopic entangled states in the chaotic regime.

Experimentally it is possible to generate precisely controllable double-well potentials for Bose-Einstein condensates (BECs) (Ref. [1] and references therein). A future goal for this system is the realization of mesoscopic entanglement [1]. When combined with a time-periodic potential difference between the two wells, a BEC in a double well could also be used to investigate quantum chaos [2,3,4,5]. Another system which is widely used to investigate quantum chaos is the quantum delta-kicked rotor [6,7,8]. Research on quantum chaos includes topics like quantum signatures of chaos [START_REF] Haake | Quantum Signatures of Chaos[END_REF], quasi-stationary distributions [START_REF] Breuer | [END_REF], entanglement [11,12] and decoherence [START_REF] Braun | Dissipative Quantum Chaos and Decoherence[END_REF].

Often, a mean-field approach within the Gross-Pitaevskii equation is applied to describe BECs. Still, there are noticeable differences between mean-field dynamics and quantum dynamics: only the latter displays the well known collapse and revival phenomenon (cf. [START_REF] Greiner | [END_REF]). By time-averaging over several of those oscillations, these differences usually disappear. However, preliminary results [15] for the periodically driven doublewell potential indicate that even under time-average, mean-field dynamics and quantum dynamics can display qualitatively different results in the regime for which the mean-field dynamics become chaotic.

In this Letter, these differences are investigated systematically. First, the N -particle Hamiltonian is introduced for which the Gross-Pitaevskii equation corresponds to a driven nonrigid pendulum. If decoherence is implemented on the N -particle level via piecewise deterministic processes, the quantum dynamics can become qualitatively similar to the mean-field dynamics. The reason for the remaining differences between both approaches is the emergence of mesoscopic entangled states.

To describe a BEC in a double well with single-particle tunneling frequency Ω and pair interaction energy 2 κ, we use the Hamiltonian in two-mode approximation [16]:

Ĥ = - Ω 2 â1 â † 2 + â † 1 â2 + κ â † 1 â † 1 â1 â1 + â † 2 â † 2 â2 â2 + µ 0 + µ 1 sin(ωt) â † 2 â2 -â † 1 â1 , (1)
where â( †) j creates (annihilates) a boson in well j; µ 0 models the tilt and µ 1 is the driving amplitude. Such Hamiltonians have been used for schemes of entanglement generation [17,18]; without the periodic driving, entanglement has been investigated in BECs [19,20]. Other applications include high precision measurements, many-body quantum coherence [21,22] and spin systems [23].

On the level of the Gross-Pitaevskii equation for the above model, a wave function is characterized by the variables θ and φ, where cos 2 [θ/2] (sin 2 [θ/2]) is the probability of finding the condensate in well 1 (well 2) and exp(iφ) is the phase between the two wells. The corresponding N -particle wave-function ("atomic coherent states" [START_REF] Mandel | Optical coherence and quantum optics[END_REF]) with all particles in this state reads (in an expansion in the Fock-basis |n, Nn with n atoms in well 1):

|θ, φ = N n=0 N n 1/2 cos n (θ/2) sin N -n (θ/2) × e i(N -n)φ |n, N -n . (2)
The mean-field dynamics can be mapped to that of a nonrigid pendulum [15,[START_REF] Smerzi | [END_REF]; including periodic driving the Hamilton function reads (z = cos 2 (θ/2)sin 2 (θ/2)):

H mf = N κ Ω z 2 -1 -z 2 cos(φ) -2z µ 0 Ω + µ 1 Ω sin ω Ω τ , τ = tΩ . (3)
The experimentally measurable [1] population imbalance z/2 can be used to characterize the mean-field dynamics. Fig. 1 shows typical Poincaré surfaces of section. The initial parameters were chosen such that tunneling in the driven, tilted double-well potential is enhanced by "photon"-assisted tunneling [26] (cf. Ref. [27]). If the interaction is not too low (N κ/Ω 0.4 . . . 0.6), regular and chaotic dynamics coexist (Fig. 1.a cf. [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF]), for low interaction the dynamics are regular (Figs. 1.b and 1.c).

For the parameters corresponding to the Poincaré surface of section in Fig. 1.a, Fig. 2.a displays the differences between N -particle and mean-field dynamics by numerically calculating (using the Shampine-Gordonroutine [START_REF] Shampine | Computer Solution of Ordinary Differential Equations[END_REF]) the time-average of the (experimentally measurable [1]) population imbalance J z /N (which corresponds to the mean-field z/2):

J z T N = 1 N T T 0 dt 1 2 ψ|â † 1 â1 -â † 2 â2 |ψ , (4)
where for J z /N = ±0.5 the entire condensate is in the left, respectively, right well. Each point represents an initial condition (2). The differences are small if the meanfield dynamics are regular (cf. Fig. 1.a) while they can be rather large in the chaotic regime (up to half the theoretical limit, max{|z/2 -J z /N |} = 1). Most of the deviations between N -particle dynamics and mean-field dynamics in Fig 2 .a lie within twice the root-mean-square (r.m.s.)-fluctuations of the N -particle dynamics. However, contrary to the preliminary results of Ref. [15], for many initial conditions in the (classically) chaotic regime the differences can be very small; they are large near the boundaries of stable regions. 2√ N ) ≤ 0.05, thus clearly indicating that more than one atomic coherent state is involved. Bose-Einstein condensates of N ≈ 100 have been realized experimentally [START_REF] Chuu | [END_REF], both the validity of the two-mode approximation will be better and life-times of mesoscopic entangled states will be longer than in larger condensates. However, even when the calculation is repeated for N = 1000 particles, the differences in the chaotic regime remain. As the (non-linear) Gross-Pitaevskii equation does not allow any superpositions, decoherence should reduce the differences between mean-field and quantum dynamics.

In this Letter, we use a piecewise deterministic process (PDP) (Ref. [START_REF] Breuer | The Theory of Open Quantum Systems[END_REF], cf. [START_REF] Dalibard | [END_REF]) to model decoherence. To avoid to have to introduce decoherence also on the mean-field level (the atomic coherent states (2) become orthogonal in the limit N → ∞), we use the projection on the atomic coherent states [START_REF] Mandel | Optical coherence and quantum optics[END_REF]:

1 = N + 1 4π dθ sin(θ) dφ |θ, φ θ, φ| . (5)
Now, the PDP simplifies to having jumps on one of the atomic coherent states (2) after time t with probability

p jump = 1 -exp(-αt) , α = const. > 0 , (6)
and Hamiltonian dynamics (1) between jumps. The state on which the wave-function is projected is determined by the probability distribution

p θ,φ dΩ = N + 1 4π | ψ|θ, φ | 2 sin(θ) dθ dφ . ( 7
)
FIG. 3: Time-averaged population imbalance Jz T /N for various driving amplitudes µ1 in a tilted driven double well (2µ0/Ω = 3.0, ω = 3Ω, T = 100/Ω). The BEC initially is in the lower well (z0 = 1). Solid line: Jz T /N for N = 1000 is a smooth curve as opposed to the mean-field results depicted in the inset, which display chaotic jumps for small changes of the driving amplitude. Dots in the main plot: If decoherence is included via the PDP-process described around Eq. ( 6) with on average ≃ 5 jumps (α = 1/20) [START_REF]A higher α leads to a larger statistical spread. α includes the influence of the environment; it would, e.g., be larger for larger BECs[END_REF], the behavior is closer to the mean-field dynamics. Many dots lie in the area defined by the curves ( Jz T ± ∆Jz T )/N (dashed lines).

Figure 3 shows that the PDP can qualitatively reproduce the results of the Gross-Pitaevskii equation [START_REF]A higher α leads to a larger statistical spread. α includes the influence of the environment; it would, e.g., be larger for larger BECs[END_REF]. Without introducing the decoherence, the qualitative difference between mean-field and quantum dynamics are quite large; averaging over several PDPs would again result in a smooth curve within the error-bars in Fig. 3. While BEC-research in quantum chaos often assumes the validity of the Gross-Pitaevskii equation [2,3,4,5,[START_REF] Martin | [END_REF], at least for the model investigated here, only decoherence can lead to the chaotic behavior predicted by mean field.

Furthermore, differences between quantum dynamics and mean-field dynamics can also occur in the regular regime: Fig. 4 shows that, at least for N = 100, the differences can even lie above the result for many initial conditions in the chaotic regime (Fig. 2.a). One way to reduce the differences is to average over the Husimi distribution (7) (see Fig. 4, cf. Refs. [2,35] and references therein). This decreases the peaks of the differences between mean-field and quantum dynamics by a factor of 2 (in the chaotic regime, the factor can be of the order of 5). A perfect agreement cannot be expected as the averaged probability distribution on the mean-field level is always added whereas in quantum mechanics also destructive interference can occur.

On the level of quantum dynamics, the differences could be due to either a distribution of many atomic coherent states -or maybe even mesoscopic superpositions. For our model all mesoscopic quantum superpositions of all N particles being either in one quantum state or in another can be expressed as a sum of two atomic coherent states (see the explanation before Eq. ( 2)):

|ψ sp = η |θ 1 , φ 1 + e iγ |θ 2 , φ 2 , 0 ≤ γ ≤ 2π (8)
If both parts hardly overlap, | θ 1 , φ 1 |θ 2 , φ 2 | ≪ 1, the normalization η ≃ 1/ √ 2 and |ψ sp is a highly entangled mesoscopic state (for finite N , the only two orthogonal atomic coherent states (2) are |0, φ 1 and |π, φ 2 ). In a two-dimensional projection (cf.

p fid = 0 : R1 and R2 overlap 1 2 √ m 1 + √ m 2 2 : else . (9)
Yet, this only indicates entanglement if p fid > 0.5. With

σ ent = m 2 m 1 p fid , σ ent ≤ p fid (10)
even values of σ ent 0.5 can identify mesoscopic superpositions (Fig. 5.c). In Fig. 5.a, the maximum value of entanglement (evaluated at τ = 5 and 10) is plotted as a function of the initial condition (z 0 , φ 0 ): within the chaotic regime (left), entanglement generation happens on faster time-scales than in the regular regime (right); for longer time-scales (Fig. 5.b) the entanglement in the entire chaotic regime is more pronounced. It reaches particularly high values near initial conditions with large differences in the time-averaged population imbalances (Fig. 2.a). We obtained qualitatively similar results also for other values of driving amplitude and interaction.

To conclude, generation of mesoscopic entangled states can be a signature of quantum chaos for a BEC in a periodically driven double well potential. We investigated the driving near multi-"photon" tunneling resonances [26] which were recently observed experimentally for a BEC in an optical lattice [36]. While decoherence can lead to a "chaotic" behavior similar to the predictions of the Gross-Pitaevskii equation, the differences between quantum dynamics and mean-field dynamics are due to the emergence of mesoscopic superpositions. If the meanfield dynamics are chaotic, the entanglement generation is accelerated and its values are enhanced.

FIG. 1 :

 1 FIG.1: Poincaré surface of section for the forced nonrigid pendulum (the mean-field dynamics (3) are plotted for various starting points at integer multiples of the oscillation period 2π/ω). Closed loops are characteristic for stable orbits whereas irregular dots represent chaotic regions. For a BEC in a double-well, the parameters correspond to: (a) a tilt of 2µ0/Ω = 3.0, a driving frequency of ω = 3Ω, an interaction of N κ/Ω = 0.8 and a driving amplitude of 2µ1/Ω = 0.9 (i.e. a one-"photon"-resonance[26]), (b) the 3/2-"photon"-resonance with N κ/Ω = 0.1, 2µ0/Ω = 3.0, ω/Ω = 2.08 and 2µ1/Ω = 1.8, (c) all parameters as in a) except for N κ/Ω = 0.3.

FIG. 2 :

 2 FIG. 2: (color online) Quantum dynamics (N = 100) versus mean-field dynamics using the parameters of Fig. 1.a. (a) The difference of the time-averaged population imbalances Jz T /N and z/2 T as a function of 101 2 initial conditions (z0, φ0) in a two-dimensional projection of the resulting threedimensional plot (T = 100/Ω). (b) The time-averaged root-mean-square (r.m.s.)-fluctuations ∆Jz T /N of the population imbalance as a function of the initial atomic coherent state (2).

FIG. 4 :

 4 FIG. 4: (color online) Time-averaged population imbalances of quantum dynamics (N = 100) versus mean-field dynamics at the 3/2-"photon"-resonance of Fig. 1.b. (a) The difference is plotted as a function of the initial condition (z0, φ0) in a two-dimensional projection (T = 100/Ω). (b) As in a) but the mean-field dynamics are replaced by an average over the distribution of initial conditions (7).

  Fig 5.c) such a state is a bimodal distribution (for N → ∞: two delta-peaks). To numerically identify if a given wave-function |ψ is in a mesoscopic superposition, we start by searching the atomic coherent state |θ 1 , φ 1 for which | ψ|θ, φ | 2 reaches its maximum, m 1 . Around (θ 1 ,φ 1 ) we define the set R1 by | θ, φ|θ 1 , φ 1 | 2 > 10 -3 (cf. Fig 5.c). As both parts of the mesoscopic superposition (8) should hardly overlap, the second maximum m 2 = | ψ|θ 2 , φ 2 | 2 is searched outside the set R1. The set R2 is defined analogously to R1 by | θ, φ|θ 2 , φ 2 | 2 > 10 -3 . The fidelity | ψ|ψ sp | 2 still is a function of γ, taking its maximum and excluding large overlaps (R1 ∩ R2 = ∅) yields:

FIG. 5 :

 5 FIG. 5: (color online) Entanglement (10) for parameters as in Fig. 1.a (left column) and as in Fig. 1.c (right column). (a), (b): Mesoscopic quantum superpositions were identified at times τ = 5, 10, 15, . . .; the maximum value of σent is displayed for 101 2 initial conditions (z0, φ0) and for (a) short times (τ = 10) and (b) longer times (τ = 100). (c) Projection of two characteristic entangled states (with maxima m1, m2) on the atomic coherent states (2). Black (grey/red) regions: | θ, φ|ψ | 2 > 0.16 (> 0.05). Left: z0 = -0.6, φ0 = -2.764601535, τ = 80, σent ≃ 72.3% . Right: z0 = -0.98, φ0 = -2.701769682, τ = 75, σent ≃ 33.5%. In the left plot, the large blue/grey circle is a typical set R around | θ, φ with | θ, φ| θ, φ | 2 > 0.001 (cf. Eq. (9)).

We thank H. P. Breuer, Y. Castin, A. Eckardt and M. Holthaus for insightful discussions. Funding by the EU is gratefully acknowledged (CW: contract MEIF-CT-2006-038407; NT: contract MEST-CT-2005-019755).