Adrien Bernhardt

André Maximo

Luiz Velho

Houssam Hnaidi

Marie-Paule Cani

Real-time Terrain Modeling using CPU-GPU Coupled Computation

Keywords: coupled computation, sketch-based modeling, terrain modeling, GPU tessellation

. Sketch-based terrain modeling example in our tool: as each stroke is drawn (left) or manipulated (right), the terrain is tessellated in the GPU to follow the stroke. In the CPU, the quadtree data structure (middle) controls the quadrilateral patches sent to the GPU.

Abstract-Motivated by the importance of having real-time feedback in sketch-based modeling tools, we present a framework for terrain edition capable of generating and displaying complex and high-resolution terrains. Our system is efficient and fast enough to allow the user to see the terrain morphing at the same time the drawing editing occurs. We have two types of editing interactions: the user can draw strokes creating elevations and crevices; and previous strokes can be interactively moved to different regions of the terrain. One interesting feature of our tool is that terrain primitives can be interactively manipulated similarly to primitives in vector-graphics tools. We achieve realtime performance in both modeling and rendering using a hybrid CPU-GPU coupled solution. We maintain a coarse version of the terrain geometry in the CPU by using a quadtree, while a fine version is produced in the GPU using tessellation shaders.

I. INTRODUCTION

Recent years have been marked with a striking evolution of sketch-based modeling tools [START_REF] Igarashi | Teddy: a Sketching Interface for 3D Freeform Design[END_REF], [START_REF] Kalnins | WYSIWYG NPR: Drawing strokes directly on 3D models[END_REF], [START_REF] Schmidt | ShapeShop: Sketch-based Solid Modeling with BlobTrees[END_REF], [START_REF] Nealen | FiberMesh: Designing Freeform Surfaces with 3D Curves[END_REF], which have the interesting property of benefiting greatly from real-time feedback while editing. As the object to be modeled increases in complexity, the modeling tool becomes more involving to develop. In the special case of terrain modeling, exemplified in figure 1, the challenge is to create in real-time realistic terrain features with intuitive strokes and guided landscape design.

Synthetic terrains are applied in many applications, such as: crowd and flight simulators, real movies using chroma-key effects, film animations, and computer games. Regardless of the application, synthetic terrains often have direct influence in the visual and artistic aspects to be accomplished. The common strategy to produce terrain models is by procedural generation [START_REF] De Carpentier | Interactive GPU-based Procedural Heightfield Brushes[END_REF], [START_REF] Yeguas | Example-based Procedural Modelling by Geometric Constraint Solving[END_REF] through texture synthesis and height-field captured examples. The idea is to achieve a high level of realism through Digital Elevation Models (DEMs) of real terrains [START_REF] Zhou | Terrain Synthesis from Digital Elevation Models[END_REF]. Although this strategy produces convincing results, it lacks controllability and restricts the result to variations of existing terrain primitives.

One way to improve controllability over the terrain generation is to use image samples from a procedural input [START_REF] Perlin | Live Paint: Painting with Procedural Multiscale Textures[END_REF]. The idea behind this technique is to allow the user to intuitively paint bumps at different resolution levels, enhancing the terrain modeling tool while maintaining realistic results. Another step further usability and control is to permit angles different than 90 • when painting details on the terrain [START_REF] Keller | Introducing ZBrush[END_REF]. As the tool increases in controllability, the terrain modeling tool becomes more user-controlled and less automatic. The zenith of this process is an interface closer to pen-and-paper sketching [START_REF] Cohen | Harold: a World Made of Drawings[END_REF], allowing total control over terrain primitives exchanging physical realism by artistic freedom.

Parallel to the question of balancing controllability and realistic results, terrain modeling tools can be analyzed by two different computational perspectives: CPU-based modeling strategies [START_REF] Gain | Terrain Sketching[END_REF], [START_REF] Smelik | Integrating Procedural Generation and Manual Editing of Virtual Worlds[END_REF] and GPU-based shading and solver techniques [START_REF] Anh | Physically based Hydraulic Erosion Simulation on Graphics Processing Unit[END_REF], [START_REF] Peytavie | Arches: a Framework for Modeling Complex Terrains[END_REF], [START_REF] Hnaidi | Feature based Terrain Generation using Diffusion Equation[END_REF]. On the one hand, the CPU plays the role of a controller and its strategies focus mainly on integrating manual editing with complex deformations and simulations. On the other hand, the GPU offers massive parallelism and its techniques focus on shading landscapes and performing computational intensive tasks, such as solving diffusion equations. The gap between one unit and the other becomes clear when analyzing the data structure involved on each computation. While the CPU allows a broad range of inter-connected data without penalizing access, the GPU requires coherent memory requests imposing data regularity.

Contribution: In this paper, we present a new modeling framework for terrain editing by coupling CPU-GPU computational perspectives. Differently from other approaches, discussed in section II, our modeling tool uses the CPU not only to control the initial level-of-detail (LOD) visualization but also to start and stop the terrain generation done by the GPU. The combination of the different computational granularities makes the data paths between the two units more complex, as explained in section III.

Our approach partitions the terrain with a view-and heightmap-dependent quadtree [START_REF] Finkel | Quad Trees a Data Structure for Retrieval on Composite Keys[END_REF], as shown in figure 1 (middle). This lightweight data structure controls the quadrilateral patches to be rendered, serving two main purposes: it allows a first coarse LOD analysis of the terrain; and it frees the CPU from the costly task of generating the entire heightmap data.

Next, we employ the GPU to solely generate and maintain the terrain data on its own texture memory. The generation is accomplished by a fast multigrid solver [START_REF] Hnaidi | Feature based Terrain Generation using Diffusion Equation[END_REF] specifically designed for terrains. The generated data is constantly used in our approach to morph the original rendering of low-resolution quadtree patches in smooth high-resolution triangles through the graphics-card tessellation shaders. The tessellator allows for a second fine LOD analysis of the terrain, balancing the decision between the units. Although the whole data is updated and used by the GPU, the management is done by the CPU.

The final result of our coupled solution is a real-time terrain modeling tool based on sketches and capable of dynamically generate and visualize multiresolution heightmaps (see section IV). The framework presented is discussed in more details in section V, and section VI concludes our work giving future research directions.

II. RELATED WORK

In many editing tools, especially sketch-based modeling, it is important to have real-time feedback to help improve the editing quality [START_REF] Olsen | Sketch-based Modeling: A Survey[END_REF], [START_REF] Igarashi | Teddy: a Sketching Interface for 3D Freeform Design[END_REF], [START_REF] Igarashi | A Suggestive Interface for 3D Drawing[END_REF]. This importance is emphasized particularly in sketch-based terrain modeling, being able to see the terrain morphing at the same time the drawing editing occurs constitutes a great user experience, an example is shown in figure 1 (left).

The quality of the modeling tool depends mainly on three aspects: controllability, intuitiveness and responsiveness. The first makes the actual modeling possible, the second makes the tool more natural and easy to learn, and the third enhances the creativity process when the tool provides real-time feedback. One interesting example combining these three aspects is the FIBERMESH [START_REF] Nealen | FiberMesh: Designing Freeform Surfaces with 3D Curves[END_REF]. In this tool, the user is able to create freeform surfaces with 3D curves in real-time with natural and unconstrained design. However, there are two drawbacks in FIBERMESH: the high pre-computational cost required when adding or removing 3D curves; and the relative small models that the tool can handle. Like FIBERMESH, our system builds the surface by interpolating curves but allowing real-time interaction with more complex models.

Terrain modeling tools also benefit from intuitiveness and natural interactions. The usage of sketch-based terrain modeling exploring this aspect was presented by Watanabe and Igarashi [START_REF] Watanabe | A Sketching Interface for Terrain Modeling[END_REF], and further inspired the work of Gain et al. [START_REF] Gain | Terrain Sketching[END_REF]. In Gain et al.'s work, the user can draw strokes defining mountain silhouettes and modify previous strokes forming the base of mountains. Their CPU-based modeling framework generates terrain primitives extracting, at the same time, noise from the user strokes to create more realistic terrains. Although we do not extract noise from the user strokes in our tool, we have two types of editing interactions, in the same spirit of Gain et al.'s work: the user can draw strokes creating mountain elevations or crevices; and previous strokes can be interactively moved to different regions of the terrain, as shown in figure 1 (right). Moreover, we use a CPU-GPU coupled method to drastically improve the performance of our tool, generating terrains two orders of magnitude faster than Gain et al.'s work.

One work dealing with static real terrains and focused on rendering is the work of Losasso and Hoppe [START_REF] Losasso | Geometry Clipmaps: Terrain Rendering Using Nested Regular Grids[END_REF]. In their approach, called geometry clipmaps, they define a regular-grid hierarchy centered about the viewer instead of a world-space quadtree. This choice simplifies both their pyramid compression scheme and inter-level continuity while the viewer is moving. The geometry clipmaps uses the CPU to (1) construct LOD regular regions (see [21, figure 2]) based on viewer distance, (2) tessellate the regions in fine-to-coarse order and using triangle strips (see [21, figure 3]) to exploit hardware occlusion culling and vertex caching, (3) update the terrain maintained in CPU memory via decompression. There is one simple data path as the GPU is only used to interpolate height values for visual continuity in the vertex shader. In a revised version of this idea [START_REF] Asirvatham | Terrain Rendering using GPUbased Geometry Clipmaps[END_REF], the height-value interpolation is improved in the GPU by elevation-map upsampling, terrain detail synthesis and normal-map computation but remaining with a single data path between the units. In our approach, we are interested in dynamically changing terrains and use different techniques: a CPU-based quadtree data structure rather than a nested regular grid; and a GPU-based runtime terrain management for both tessellation and terrain modification.

Realistic terrains can also be created through modeling by employing simulation of natural phenomena. Anh et al. [START_REF] Anh | Physically based Hydraulic Erosion Simulation on Graphics Processing Unit[END_REF] presented a GPU-based method for hydraulic erosion simulation capable of simulating water transportation and erosion one order of magnitude faster than CPU-based methods. The work of Št'ava et al. [START_REF] Št'ava | Interactive Terrain Modeling using Hydraulic Erosion[END_REF] presents a comprehensive set of GPU-accelerated erosion tools for real-time interactive terrain modeling, indicating the usage of graphics hardware for the solution of physically-based equations. Although in our modeling system we do not have simulation effects, we use the GPU computational power to generate the terrain by solving a different type of physically-based equation.

The terrain modeling in our approach is accomplished by combining the multigrid GPU solver of Hnaidi et al. [START_REF] Hnaidi | Feature based Terrain Generation using Diffusion Equation[END_REF] with an adaptive tessellation-based rendering shader capable of handling dynamic heightmaps. The main contribution of Hnaidi et al.'s work is to introduce a GPU-based multigrid diffusion equation solver for terrains relying on Laplacian equations, which interpolates not only heights but also amplitude and frequency noise parameters. Our modeling tool uses Hnaidi et al.'s solver to allow an interactive edition and manipulation of complex terrain primitives.

III. REAL-TIME TERRAIN MODELING

The main contribution of this paper is a real-time terrain modeling system, combining a novel CPU-GPU coupled solution with a natural sketch-based interaction. In this section, we describe both the coupled computation and this interaction.

A. CPU-GPU Coupled Computation

Creation of terrain models in real-time involves dealing with dynamically changing data that increases exponentially depending on the terrain resolution. However, high resolution is only interesting in parts of the terrain with a high level of detail. In our terrain modeling framework, presented in figure 2, we make use of this observation dividing the terrain generation process over a multiresolution pyramid. The base level of this pyramid, i.e. coarse resolution, is used whenever the details are low or do not appear in the current point-ofview. The top level of the pyramid, i.e. fine resolution, is used whenever the details are high and dominant in the current shading of the terrain model. Fig. 2. The terrain modeling framework of our tool. While the CPU controls the multigrid solver depending on drawings and sends terrain patches to be rendered, the GPU generates the terrain in texture memory and tessellates its own produced heightmap. This setting requires three different data paths: adaptive rendering using a quadtree (red); solver iterations control (blue); and updating coarse-resolution details on the quadtree (green).

There are two main user interactions related to the multiresolution pyramid: draw of strokes controlling the terrain generation process; and camera movement changing the rendering process. These two interactions guide two complementary approaches of our terrain modeling tool. First, a coarse version of the terrain geometry is maintained in the CPU using a quadtree, as shown in figure 1 (middle), where regions closer to the viewer are subdivided more than far regions. This simple and lightweight data structure fits the CPU main role of data control, while allowing it to send adaptive quadrilateral (or quad) patches to the GPU (illustrated by the camera-thenrender path in red in figure 2). Second, a fine version of the terrain is produced in the GPU using two programmable shaders within the tessellator stage. One of them is responsible to subdivide regularly each quad patch, i.e. the quad leaf node sent by the CPU, while the other reads the height values from a texture that changes interactively based on manipulation and drawing of strokes (illustrated by the draw-then-control path in blue in figure 2). It is important to note that the different data paths in our framework have an influence in each other, e.g. drawing changes the heightmap that updates the quad patches that changes rendering.

The first phase of our algorithm is to update the quadtree data structure using a level-of-detail mechanism. Our LODbased approach considers the viewing projection of the bounding box of each quadtree node (see figure 3). This bounding box is constructed by reading the minimum and maximum height values that fall inside the node (illustrated by the update-bounding-box path in green in figure 2). We use a linear metric on the number of pixels of the bounding-box projection in order to make our LOD approach both view and height-values dependent. Using this metric, the quadtree node is recursively subdivided until a given criterion is met, effectively increasing the resolution of the current terrain and moving upward in the multiresolution pyramid. In our framework, we have used the rectangular shape of the bounding-box projection to determine the quadtree node subdivision. The recursive creation of child nodes, four at a time, is done in order to guarantee that each projection shape has about the same size (in pixels). Moreover, the screen rectangles are used to determine if a quadtree node is completely occluded and can be discarded, i.e. not sent for rendering as in the viewfrustum culling technique. The result is an adaptive quadtree controlling the minimum resolution level imposed to the GPU (depending on each quad node subdivision), which changes interactively depending on camera properties. In the GPU, we use regular subdivision of patches (each patch corresponds to a quadtree node sent for rendering) but a more refined LOD technique could also be applied, allowing a second form of LOD in the top level of the multiresolution pyramid.

The second phase of our algorithm comprehends the following modules (depicted in figure 2): the multigrid solver; and the tessellation shaders module. The former provides the height values of the terrain to be used by the latter module.

The multigrid solver module, based on the work of Hnaidi et al. [START_REF] Hnaidi | Feature based Terrain Generation using Diffusion Equation[END_REF], computes the heightmap of the terrain through a sequence of increasing resolution textures, up to an arbitrary size, using a two-step algorithm. The first step solves a Poisson-based diffusion system to interpolate gradient constraints. The resulting gradients and height constraints are then used by the second step to compute the final heightmap employing a similar diffusion system. Both systems are solved using a multigrid strategy, detailed in figure 5, which first solves a coarse-resolution system to then increase the resolution and start solving again. While the dyadic upsampling is done once per level, the number of solver iterations is normally the same for each level and depends on the desired accuracy. The input of the multigrid solver module is given by stroke drawing and is explained in the next section. The output of this module is a multiresolution terrain stored as a mipmap texture pyramid (see figure 4). While the GPU is responsible for this module, the CPU is able to start and stop solver iterations within drawing and rendering events, keeping the GPU busy. The multiresolution texture scheme of our tool. The terrain is stored as a mipmap texture pyramid that changes over time, depending on the multigrid solver iterations done so far.

The tessellation shaders module uses the mipmap-pyramid texture generated by the solver module in a producer-consumer strategy, synchronized and controlled by the CPU. After each new resolution level is constructed by the solver, a shader uniform value is updated to the current maximum level. The CPU also controls the minimum level through the quadtree refined nodes sent as quad patches to the GPU.

The tessellation shaders module employs two stages of the graphics-hardware pipeline: the tessellation control shader, used to determine the subdivision level done by the tessellator fixed functionality; and the tessellation evaluation shader, used to specify properties of the vertices created by the tessellator. In our terrain modeling tool, the control shader uses a constant uniform value determined by the CPU to subdivide regularly each quad patch. Depending on the subdivision, the resolution of the patch can increase or remain the same. The evaluation shader reads the appropriate level of the mipmap-pyramid texture to place each generated vertex at the correct height position. The multiresolution texture is also used by the CPU, but only a small resolution of it (we use 64 × 64 from the 6th solver iteration) since the quadtree minimum leaf size is much bigger than the texel from the highest resolution texture.

In the tessellation shaders, we use a straightforward strategy for subdivision that does not consider different tessellation levels across patch boundaries. This may lead to T-junctions, i.e. edges connected forming a T. In our experiments, we generate a large number of triangles through the tessellation, making this possible T-junctions imperceptible. However a more elaborated strategy for low-level tessellation can use irregular subdivision ac-cross patches to avoid T-junctions. Our CPU-based LOD guarantees that visible quad patches have at most 2 neighbor patches (invisible patches can be arbitrarily big) and this fact can be used to construct a bit-mask for each patch, to be sent to the GPU, imposing a twice coarser subdivision for patches with bigger neighbors. This approach can be thought as a second form of LOD done by the GPU that is local-geometry dependent.

The integration of the first CPU phase and the second GPU phase is done after the quadtree is refined and before rendering. The leaf nodes of the current quadtree are transformed into quad patches, where each patch contains its corresponding mipmap level and texture position. From this patch, the GPU is able to further increase (but not decrease) the terrain resolution in the tessellation control shader, up to a specified resolution depending on the last solver iteration. This coupled solution defines the entire multiresolution pyramid, where the base level is determined by the CPU and the top level is computed and generated by the GPU.

B. Sketch-based and Interactive Edition

Interaction with terrain models in real-time requires primitives that are easy to create, manipulate and understand. Creating mountains or plains by stipulating each small piece of the terrain can be a cumbersome and slow work. Our terrain modeling system allows intuitive and natural interaction through simple sketches drawn over the current camera viewing plane. These sketches act by pulling and pushing the terrain surface. The sketches can be created, moved or deleted; effectively changing the terrain. Depending on the desired primitive to create or modify, the user can move the camera to choose the best view of the terrain. The camera zoom can also be used to interact with terrain primitives at different scales.

There are two main primitives supported by our sketchbased modeling tool: mountains and crevices primitives (drawing up or down the terrain level), and plateaux primitives (drawing in the same terrain level). These primitives can be seen in a terrain example shown in figure 6.

To interactively morph the terrain as each primitive is drawn or manipulated, the modeling strokes are converted into constraints used by the GPU multigrid solver in the second phase of our algorithm (explained in the previous section). The 3D strokes are projected and rasterized on top of the terrain creating an input 2D image used by the multigrid solver to specify its constraints. This process is further detailed in figure 5. The 3D curves before rasterization represent terrain primitives and are similar to curves in a vector-graphics tool.

The projection of 3D curves depends on the primitive been drawn. In the case of a plateaux primitive, the stroke is projected on the horizontal plane defined by the first point of the drawing that touches the terrain. In the case a mountain Fig. 5. The landscape primitives specified by 3D strokes are converted to constraints through rasterization and stored on a 2D constraint texture. This texture is downsampled to be used in a solve-then-upsample process, i.e. the multigrid solver, where it first solves for the small resolutions before extrapolating the solution for finer resolution systems. The final solution is the heightmap of the terrain in a multiresolution pyramid stored as a mipmap texture.

primitive, the first and last point of the drawing touching the terrain determine a vertical plane where the stroke is projected. However, in our real-time modeling scenario, the last point is unknown while the drawing editing occurs. One simple solution is to project the extremity of the stroke been drawn onto the terrain to obtain a dynamic last point. Although this solution successfully changes the terrain, the result is confusing to the user since both terrain and drawing are changing dynamically (one stroke point updates the terrain that updates where the next stroke point is projected). Instead of projecting the stroke extremity to the dynamic terrain, we use the depthmap of the scene when the user starts drawing as an approximation of the current terrain. By projecting onto this depthmap, with the appropriate transformations, the moving last point is now based on a fixed terrain, matching what the user expects to modify.

The manipulation of 3D curves can be done independently from each other and from new drawings. The first or last points used to compute the projection plane for primitives can be manipulated with drag-and-drop interaction changing the projection and thus the primitive. These points act as handles that can be moved morphing the terrain to a different configuration.

The creation and manipulation of different primitives are handled by the CPU. The solver in the GPU receives the final projection by reading from the 2D constraint texture (see figure 5). The different type of strokes and the number of terrain primitives affect little the performance of the solver.

IV. RESULTS

We have tested our modeling system with different heightmap resolutions and number of solver iterations. Table I presents the terrain creation solver and tessellation-based rendering timings and GPU memory consumption for the tested terrain size and respective number of iterations. The timings are given using a 1024 × 768 pixel viewport and considering the camera constantly moving. All timings were performed in an Intel Core2Quad CPU and an nVidia 480 GTS GPU using OpenGL 4.1. For these results we have drawn several strokes, changing the terrain shape and inducing the GPU multigrid solver, and considered the timing to create and tessellate the maximum resolution. The timings show that our system can compute the solver and tessellate high-resolution terrains in real-time (at about 20 frames per second) allowing for still higher resolutions when working with adaptive tessellation.

Analyzing the performance results further, the ability to control incrementally the creation of the terrain constitutes another important characteristic of our approach. The solver iterations can be used to specify the resolution levels and terrain features precision, which affects both timing and memory consumption. More iterations generate improved and higher resolution terrains, at the cost of more GPU time and memory. The memory usage in our system refers to the mipmappyramid texture stored in the GPU.

One terrain example modeled with our tool can be seen in figure 6. In this example appears different terrain primitives drawn: mountains, plateaux and crevices. In addition to terrain primitives, our tool allows several user interactions, including sketching and manipulation of previous strokes and camera movement.

Adding the GPU to the process of creating and tessellating the terrain (not only for rendering) makes our approach to run at 27.8 ms for a 512 × 512 grid. While Gain et al. [START_REF] Gain | Terrain Sketching[END_REF] report 2.3 seconds on a similar CPU to generate the same terrain grid. That is, our approach generates terrains two order of magnitude faster than Gain et al.'s work employing one single unit. Fig. 6. Terrain example generated by our modeling system (north and south views). The blue strokes were drawn to create mountains and crevices, and the red strokes were drawn to create plateaux in the midst.

Compared against Hnaidi et al.'s approach [START_REF] Hnaidi | Feature based Terrain Generation using Diffusion Equation[END_REF], our coupled solution introduces the usage of tessellation shaders and a CPU-based method to control and synchronize between drawing and rendering events. One direct advantage of this approach comes from the high-cost memory bandwidth between the CPU and the GPU; by using the CPU only for control and synchronization, the communication between the units becomes small (few kilobytes in our system) even though the GPU is handling large amounts of data (hundreds of megabytes for high-resolution terrains). Without the CPU-GPU coupled solution, the modeling framework changes from computationbound (regarding FLOP count) to I/O-bound (regarding data transfers), compromising the real-time performance of the terrain modeling tool.

Another advantage is that while the multigrid GPU solver acts in a global scope, generating and improving the entire terrain, the CPU quadtree adaptation and the GPU tessellation (see an illustrative example in figure 7) act in a local scope, refining the terrain on demand for rendering. This strategy alleviates the burden on the CPU, in terms of editing and visualization, and on the GPU, in terms of control and synchronization, providing a comprehensive sketch-based terrain modeling tool. Without one unit or the other, the tool is less capable of dealing with high-resolution and dynamically changing terrains.

V. DISCUSSIONS

The majority of procedural editing and feature extraction tools for terrain modeling are based on a top-view interaction [START_REF] Zhou | Terrain Synthesis from Digital Elevation Models[END_REF], [START_REF] Cohen | Harold: a World Made of Drawings[END_REF], [START_REF] Olsen | Sketch-based Modeling: A Survey[END_REF]. Although interesting, this approach limits the artistic freedom of the tool. Our modeling system is based on a 3D camera with a first-person perspective and sketchbased drawing interaction over the moving view plane. This method allows the user to draw landscape silhouettes at any distance and with more freedom.

The drawing interaction is only convincing when the tool responds in real-time. Inspired by FIBERMESH [START_REF] Nealen | FiberMesh: Designing Freeform Surfaces with 3D Curves[END_REF], our tool has controllability and real-time responsiveness. However, to deal with high-resolution terrains, we use a different type of solver to generate the model. While FIBERMESH uses a direct solver with a high pre-computational cost, we use a fast multigrid solver in the GPU that is suitable for our goal.

Our GPU-based solver hints for a multiresolution hierarchy of terrains in the form of a texture pyramid. This pyramid is constructed top-down and can be used profitably for level-ofdetail control on both the CPU and the GPU (as explained in section III). Aiming at only visualizing static and real terrains, the technique from the geometry clipmaps [START_REF] Losasso | Geometry Clipmaps: Terrain Rendering Using Nested Regular Grids[END_REF], [START_REF] Asirvatham | Terrain Rendering using GPUbased Geometry Clipmaps[END_REF] uses a better solution for LOD control, building a compressed multiresolution pyramid that is accessed bottom-up. However, we are interesting in a flexible modeling interface for dynamic terrain generation and the geometry clipmaps' technique is not appropriate.

Between creating and tessellating the terrain in the GPU, the creation is more costly (see table I) since we use a simple regular subdivision approach in the tessellation control shader. In the case of a complex second form of LOD in the GPU, to allow even higher terrain resolutions, the bottleneck may become the terrain tessellation.

The two modules of the CPU, controlling the solver iterations and building the quadtree data structure, depend on the user input, drawing of strokes and moving the camera, and are intrinsic serial computations. Although the two modules can not be implemented by the GPU, they could be suppressed, making the modeling tool more simple to implement but less interactive.

VI. CONCLUSIONS

We have presented a real-time terrain modeling tool combining two strategies in the GPU with a lightweight CPU-based data structure. Our tool is capable of dynamically generate heightmaps with adaptive resolution, that is, our tessellation shaders are able to generate on-the-fly different parts of the terrain at different resolutions.

One interesting aspect of our CPU-GPU coupled computation solution is the usage of the CPU to control the GPU solver iterations and stop at a certain resolution and then resume computing when the GPU is idle. Another interesting feature of our method is the balance between terrain generation in the CPU and in the GPU -we can control this balance by simply changing the quadtree refinement. With these features, the user can draw strokes and see at the same time the terrain morphing to the drawing. Terrain primitives, such as mountains and clefts, are controlled seamlessly in our framework.

Although with interesting features, our modeling framework has a wide variety of future work directions. First, there is room for optimizations in the GPU solver: instead of performing the initial iterations in the GPU, the CPU is more adequate to compute these low-resolution mipmaps and use them in its quadtree data structure (avoiding transfering data back), while the GPU is not underutilized with small images. Second, the current usage of the mipmap-pyramid texture in the tessellation evaluation shader indicates the possibility of using anisotropic filtering with a gradient vector depending on the solver. Third, the tessellation shaders allow a second fine LOD analysis of the terrain that is interesting to explore. Finally, adding normal and fractal texture painting or simulation effects in our sketch-based modeling tool can increase the realism of the generated terrains. In short, we believe this work has the potential to become an effective terrain modeling tool in practical applications.

Fig. 3 .

 3 Fig. 3.Illustration of a quadtree node and its corresponding bounding box. The minimum and maximum height values are used to specify the limits inside the current node. The quadtree is refined depending on the projection shape of each node.

  Fig. 4.The multiresolution texture scheme of our tool. The terrain is stored as a mipmap texture pyramid that changes over time, depending on the multigrid solver iterations done so far.

Fig. 7 .

 7 Fig. 7. Illustration of coarse (top) and fine (down) version of a terrain tessellated by the GPU. The coarser resolution corresponds to the quadtree (shown on the right), that is the minimum resolution level imposed by the CPU, and the finer resolution is the maximum resolution level.

ACKNOWLEDGMENTS

We would like to acknowledge the grant of the first author provided by the PPF Interactions Multimodales, and the grant of the second author provided by Brazilian agency CNPq (National Counsel of Technological and Scientific Development).

We also would like to thank the anonymous reviewers for their valuable comments and suggestions.


Alexander Mattioli Pasqual
email: ampasqual@gmail.com

Vincent Martin

INVERSE SOUND SOURCE RECONSTRUCTION BY EXTERIOR SPHERICAL ACOUSTICAL HOLOGRAPHY WITH MODEL ADAPTATION

Spherical acoustical holography (SAH) can be used to reconstruct the 3D sound field from acoustic data on a spherical surface, called the hologram sphere. In exterior SAH, the hologram sphere encloses all the sound sources and free-field conditions are assumed. Then, the acoustic variables can be reconstructed at any point exterior to the smallest sphere surrounding the sources (called the source sphere) by using a radiation model written in terms of spherical Hankel functions and spherical harmonics. This work deals with the inverse problem of reconstructing the radial particle velocity on the source sphere from the sound pressure field on the hologram sphere. This requires prior knowledge of the position of the source sphere inside the hologram sphere; as far as possible, these spheres must be made concentric to simplify the reconstruction equations. If the relative position of the hologram and source spheres is not known, an estimated position must be used. This may yield an inadequate propagator, so that the acoustic velocity field is poorly reconstructed. This contribution discusses the role of the relative position of the spheres in inverse source reconstruction. In addition, some strategies to improve the radiation model by adapting the source position from the hologram data are investigated.

Introduction

Spherical acoustical holography (SAH) is an imaging technique to reconstruct a 3-D sound field from available data on a spherical surface (the "hologram sphere"), which requires a wave propagation model [START_REF] Williams | Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography[END_REF]. The reconstruction equations differ essentially whether the problem of interior or exterior SAH is dealt with. In interior SAH, the sound sources surround the acoustic domain, and thus they are located outside the hologram sphere. In exterior SAH, the sources are located inside the hologram sphere and free-field conditions are assumed, so that the acoustic field can be reconstructed at any point exterior to the smallest imaginary sphere surrounding the sources (the "source sphere"). This work deals with exterior SAH. In particular, we focus on the inverse problem of reconstructing the radial particle velocity on the source sphere from the sound pressure on the hologram sphere.

Since the sound pressure is known on a spherical surface, it is naturally subjected to spherical harmonic analysis [START_REF] Williams | Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography[END_REF][START_REF] Arfken | Mathematical Methods for Physicists[END_REF]. The spherical harmonics representation of a sound field depends on the position and orientation of the coordinate system, i.e., translations and rotations have significant effects on the expansion coefficients [START_REF] Gumerov | Recursions for the computation of multipole translation and rotation coefficients for the 3-D Helmholtz equation[END_REF]. Moreover, unlike rotation of the coordinate axes, translating the origin also changes the number of harmonics that must be retained to obtain an accurate representa-tion of the acoustic field. Some recent studies [START_REF] Rafaely | Spatial alignment of acoustic sources based on spherical harmonics radiation analysis[END_REF][START_REF] Deboy | Acoustic centering and rotational tracking in surrounding spherical microphone arrays[END_REF][START_REF] Deboy | Acoustic center and orientation analysis of sound-radiation recorded with a surrounding spherical microphone array[END_REF] have described the effects of the translation on the spherical harmonics decomposition of sound pressure fields, and proposed some cost functions to obtain an optimal decomposition center that minimizes the high-order spherical harmonics content, so that a more compact description of the acoustic field is obtained. However, these works do not deal explicitly with inverse sound source reconstruction.

In this paper, the origin of the coordinate system is the center of the source sphere. Clearly, its position inside the hologram sphere is a parameter of the radiation model, which is assumed to be known in conventional exterior SAH; as far as possible, the spheres must be made concentric in order to simplify the reconstruction equations. However, in practice, the source may not be precisely positioned inside the hologram sphere. If so, a nominal position can be used, which may degrade the holography results. This work investigates the effects of the source position on the inverse reconstruction of its radial particle velocity. In addition, we discuss some strategies to estimate the source position from the sound pressure field it radiates on the hologram sphere. The paper is organized as follows: Section 2 provides a description of the problem to be addressed. Section 3 deals with the spherical harmonic analysis of the pressure field on the hologram sphere with a given but arbitrary position of the decomposition center. Section 4 focuses on the reconstruction equations and on two cost functions that can be used to estimate the source position. Finally, Sec. 5 gives simulation results in order to illustrate the ideas presented throughout the paper.

From a general point of view, the problem addressed in this paper consists in adapting a parameter of a propagation model in order to minimize a given cost function. Here, the model describes the sound radiation from the source, and the parameter that is to be adapted is its position. If an adequate cost function is used, the optimal position leads to the lowest reconstruction error, which corresponds obviously to the real position of the source. It is worth mentioning that Martin et al. [START_REF] Martin | Numerical simulation of acoustic holography with propagator adaptation: Application to a 3D disc[END_REF] applied this approach to a related acoustical problem, that of planar near-field acoustical holography with generalized impedance boundary conditions.

Problem description

Throughout this paper, a harmonic time dependence e -ιωt is assumed, where ι ≡ √ -1, ω is frequency, and t is time. Lower and upper case bold letters indicate vectors and matrices, respectively.

Let Γ h and Γ s be imaginary spheres with radius r h and r s , centered at C and O, respectively. They enclose all the sources in the acoustic domain, such as the hermetic compressor depicted in Fig. 1a. Γ s is the smallest sphere surrounding the sources, which will be referred to as "source sphere". Also, let P be the space of continuous functions defined on Γ h and p(x) : |x| ≥ r s be the free-field sound pressure produced by the sources, where x is the position vector in relation to the origin O. The sound pressure is given only at x h ∈ Γ h , which defines the "hologram sphere". Exterior SAH can be used to reconstruct the sound field at any point exterior to or on Γ s from p h ≡ p(x h ) ∈ P. In particular, we are interested in the inverse problem of determining the radial particle velocity v s ≡ v(x s ) : x s ∈ Γ s , as shown in Fig. 1b, where d ≡ -→ CO. This requires a radiation model that satisfies the Helmholtz equation and the geometrical configuration of the problem. Since d is clearly a parameter of the model, it must be known a priori. As far as possible, the sources position must be set so that Γ h and Γ s are concentric (d = 0) to simplify the reconstruction equations. However, in practice, the sources may not be precisely positioned inside the hologram sphere. If so, a nominal position d ≡ -→ C Õ can be used in the holographic reconstruction, which leads to the radial particle velocity (ṽ s ) on a translated source sphere ( Γs ), as illustrated in Fig. 1c. Γs will be referred to as "equivalent source sphere". This work investigates the effects of δd ≡ d -d on the inverse sound source reconstruction, as well as some strategies to estimate d from a given p h .

Spherical harmonic analysis of the sound pressure field

Let r ≥ 0, 0 ≤ θ ≤ π, and 0 ≤ ϕ < 2π be spherical coordinates. Accordingly, the position vector is x = r sin θ cos ϕe x + r sin θ sin ϕe y + r cos θe z , where e x , e y , and e z are the versors of the Cartesian coordinate system centered at the origin Õ. If the equivalent source sphere in Fig. 1c has a number sufficiently large of degrees of freedom, it is able to produce a given sound pressure field p h on Γ h regardless of d. The sound pressure it radiates, p(x) : |x| ≥ r s , can be written as [START_REF] Williams | Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography[END_REF]

p(x) = N ∑ n=0 n ∑ m=-n C m n (ω, d)h n (kr)Y m n (θ, ϕ) = N ∑ n=0 n ∑ m=-n γ m n (r, d)Y m n (θ, ϕ), ( 1
)
where k is the wave number, C m n are expansion coefficients, h n are the spherical Hankel functions of the first kind, Y m n are the spherical harmonics, and

γ m n (r, d) = C m n (ω, d)h n (kr) (2)
is the so-called "spherical wave spectrum" of p on a sphere defined by |x| = r [START_REF] Williams | Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography[END_REF]. Note that C m n are obtained from p h , and thus they depend on d.

Here, we only deal with band-limited functions, i.e., C m n = 0 if n > N , so that the series has N ′ ≡ (N + 1) 2 terms. In fact, the order N required to ensure that p(x h ) = p h : xh ∈ Γ h depends on d in a difficult to predict way [START_REF] Gumerov | Recursions for the computation of multipole translation and rotation coefficients for the 3-D Helmholtz equation[END_REF][START_REF] Rafaely | Spatial alignment of acoustic sources based on spherical harmonics radiation analysis[END_REF][START_REF] Deboy | Acoustic centering and rotational tracking in surrounding spherical microphone arrays[END_REF][START_REF] Deboy | Acoustic center and orientation analysis of sound-radiation recorded with a surrounding spherical microphone array[END_REF]. In this work, N is made large enough so that p(x h ) ≈ p h for the whole d range considered. Clearly, if d = d, one has p(x) = p(x) ≈ p(x) : |x| ≥ r s .

Concentric spheres

It is known that Y m n are orthonormal functions over the unit sphere (see, e.g., Ref. [START_REF] Arfken | Mathematical Methods for Physicists[END_REF]). Therefore, if Γ h and Γs are concentric ( d = 0), the expansion coefficients in Eq. ( 1) are given by

C m n (ω, 0) = 1 h n (kr h ) ∫ π 0 ∫ 2π 0 p h Y m n (θ, ϕ) * sin θdθdϕ, ( 3
)
where the asterisk indicates complex conjugate. Since C m n depends on frequency and |h n | increases rapidly with n for small arguments, it is helpful to use the set of coefficients γ m n (r h , 0), which is the spherical wave spectrum of p h . In this case, the Parseval's theorem,

∫ π 0 ∫ 2π 0 |p h | 2 sin θdθdϕ = ∞ ∑ n=0 n ∑ m=-n |γ m n (r h , 0)| 2 ≈ N ∑ n=0 n ∑ m=-n |γ m n (r h , 0)| 2 , (4)
shows that γ m n rather than C m n gives the contribution of the mode nm to the total energy of p h .

Non-concentric spheres

For the sake of notation convenience, it is useful to combine the coefficients C m n into the column vector c ∈ C N ′ , and the functions h n (kr)Y m n (θ, ϕ) into the column vector s(x), so that c i ≡ C m n and s i (x) ≡ h n (kr)Y m n (θ, ϕ), with i = n 2 + n + m + 1 for linear indexing. Then, Eq. (1) becomes

p(x) = s(x) T c. (5)
If d ̸ = 0, the sphere Γ h is described by a radial coordinate, r, that depends on (θ, ϕ). As a consequence, unlike for d = 0, there is not a simple integral to evaluate c i because s i are not orthogonal over Γ h . In this case, c can be obtained by solving the least-squares minimization problem,

min ĉ R(ĉ) = min ĉ ||s(x h ) T ĉ -p h ||, ( 6
)
where

ĉ ∈ C N ′ , || • || = ⟨• , •⟩ 1/2 is the 2-norm, ⟨• ,
•⟩ is the inner product, and R is the residual norm. Throughout the paper, the norms and inner products of P and C N ′ will be denoted by || • || and ⟨• , •⟩.

We omit indices indicating the space because it will become clear from the context. Problem ( 6) is convex and has an analytical solution given by

c = A -1 b, (7)
where

A ∈ C N ′ ×N ′ with A ij = ⟨s j (x h ), s i (x h )⟩, and b ∈ C N ′ with b i = ⟨p h , s i (x h )⟩.
Notice that, for low frequencies and large N , the matrix A ij = ⟨s j (x h ), s i (x h )⟩ is ill-conditioned due to the very different decay rates of h n . In order to overcome this difficulty, γ m n (r h , d) can be used instead of C m n . In matrix/vector notation, Eq. ( 2) becomes

γ = Hc, (8)
where

H ∈ C N ′ ×N ′ is a diagonal matrix with H ii ≡ h n (kr h ), and
γ i ≡ γ m n (r h , d), i = n 2 + n + m + 1.
Then, by letting γ ≡ Hĉ, substitution of Eq. ( 8) into problem (6) leads to

min γ R(γ) = min γ ||s(x h ) T H -1 γ -p h ||, ( 9
)
whose solution is given in Eq. ( 7), with

γ = γopt = A -1 b, A ij = ⟨ s j (x h )H -1 jj , s i (x h )H -1 ii ⟩ , and b i = ⟨ p h , s i (x h )H -1
ii ⟩ . For r h >> | d|, the arguments of the spherical Hankel functions in H and s(x h ) are nearly the same. If so, d has no appreciable effects on the conditioning of problem (9), which is well-posed.

Reconstruction of the radial particle velocity 4.1 Reconstruction equation

The solution of problem (9) leads to γ m n (r h , d). Then, ṽs can be calculated by [START_REF] Williams | Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography[END_REF] ṽs = 1 ιρc

N ∑ n=0 n ∑ m=-n h ′ n (kr s ) h n (kr h ) γ m n (r h , d)Y m n (θ, ϕ) = N ∑ n=0 n ∑ m=-n ψ m n (r s , d)Y m n (θ, ϕ), ( 10
)
where ρc is the characteristic impedance of the fluid, h ′ n is the first derivative of h n ,

ψ m n (r s , d) = 1 ιρc h ′ n (kr s ) h n (kr h ) γ m n (r h , d) ≡ 1 ιρc G n ( kr s , r h r s ) γ m n (r h , d) (11)
is the spherical wave spectrum of ṽs , and G n is the inverse propagator.

Figure 2 shows |G n | as a function of kr s for r h /r s = 2, 4, and 6. It can be noticed that |G n | increases significantly as kr s decreases, and n and r h /r s increase. Therefore, inverse sound source reconstruction is usually an ill-posed problem, i.e., small measurement errors on p h (and thus on γ m n (r h , d)) can be greatly amplified, so that large errors may occur on the reconstructed velocity. This can be dealt with by using a low-pass filter in the spherical harmonic domain [START_REF] Williams | Regularization methods for near-field acoustical holography[END_REF], which, however, reduces the spatial resolution of the reconstructed velocity field. In this work, it is assumed that p h is free from errors and that G n can be accurately evaluated in a digital computer. Thus, if d = d, one has ṽs = v s with no need of low-pass filtering. It is worth noting that, unlike measurement errors, the use of a nominal d ̸ = d leads to biased errors on the reconstructed field, so that they cannot be reduced by low-pass filtering. In the following, we discuss two cost functions that can be used to estimate d from p h , and thus to reduce such biased errors.

r h /r s = 6 n = 0 n = 1 n = 2 n = 3 n = 4 n = 0 n = 1 n = 2 n = 3 n = 4 n = 0 n = 1 n = 2 n = 3 n = 4

Model parameter estimation

As stated before, the number of spherical harmonics required to accurately describe p h depends on d. As |δd| = |d-d| increases, more terms must be retained in the spherical harmonic expansion [START_REF] Rafaely | Spatial alignment of acoustic sources based on spherical harmonics radiation analysis[END_REF][START_REF] Deboy | Acoustic center and orientation analysis of sound-radiation recorded with a surrounding spherical microphone array[END_REF]. Thus, the minimal point of a cost function that penalizes high-order harmonics can be a good estimation of d. This paper investigates the following two cost functions:

J 1 ( d) = ∑ N n=0 n ∑ n m=-n |ψ m n (r s , d)| 2 ∑ N n=0 ∑ n m=-n |ψ m n (r s , d)| 2 , ( 12
)
and

J 2 ( d) = 1 - ∑ N n=0 ∑ n m=-n |ψ m n (r s , d)| 2 ∑ N n=0 ∑ n m=-n |ψ m n (r s , d)| 2 . ( 13
)
Because of the factor n in the summation, J 1 increases significantly as the contribution of high-order spherical harmonics increases. Unlike J 1 , J 2 requires a cut-off order N < N , so that this function measures the contribution of the spherical harmonics of orders n > N ; note that 0 ≤ J 2 ≤ 1. J 2 would benefit from prior knowledge of the spatial resolution of the source velocity pattern, so that N could be easily determined. In this case, J 2 = 0 if d = d.

These cost functions were proposed by Rafaely [START_REF] Rafaely | Spatial alignment of acoustic sources based on spherical harmonics radiation analysis[END_REF] in order to evaluate the spectral changes in the pressure field due to translation of the coordinate system, and thus this author formulated them in terms of γ m n instead of ψ m n . Since we are interested in reconstructing the velocity field on a source sphere whose radius is known, it is more natural to use ψ m n for the purposes of this work.

Simulation results

In order to illustrate the ideas discussed throughout the paper, this section presents simulation results corresponding to a source sphere (Γ s ) of radius r s placed at d = (0, 0, 0.1r s ). We investigate four axisymmetric velocity fields (v s ) on Γ s with the following wave spectra (ψ m n = 0 if not specified):

(a) ψ 0 0 = 1;

(b) ψ 0 1 = 1; (c) ψ 0 0 = 1, ψ 0 1 = 0.3, ψ 0 2 = 0.1; (d) ψ 0 0 = 0.2, ψ 0 1 = 1, ψ 0 4 = 0.1.
Figure 3 shows these radial velocity fields, where the shape and the color map indicate, respectively, the amplitude and phase (in degrees) of v s . Notice that the vertical axis, z, is the symmetry axis.

- To obtain the "measured" pressure fields (p h ) on Γ h corresponding to the velocity fields shown in Fig. 3, we let ρc = 415kg.m -2 .s -1 and r h = 2r s . Two different kr s values are considered: kr s = 0.1 and kr s = 5.

For the inverse problem of determining ṽs from p h , we only deal with points Õ located on the z axis due to the axisymmetry of p h . Therefore, only positions given by d = (0, 0, z s ) are considered here. In addition, we let N = 11 in the spherical harmonics decomposition of p h (problem (9))we have observed that this value ensures that ||p(x h ) -p h || < 0.001||p h || for the whole d range considered, namely, -0.5r s ≤ z s ≤ 0.5r s . Since the spherical harmonics with m ̸ = 0 are not axisymmetric, only the harmonics with m = 0 are taken into account, which greatly reduces the computational cost because we deal with N + 1 instead of (N + 1) 2 functions. Finally, the usual inner product of P is approximately evaluated from 30 samples, which correspond to the co-latitudinal circumferences shown in Fig. 3. Area weighted factors are applied to the samples in order not to privilege densely sampled regions on the sphere.

Figure 4 shows the normalized least-square error (||ṽ s -v s ||/||v s ||) in velocity reconstruction, as well as the cost functions J 1 and J 2 (cut-off order N = 2) as functions of z s . These results correspond to the velocity fields shown in Fig. 3 for kr s = 0.1 and kr s = 5. The vertical dotted line indicates the real position of the source sphere.

The results reveal that even a small δd leads to a large error in the reconstructed velocity, so that d is a very important parameter of the radiation model. For kr s = 0.1, the global minimal point of J 1 occurs at: (a) z s = 0.1r s , (b) z s = 0.1r s , (c) z s = 0.34r s , and (d) z s = 0.1r s . Similarly, the global minimal point of J 2 occurs at: (a) z s = 0.1r s , (b) z s = 0.1r s , (c) z s = 0.1r s , and (d) z s = 0.09r s . Thus, J 1 and J 2 give the correct source position for the cases (a) and (b). For the case (c), unlike J 2 , 18 th International Congress on Sound and Vibration, 10-14 July 2011, Rio de Janeiro, Brazil J 1 fails significantly. For the case (d), J 1 is able to correctly identify the source position, whereas J 2 leads to a small error.

z s /r s J 2 (a) (b) (c) (d) (a) (b) (c) (d) (a) (b) (c) (d) kr s = 5
For kr s = 5, the global minimal point (or points) of J 1 occurs at: (a) z s = 0.1r s , (b) z s = (-0.21, 0.41)r s , (c) z s = 0.13r s , and (d) z s = 0.35r s . For J 2 , it occurs at z s = 0.1r s in the four cases under study. The source position is now easily identified for the case (a) by J 1 , but this cost function fails for the remaining cases. Nevertheless, for the case (c), the "high-frequency" (kr s = 5) simulations lead to an optimal position that is closer to the real position of the source sphere when compared to the "low-frequency" simulations (kr s = 0.1). For the case (b) (dipole source), the real position of the source corresponds to a local maximum point of J 1 , which is located between its two symmetrical global minima. This has also been observed by Deboy [START_REF] Deboy | Acoustic centering and rotational tracking in surrounding spherical microphone arrays[END_REF]. Unlike J 1 , J 2 presents a similar behavior at low and high frequencies. Note that J 2 (0.1r s ) = 0 for the cases (a), (b), and (c). This is due to the fact that the chosen cut-off order ( N = 2) is larger than or equal to the spatial resolution of the source, whereas J 2 (0.1r s ) > 0 for the case (d).

The results show that using J 1 to identify the source position might lead to very large errors in inverse source reconstruction. On the other hand, if N is properly set, the real position of the source is ensured to be a minimal point of J 2 . However, J 2 is less sensitive to z s than J 1 in the neighborhood of the real position. Hence, in practice, J 2 cannot be sufficiently robust to identify the source position.

Conclusion

This paper has addressed the issue of inverse sound source reconstruction by exterior SAH. In particular, the effects of the source position on the reconstructed velocity field have been studied. In addition, we have investigated two cost functions to estimate the position of the source sphere from the sound pressure it radiates on the hologram sphere.

It has been shown that using a nominal position with a small deviation from the real position leads to large errors in the reconstructed velocity field. Thus, the position of the source sphere inside the hologram sphere must be precisely known in order to produce meaningful results. However, the cost functions investigated here are not always able to accurately identify the source position. Hence, the discussion of source position identification based on spherical harmonic analysis is left open for future work.

Figure 1 :

 1 Figure 1: Illustration of the exterior SAH: (a) Two imaginary spheres Γ h and Γ s enclose a source that radiates under free-field conditions; (b) d known: The radial particle velocity (v s ) on Γ s is reconstructed from the sound pressure (p h ) on Γ h ; (c) d unknown: The radial particle velocity (ṽ s ) on a translated sphere Γs is reconstructed from p h and a nominal source position ( d).

Figure 2 :

 2 Figure 2: Inverse propagator up to order n = 4 as a function of kr s for r h /r s = 2, 4, and 6.

Figure 3 :

 3 Figure 3: Axisymmetric radial velocity fields (v s ) on Γ s . The shape and the color map indicate, respectively, the amplitude and phase (in degrees) of v s .

Figure 4 :

 4 Figure 4: Error in velocity reconstruction (||ṽ s -v s ||/||v s ||), and cost functions J 1 and J 2 (cut-off order N = 2) as functions of z s . Results corresponding to the real velocity fields shown in Fig. 3 for kr s = 0.1 and kr s = 5. The vertical dotted line indicates the real position of the source sphere.


Lars Andersson

AND Pieter Blue
email: p.blue@ed.ac.uk

Jean-Philippe Nicolas
email: jean-philippe.nicolas@univ-brest.fr

A DECAY ESTIMATE FOR A WAVE EQUATION WITH TRAPPING AND A COMPLEX POTENTIAL

In this brief note, we consider a wave equation that has both trapping and a complex potential. For this problem, we prove a uniform bound on the energy and a Morawetz (or integrated local energy decay) estimate. The equation is a model problem for certain scalar equations appearing in the Maxwell and linearised Einstein systems on the exterior of a rotating black hole.

Introduction

We consider the Cauchy problem

-∂ 2 t + ∂ 2 x + V (∆ ω -N ) + iǫW u = 0, (1)
u(0, x, ω) = ψ 0 (x, ω), ∂ t u(0, x, ω) = ψ 1 (x, ω),

on (t, x, ω) ∈ M = R × R × S 2 with smooth, compactly supported initial data. Here u is a complex function u = v + iw,

V = 1 x 2 + 1 ,
W is a smooth, real-valued, compactly supported function which is nonvanishing at x = 0 and uniformly bounded by 1, and ǫ > 0 is a small parameter. Finally, ∆ ω is the Laplacian in the angular variables and N is a number chosen to be sufficiently large to allow us to avoid certain technical problems. The equation (1) has both trapping, which occurs at x = 0, and a complex potential, which does not vanish at the trapped set. The interaction of these creates problems, which appear to frustrate the use of energy and Morawetz estimates at the classical level. By adapting known, pseudodifferential methods, we show how to overcome these problems. We now state our main result in terms of the energy

E(t) = 1 2 {t}×R×S 2 |∂ t u| 2 + |∂ x u| 2 + V (|∇ ω u| 2 + N |u| 2 ) dxd 2 ω, Theorem 1.
There is a constant C such that, if ψ 0 and ψ 1 are such that E(0) is finite then

∀t ∈ R : E(t) ≤ CE(0), (3a)
M |∂ x u| 2 x 2 + 1 + x 2 1 + x 2 |∇ ω u| 2 1 + |x| 3 + |∂ t u| 2 x 2 + 1 + |u| 2 1 + |x| 3 dtdxd 2 ω ≤ CE(0), (3b) M |u||∂ t u| 1 + |x| 3 dtdxd 2 ω ≤ CE(0). (3c)
Date: July 29, 2011 File:ModelProblem.tex.

Since the equation (1) has t independent coefficients, one might naively think that Noether's theorem provides a positive conserved energy. However, for the Lagrangian L 1 [u, ∂u] = -(∂ t u) 2 + (∂ x u) 2 + V (∇ ω u • ∇ ω u + N u 2 ) -iǫW u 2 , which has the wave equation ( 1) as its Euler-Lagrange equation, the conserved quantity associated to the time translation symmetry is indefinite, being approximately the energy of the real component of u minus the energy of the imaginary component (plus ǫ times a term involving W vw). On the other hand, a Lagrangian of the form

L 2 [u, ∂u] = -|∂ t u| 2 + |∂ x u| 2 + V (|∇ ω u| 2 + N |u| 2 )
, which corresponds to the energy expression considered above, does not yield equation ( 1) as its Euler-Lagrange equation.

The wave equation ( 1) is a model for equations arising in the study of the Maxwell and linearised Einstein equations outside a Kerr black hole. The Kerr black holes are a family of Lorentzian manifolds arising in general relativity, and they are characterised by a mass parameter M and an angular momentum parameter a. Black holes are believed to be the enormously massive objects at the center of most galaxies. The case |a| ≤ M is the physically relevant. The case a = 0 is the Schwarzschild class of black holes.

It is expected that every uncharged black hole will asymptotically approach a Kerr solution under the dynamics generated by the Einstein equations of general relativity. The wave, Maxwell, and linearised Einstein equations on a fixed Kerr geometry are a sequence of increasingly accurate models for these dynamics. By projecting on a null tetrad, the Maxwell and linearised Einstein fields can be decomposed into sets of complex scalars, the Newman-Penrose (NP) scalars [START_REF] Newman | An Approach to Gravitational Radiation by a Method of Spin Coefficients[END_REF][START_REF] Newman | Errata: An Approach to Gravitational Radiation by a Method of Spin Coefficients[END_REF][START_REF] Geroch | A space-time calculus based on pairs of null directions[END_REF]. It is well-known that the NP scalars with extreme spin weights satisfy decoupled wave equations, known as the Teukolsky equations, and that the solutions to these reduced equations can be used to reconstruct the full system [START_REF] Teukolsky | Rotating black holes: Separable wave equations for gravitational and electromagnetic perturbations[END_REF].

For the Maxwell field on the Kerr background, the spin weight 0 NP scalar can be treated in the same way, and the resulting equation is known as the Fackerell-Ipser equation [START_REF] Fackerell | Weak electromagnetic fields around a rotating black hole[END_REF]. For linearized gravity on the Schwarzschild background, it is also well known that the imaginary part of the spin weight 0 NP scalar is governed by a wave equation, the Regge-Wheeler equation [START_REF] Regge | Stability of a Schwarzschild singularity[END_REF][START_REF] Price | Nonspherical Perturbations of Relativistic Gravitational Collapse. II. Integer-Spin, Zero-Rest-Mass Fields[END_REF]. The corresponding equation for the real part is more complicated, cf. [START_REF] Zerilli | Effective Potential for Even-Parity Regge-Wheeler Gravitational Perturbation Equations[END_REF][START_REF] Moncrief | Gravitational perturbations of spherically symmetric systems. I. The exterior problem[END_REF], see also [START_REF] Aksteiner | Linearized gravity and gauge conditions[END_REF].

It was recently shown [START_REF] Aksteiner | Linearized gravity and gauge conditions[END_REF] that in the general (|a| < M ) Kerr case, by imposing a gauge condition related to the wave coordinates gauge, the equation for both the real and imaginary parts of the spin weight 0 NP scalar of the linearized gravitational field may be put in a form analogous to the Regge-Wheeler and Fackerell-Ipser equations. Explicitely (in the Kerr spacetime with signature -+++, working in Boyer-Lindquist coordinates) these take the form

∇ α ∇ α + 2s 2 M (r -ia cos θ) 3 u = 0, (4)
where s = 0 corresponds to the free scalar wave equation, s = 1 corresponds to the Maxwell (Fackerell-Ipser) case, while s = 2 corresponds to the linearized gravity (generalized Regge-Wheeler) case. In particular, for the a = 0 cases, the analogues of the Regge-Wheeler equations have complex potentials, with the imaginary part depending continuously on a.

For the wave equation in the Schwarzschild case, the use of the energy estimate [START_REF] Wald | Note on the stability of the Schwarzschild metric[END_REF], like (3a) with C = 1, and Morawetz estimates 1 are well established [START_REF] Laba | Global existence and scattering for the nonlinear Schrödinger equation on Schwarzschild manifolds[END_REF][START_REF] Blue | Semilinear wave equations on the Schwarzschild manifold. I. Local decay estimates[END_REF][START_REF] Blue | Errata for "Global existence and scattering for the nonlinear Schrodinger equation on Schwarzschild manifolds[END_REF][START_REF] Blue | Uniform decay of local energy and the semi-linear wave equation on Schwarzschild space[END_REF][START_REF] Dafermos | The red-shift effect and radiation decay on black hole spacetimes[END_REF]. In the Morawetz estimate (3b), there is a loss of control of time and angular derivatives near x = 0, in the sense that the integrand cannot control |x| p (|u| qt |∂ t u| 2-qt + |u| qω |∇ ω u| 2-qω ) with both p = 0 and either q t = 0 or q ω = 0. The presence of trapping makes some loss unavoidable [START_REF] Ralston | On the construction of quasimodes associated with stable periodic orbits[END_REF]. By applying "angular modulation" and "phase space analysis", the range for the angular parameter q ω can be refined to p = 0 and q > 0 [START_REF] Blue | Phase space analysis on some black hole manifolds[END_REF]. This type of refinement is crucial in the current paper. Alternatively, certain pseudodifferential operators have been used to obtain refinements near x = 0, to p > 0 and q t = q ω = 0 [START_REF] Marzuola | Strichartz estimates on Schwarzschild black hole backgrounds[END_REF].

For the wave equation in the general (|a| < M ) Kerr case, it is possible to apply Fourier transforms first in the φ and t variables 2 and then the remaining variables. The individual φ modes decay pointwise [START_REF] Finster | Decay of solutions of the wave equation in the Kerr geometry[END_REF]. Although the problem has a timetranslation symmetry, because the generator of time translations fails to be a timelike vector with respect to the Lorentzian inner product of the Kerr geometry, there is no positive, conserved energy. A major advance was the proof that, in the slowly rotating case |a| ≪ M , there is a uniform energy bound, like estimate (3a). The first proof used an estimate similar to (3b), but with additional restrictions on the support of the Fourier transform [START_REF] Dafermos | A proof of the uniform boundedness of solutions to the wave equation on slowly rotating Kerr backgrounds[END_REF]. Independent work [START_REF] Tataru | Local energy estimate on Kerr black hole backgrounds[END_REF] established estimates similar to (3a) and (3b), but with no restriction on the Fourier support, and there were subsequent pseudodifferential refinements [START_REF] Tohaneanu | Strichartz estimates on kerr black hole backgrounds[END_REF]. Also, the first two authors have proved similar results using methods which require two additional levels of regularity but which completely avoid the use of Fourier transforms. Morawetz estimates and refinements are a crucial step in proving pointwise decay estimates [START_REF] Blue | Phase space analysis on some black hole manifolds[END_REF][START_REF] Blue | Uniform decay of local energy and the semi-linear wave equation on Schwarzschild space[END_REF][START_REF] Dafermos | The red-shift effect and radiation decay on black hole spacetimes[END_REF][START_REF] Dafermos | The black hole stability problem for linear scalar perturbations[END_REF][START_REF] Luk | Improved decay for solutions to the linear wave equation on a Schwarzschild black hole[END_REF][START_REF] Luk | A vector field method approach to improved decay for solutions to the wave equation on a slowly rotating kerr black hole[END_REF] and Strichartz estimates [START_REF] Marzuola | Strichartz estimates on Schwarzschild black hole backgrounds[END_REF][START_REF] Tohaneanu | Strichartz estimates on kerr black hole backgrounds[END_REF], including the long-conjectured, inverse-cubic, Price law [START_REF] Metcalfe | Price's law on nonstationary spacetimes[END_REF][START_REF] Tataru | Local decay of waves on asymptotically flat stationary space-times[END_REF].

The study of the Maxwell and linearised Einstein systems in the Kerr geometry is still in its infancy. For the general Kerr case, a certain transformed, separated version of the Teukolsky system has no exponentially growing modes [START_REF] Whiting | Mode stability of the Kerr black hole[END_REF]. In the Schwarzschild case, the φ modes of the Teukolsky equation decay pointwise [START_REF] Finster | Decay of solutions of the Teukolsky equation for higher spin in the Schwarzschild geometry[END_REF]. Recently, improved decay estimates for the Regge-Wheeler type equation (4) on the Schwarzschild background, giving decay rates of t -3 , t -4 and t -6 , respectively, for s = 0, 1, 2, have been proved [START_REF] Donninger | On pointwise decay of linear waves on a Schwarzschild black hole background[END_REF]. The Regge-Wheeler equation has been used with the full system to prove energy and Morawetz (and pointwise decay) estimates for the Maxwell system [START_REF] Blue | Decay of the Maxwell field on the Schwarzschild manifold[END_REF] and, more recently, under assumptions on the asymptotic behaviour, for the full (not merely linearised) Einstein equation [START_REF] Holzegel | Ultimately schwarzschildean spacetimes and the black hole stability problem[END_REF].

For the spin-weight 0 equations arising from the Maxwell and linearised Einstein equations with a = 0, the presence of complex potentials in the reduced equations (4) prevents the existence of a positive conserved energy and means that an unrefined Morawetz estimate, such as (3b), is insufficient to control the growth of the energy, which is why we also prove estimate (3c). In the Kerr geometry, there is no positive, conserved energy because the generator of the time-translation symmetry fails to be timelike everywhere. In contrast, for the equation (1), there is no positive, conserved energy because the complex potential prevents any variational approach from providing an energy-momentum tensor that satisfies the dominant 1 These are also called integrated local energy decay estimates. 2 Here φ is the azimuthal angle, which would be one component of ω in the notation of this paper.

energy condition. The method of this paper combines a Fourier-transform-in-time technique (as in [START_REF] Dafermos | A proof of the uniform boundedness of solutions to the wave equation on slowly rotating Kerr backgrounds[END_REF][START_REF] Tataru | Local energy estimate on Kerr black hole backgrounds[END_REF]) with a "modulation" (or Fourier-rescaling) technique (from [START_REF] Blue | Phase space analysis on some black hole manifolds[END_REF]).

As is common, C will be used to denote a constant which may vary from line to line, but which is independent of the choice of u or T . The notation A B is used to denote that there is some C such that A < CB, with C independent of u and T , and similarly for .

A preliminary energy estimate

We derive an estimate for an energy for the wave equation ( 1) by integrating by parts against ∂ t ū and following the standard procedure for getting an energy estimate:

0 = Re (∂ t ū)(-∂ 2 t + ∂ 2 x + V (∆ ω -N ) + iǫW )u = - 1 2 ∂ t |∂ t u| 2 + ∂ x Re((∂ t ū)∂ x u) - 1 2 ∂ t |∂ x u| 2 + ∇ ω • Re((∂ t ū)∇ ω u) - 1 2 ∂ t V (|∇ ω u| 2 + N |u| 2 ) -ǫW Im((∂ t ū)u).
Introducing an energy which we denote by

E(t) = 1 2 {t}×R×S 2 |∂ t u| 2 + |∂ x u| 2 + V (|∇ ω u| 2 + N |u| 2 ) dxd 2 ω,
assuming that u decays sufficiently rapidly as |x| → ∞, and integrating the previous formula over a region [t 1 , t 2 ] × R × S 2 , we find

E(t 2 ) -E(t 1 ) = [t1,t2]×R×S 2 -ǫW Im((∂ t ū)u) dtdxd 2 ω. (5)
In particular, note that the energy fails to be conserved and that an estimate of the form (3b) would be insufficient to control the right-hand side. There is, however, a trivial exponential bound:

E(t 2
) ≤ e ǫ(t2-t1) E(t 1 ).

The Morawetz estimate

Following the standard procedure for investigating the wave equation, we derive a Morawetz estimate by multiplying the wave equation by (f (x)∂ x ū+q(x)ū), where f and q are real-valued functions.

In performing this calculation, it is useful to observe that

q ′ (x)Re(ū∂ x u) = ∂ x q ′ 2 ūu - 1 2 q ′′ ūu.
Using this and applying the product rule term-by-term, one finds

Re (f ∂ x ū + qū) -∂ 2 t u + ∂ 2 x u + V (∆ ω -N )u + iǫW u = ∂ t p t + ∂ x p x + ∇ ω p ω + - 1 2 f ′ + q |∂ t u| 2 - 1 2 f ′ + q |∂ x u| 2 + 1 2 f ′ -q V + 1 2 f (∂ x V ) |∇ ω u| 2 + N 1 2 f ′ -q V + 1 2 f (∂ x V ) + 1 2 q ′′ |u| 2 -ǫf W Im((∂ x ū)u), (6)
where

p t = p t (f, q; u) = -Re((f (∂ x ū) + qū)(∂ t u)), p x = p x (f, q, u) = 1 2 f |∂ t u| 2 + 1 2 f |∂ x u| 2 - 1 2 f V |∇ ω u| 2 + qRe(ū∂ x u) - 1 2 (N f V + q ′ )|u| 2 , p ω = p ω (f, q; u) = f V Re((∂ x ū)(∇ ω u)) + qV Re(ū∇ ω u).
We take f = -arctan(x), for which

f ′ = -(x 2 + 1) -1 = -V , f ′′ = 2x(x 2 + 1) -2
, and f ′′′ = -2(3x 2 -1)(x 2 + 1) -3 . We take q = f ′ /2 + δ(1 + x 2 ) -1 arctan(x) 2 for some sufficiently small δ.

We use the notation

E f ∂x+q (t) = {t}×R×S 2 Re(f (∂ x ū)∂ t u) + Re(qū(∂ t u)) dxd 2 ω,
and observe that, by a simple Cauchy-Schwarz argument, there is the estimate

f ∂x+q | ≤ CE.
Observing that the left-hand side of ( 6) vanishes, we have

0 = ∂ t p t + ∂ x p x + ∇ ω p ω + δ arctan(x) 2 1 + x 2 |∂ t u| 2 + 1 1 + x 2 (1 -δ arctan(x) 2 )|∂ x u| 2 + x arctan(x) -δ arctan(x) 2 (1 + x 2 ) 2 |∇ ω u| 2 + N x arctan(x) -δ arctan(x) 2 (1 + x 2 ) 2 + 1 2 q ′′ |u| 2 -ǫf W Im((∂ x ū)u)
Taking ǫ sufficiently small, N sufficiently large, and δ sufficiently small, the factors in front of |∂ x u| 2 and |u| 2 are nonnegative and one can dominate the term involving W using these two terms. (These estimates are uniform, in the sense that, if the estimate holds for choices of ǫ 0 , N 0 , and δ 0 , then it remains valid for ǫ < ǫ 0 , N = N 0 , and δ = δ 0 .) Thus, by integrating over a time-space slab

M [t1,t2] = [t 1 , t 2 ] × R × S 2 ,
one can conclude that there is a constant C such that Since χ 1 is smooth, there is a uniform bound on its derivative and second derivative, each of which are supported on [0, 1] ∪ [T, T + 1], so that there is a constant

E(t 2 ) + E(t 1 ) M [t 1 ,t 2 ] |∂ x u| 2 x 2 + 1 + | arctan(x)| 2 |∇ ω u| 2 1 + |x| 3 + |∂ t u| 2 x 2 + 1 + |u| 2 1 + |x| 3
C such that |∂ t χ 1 | + |∂ 2 t χ 1 | ≤ Cχ 2 . The functions u 1 = χ 1 χ |x|≤2 u, u 2 = χ 2 χ |x|≤2 u, u 3 = χ 1 u, satisfies the equation -∂ 2 t + ∂ 2 x + V (∆ ω -N ) + iǫW u 1 = F (u 2 , ∇u 2 , t, x) + G(u 3 , ∇u 3 , t, x), (8) where
F (u 2 , ∇u 2 , t, x) = -2(∂ t χ 1 )(∂ t u 2 ) -(∂ 2 t χ 1 )u 2 G(u 3 , ∇u 3 , t, x) = 2(∂ x χ |x|≤2 )(∂ x u 3 ) + (∂ 2
x χ |x|≤2 )u 3 . Since all functions of t in this equation are smooth and supported in t ∈ [-2, T + 2], they are Schwartz class in t, so we may take the Fourier transform in t and remain in the Schwartz class. We will use to denote the Fourier transform in t, and τ for the argument of such functions. We will typically use the word "functions" to describe u, u 1 , u 2 , and u 3 and the words "Fourier transforms" to describe their Fourier transforms. We will use L 2 to denote L 2 (dωdxdt) for functions and to denote L 2 (dωdxdτ ) for Fourier transforms. We will use • for • L 2 unless otherwise specified.

We introduce the following space-time integrals

I(T ) = T +2 -2 2 -2 S 2 x 2 |∂ t u| 2 + |∂ x u| 2 + |u| 2 dωdxdt, J(T ) = R×R×S 2 |τ | 6/5 | u 1 | 2 dωdxdτ.
The dependence of J upon T is through the smooth cut-off χ 1 in u 1 . Typically, the argument T will be clear from context and will be omitted. From the Morawetz estimate [START_REF] Dafermos | A proof of the uniform boundedness of solutions to the wave equation on slowly rotating Kerr backgrounds[END_REF] and the exponential bound on the energy, it follows that I E(T ) + E(0). We now aim to prove a Morawetz estimate using the Fourier transform. We take

f = -arctan(|τ | α x), q = f 2 = 1 2 |τ | α 1 + |τ | 2α x 2 ,
with α ∈ [0, 1/2]. We multiply the Fourier transform of equation ( 8) by (f ∂ x +q) ¯ u 1 , and integrate the real part over R × R × S 2 . This integral is convergent because all the functions are compactly supported in time, so the Fourier transforms are Schwartz class.

4.2.

Controlling the terms arising from the cut-off. We consider first the integral arising from the right-hand side of (8). This is

R×R×S Re (f ∂ x + q) ¯ u 1 F + G dωdxdτ ≤ (f ∂ x + q) ¯ u 1 F + G .
The terms on the right can be estimated by

(f ∂ x + q) ¯ u 1 ≤ f ∂ x u 1 + q u 1 , f ∂ x u 1 ∂ x u 1 I 1/2 , q u 1 |τ | α u 1 u 1 + |τ | 1/2 u 1 I 1/2 + J 1/2 ,
and

F + G ≤ F + G .
Because G is supported only for t ∈ [-1, T + 1] and x ∈ [-2, 2], we have

G I 1/2 .
Similarly, because F is supported only for

t ∈ [-1, 0] ∪ [T, T + 1] and x ∈ [-2, 2],
we have that at each instant in t, the function F is bounded in L 2 (dxdω) by either CE(0) 1/2 or CE(T ) 1/2 . Since we are considering two intervals in t of length 1, we have

F = F C(E(0) 1/2 + E(T ) 1/2 ).
Thus, the terms on the right-hand side of the Fourier transform of (8) are bounded by

R×R×S 2 Re (f ∂ x + q) ¯ u 1 F + G dωdxdτ ≤ C(E(0) 1/2 + E(T ) 1/2 + J 1/2 )(E(0) 1/2 + E(T ) 1/2 ). (9)
4.3. The Morawetz estimate for the approximate solution. If we multiply the left-hand side of the Fourier transform of the wave equation ( 8) by (f ∂ x + q) ¯ u 1 and take the real part, then we have the analogue of ( 6)

Re f ∂ ¯ u 1 + q ¯ u 1 τ 2 u 1 + ∂ x u 1 + V (∆ ω -N ) u 1 + iǫW u 1 = ∂ x p x + ∇ ω p ω + - 1 2 f ′ + q |τ u 1 | 2 - 1 2 f ′ + q |∂ x u 1 | 2 + 1 2 f ′ -q V + 1 2 f (∂ x V ) |∇ ω u 1 | 2 + N 1 2 f ′ -q V + 1 2 f (∂ x V ) + 1 2 q ′′ | u 1 | 2 -ǫf W Im((∂ x ¯ u 1 ) u 1 ), (10)
where

p x = p x (f, q, u 1 ) = 1 2 f |∂ t u 1 | 2 + 1 2 f |∂ x u 1 | 2 - 1 2 f V |∇ ω u 1 | 2 + qRe( ¯ u 1 ∂ x u 1 ) - 1 2 (N f V + q ′ )| u 1 | 2 , p ω = p ω (f, q; u 1 ) = f V Re((∂ x ¯ u 1 )(∇ ω u 1 )) + qV Re( ¯ u 1 ∇ ω u 1 ).
Note that there is no p t term because, for Fourier transforms, the analogue of the product rule is simply -iτ u 1 u 1 = ¯ u 1 iτ u 1 . When this equality is integrated over a space-time slab, the p x and p ω terms integrate to zero, and the remaining terms are all non-negative except for those arising from q ′′ and from W . The integral of the term involving W is bounded by I.

We now consider the term involving q ′′ :

1 2 q ′′ | u 1 | 2 = |τ | 3α 1 -3|τ | 2α x 2 (1 + |τ | 2α x 2 ) 3 | u 1 | 2 .
From the positivity of the remaining terms and the bound (9) of the terms coming from the right-hand side of the wave equation ( 8) for u 1 , we have that

R×R×S 2 |τ | 3α 1 -3|τ | 2α x 2 (1 + |τ | 2α x 2 ) 3 | u 1 | 2 dωdxdτ ≤ C(E(0) 1/2 + E(T ) 1/2 + J 1/2 )(E(0) 1/2 + E(T ) 1/2 ).
This can be combined with an additional factor of M I, where M is a large constant.

(700 is sufficient.) The integral I dominates the integral of (|τ | 2 + 1)x 2 | u 1 | 2 and is bounded by C(E(T ) + E(0)). Thus, we have

R×R×S 2 |τ | 3α 1 -3|τ | 2α x 2 (1 + |τ | 2α x 2 ) 3 + M (|τ | 2 + 1)x 2 | u 1 | 2 dωdxdτ ≤C(E(0) 1/2 + E(T ) 1/2 + J 1/2 )(E(0) 1/2 + E(T ) 1/2
|τ | 3α 1 -3|τ | 2α x 2 (1 + |τ | 2α x 2 ) 3 + M | 2 + 1)x 2 χ |x|≤2 ≥ C|τ | 3α χ |x|≤2
and therefore we find The integrand is now controlled by I + J, which, from estimates [START_REF] Dafermos | A proof of the uniform boundedness of solutions to the wave equation on slowly rotating Kerr backgrounds[END_REF] and [START_REF] Fackerell | Weak electromagnetic fields around a rotating black hole[END_REF], we know can be estimated also by the sum of the initial and final energies. This leaves the estimate E(T ) -E(0) ≤ Cǫ(E(T ) + E(0)).

C(E(0) 1/2 + E(T ) 1/2 + J 1/2 )(E(0) 1/2 + E(T ) 1/2 ) ≥ R×R×S 2 |τ | 3α 1 -3|τ | 2α x 2 (1 + |τ | 2α x 2 ) 3 + M (|τ | 2 + 1)x 2 | u 1 | 2 dωdxdτ ≥ J, which implies J ≤C(E(T ) + E(0)). (11
By taking ǫ sufficiently small relative to the constant, we obtain a uniform bound on the energy E(T ) ≤ CE(0).

We note that, since all the constants were independent of T , the estimate holds uniformly in T . This proves the first statement (3a), in theorem 1. Combining this with estimate (7) (and estimating x 2 (1 + x 2 ) arctan(x) 2 ) gives the second, (3b). Finally, the arguments of this section and the bound on I + J, from estimates [START_REF] Dafermos | A proof of the uniform boundedness of solutions to the wave equation on slowly rotating Kerr backgrounds[END_REF] and [START_REF] Fackerell | Weak electromagnetic fields around a rotating black hole[END_REF], give the third result (3c).

Remark 2. Using the method given in [START_REF] Blue | Phase space analysis on some black hole manifolds[END_REF], the stronger Morawetz estimate [START_REF] Fackerell | Weak electromagnetic fields around a rotating black hole[END_REF] can be improved to control the integral of |τ | 2-ε | u| 2 for any ε > 0. Because of the presence of trapping, it is not possible to improve this to |τ | 2 | u| 2 .

1 .

 1 dtdxd 2 ω. (7) 4. Pseudodifferential refinements 4.The wave equation for an approximate solution. We define a smooth characteristic function of an interval [a, b] to be a function which is identically 1 on [a, b], which is supported on [a -1, b + 1], and which is monotonic on each of the [a -1, a] and [b, b + 1]. A smooth characteristic function of a collection of intervals, each of which are separated by distance at least two, is defined to be the sum of the smooth characteristic functions of each interval. Let T > 0 be a large constant. (Here large means larger than -log |ǫ| and 2.) Let χ 1 be a smooth characteristic function on [0, T ], and let χ 2 be a smooth characteristic function of [-1, 0] ∪ [T, T + 1]. Let χ |x|≤2 be a smooth characteristic function of [-1, 1]. We will use χ 1 , χ 2 , and χ |x|≤2 to denote χ 1 (t), χ 2 (t), and χ |x|≤2 (x) respectively.

) 4 . 4 . 2 W 2 W

 4422 Closing the energy estimate. It is now possible to estimate the integral on the right-hand side of the energy estimate[START_REF] Blue | Phase space analysis on some black hole manifolds[END_REF]. For |x| ≥ 1, the right-hand side (using the compact support of W and the Cauchy-Schwarz estimate) can be dominated by I ≤ C(E(T ) + E(0)). For |x| ≤ 1, we would like to dominate the integral over R × R × S 2 of |W Im( ū1 ∂ t u 1 )| by the integral J. However, this is not entirely correct, because in J there is a contribution arising from the support of u in the region t ∈ [-1, 0] ∪ [T, T + 1]. The error in this approximation is bounded by C(E(T ) + E(0)). Thus, we haveE(T ) -E(0) ≤ Cǫ E(T ) + E(0) + R×R×S Im( ū1 ∂ t u 1 )dωdxdt .We can also take the Fourier transform to obtain an estimate byCǫ E(T ) + E(0) + R×R×S Im( ¯ u 1 ∂ t u 1 )dωdxdt ≤ Cǫ E(T ) + E(0) + R×R×S 2 |τ || u 1 | 2 dωdxdt .

  ). By considering the two cases |τ | α |x| < M -1/2 and |τ | α |x| ≥ M -1/2 , one can see that if 2 -2α = 3α (i.e. α = 2/5), then