text:
# Dynamical geometry for multiscale dissipative particle dynamics

## 1 Introduction

The non-equilibrium behaviour of complex fluids continues to present a major challenge for both theory and numerical simulation. Such fluids include multi-phase flows, particulate and colloidal suspensions, polymers, amphiphilic fluids, including emulsions and microemulsions, and other fluids where Brownian motion is important. Over the last decade several strategies have been introduced both from a microdynamical point of view and from a continuum or macro-dynamical point of view. This article is particularly focused on mesoscopic/particulate approaches, but the Voronoi tessellation has been employed for similar applications to model colloidal particles and polymers by solving conventional Navier-Stokes like equations with moving boundary conditions (see  and references therein).

A recent contribution to the family of bottom-up approaches is the dissipative particle dynamics (DPD) method introduced by Hoogerbrugge and Koelman in 1992 . Successful applications of the technique have been made to colloidal suspensions , polymer solutions  and binary immiscible fluids . Dissipative particle dynamics has been shown to produce the correct macroscopic (continuum) theory; that is, for a one-component DPD fluid, the Navier-Stokes equations emerge in the large scale limit, and the fluid viscosity can be computed . However, even though dissipative particles have generally been viewed as clusters of molecules, the cluster of molecules are not considered as a thermodynamic system itself. In a following development of the DPD model, Flekkøy and Coveney  and Serrano and Español  derived a new model which provides a precise definition of the term ‘mesoscale’ and enables us to retrieve numerically the specific thermodynamics of the solvent in the sense that the exact Gibbs equilibrium is approached, the entropy is an increasing function in time, and the total energy is exactly conserved . In this approach, the fluid dissipative particles are defined as cells on a Voronoi lattice with variable sizes and masses; the Voronoi tessellation must be updated at each time step under the dynamics of the model. The scheme enables one to select any desired local scale for the mesoscopic description of a given problem. Indeed, the method may be used to deal with situations in which several different length scales are simultaneously present.

In this article, we discuss the computational aspects of the model and its computational overhead due to the Voronoi tessellation; in particular, we present the development of an efficient C++ code for the implementation of the new fluid particle model. Our discussion concentrates on the challenging aspects of maintaining the Voronoi tessellation over time, with various kinds of boundary conditions.

## 2 Multiscale dissipative particle dynamics

In our case, we follow the model described in detail in , but we apply its isothermal version, which has the advantage of greater simplicity. In fact, if heat diffusion is not important for the system under consideration, then the isothermal model furnishes an efficient alternative, because several simplifications are possible in the equations (for instance, the Itô-Stratonovich terms all vanish). Defining the state of a fluid particle by its position, mass and momentum $`𝐱=\{𝐫_i,M_i,𝐏_i\}`$, the resulting equations for the isothermal model are

$`d𝐫_i`$ $`=`$ $`𝐯_idt,`$
$`dM_i`$ $`=`$ $`{\displaystyle \underset{j}{}}{\displaystyle \frac{\rho _i+\rho _j}{2}}𝐂_{ij}𝐯_{ij}dt,`$
$`d𝐏_i`$ $`=`$ $`{\displaystyle \underset{j}{}}{\displaystyle \frac{A_{ij}}{2}}\left(p_{ji}\mathrm{𝟏}+𝚷_i+\mathrm{\Pi }_i\mathrm{𝟏}\right)𝐞_{ij}dt+d\stackrel{~}{𝐏}_i`$ (1)
$`+`$ $`{\displaystyle \underset{j}{}}{\displaystyle \frac{\rho _i+\rho _j}{2}}{\displaystyle \frac{𝐯_i+𝐯_j}{2}}𝐂_{ij}𝐯_{ij}dt.`$
where $`𝐯_i=𝐏_i/M_i`$ is the velocity, $`\rho _i=M_i/V_i`$ is the mass density, and $`d`$ is the spatial dimension. For a general variable $`x`$, we write $`x_{ij}=x_ix_j`$. The pressure $`p_i`$ is given through the equilibrium equations of state once the temperature $`T`$ is assigned. We have also introduced geometric quantities arising from the Voronoi construction: $`A_{ij}`$ is the area (length in 2D) of the face between cells $`(i,j)`$, $`𝐞_{ij}=(𝐫_i𝐫_j)/r_{ij}`$ with $`r_{ij}=|𝐫_i𝐫_j|`$ is the unit vector normal to the face $`(i,j)`$ and, finally, $`𝐂_{ij}=\frac{A_{ij}}{r_{ij}}\left(𝐀_{cm}\frac{𝐫_i+𝐫_j}{2}\right)`$ is a geometrical vector parallel to the face $`(i,j)`$. The dissipative stress tensor is

$$𝚷_i^{\alpha \beta }=\frac{\eta _i}{V_i}\left[\frac{1}{2}\underset{j}{}A_{ij}[𝐞_{ij}^\alpha 𝐯_j^\beta +𝐞_{ij}^\beta 𝐯_j^\alpha ]\frac{1}{D}\delta ^{\alpha \beta }\underset{j}{}A_{ij}𝐞_{ij}𝐯_j\right],$$
(2)
while the fluctuations in momentum are

$$d\stackrel{~}{𝐏}_i=\underset{j}{}\frac{1}{2}A_{ij}𝐝\stackrel{~}{\sigma }_j𝐞_{ij},$$
(3)
where the stress tensor is

$$𝐝\stackrel{~}{\sigma }_i=\left(4k_BT_i\frac{\eta _i}{V_i}\right)^{\frac{1}{2}}d\overline{𝐖}_i^S+\left(2Dk_BT_i\frac{\xi _i}{V_i}\right)^{\frac{1}{2}}tr[d𝐖_i].$$
(4)
Here, $`d𝐖_i`$ is a matrix $`(d\times d)`$ of independent Wiener increments, i.e. $`dW_i^{(\alpha ,\beta )}`$ has zero mean and variance $`dt`$ and is uncorrelated to any other random variables, $`d\overline{𝐖}_i^S=\frac{1}{2}\left(d𝐖_i+d𝐖_i^t\right)\frac{1}{d}tr\left[d𝐖_i\right]\mathrm{𝟏}`$ is a symmetric traceless matrix.

The SDEs (1) can be integrated numerically using any stochastic integrator, such as Euler, Heun or higher order schemes . For our simulations, we used an integrator based on a generalisation of the Trotter formula to the stochastic case  which proved to be simple to apply and very efficient.

## 3 Voronoi tessellation

The model is centrally based on the Voronoi tessellation. The Voronoi tessellation is simply defined by a set of points $`P=\{p_1,\mathrm{},p_N\}`$ in Euclidean space; a partition of the space assigning every point to its nearest site is called a Voronoi tessellation. The Voronoi cell $`V(p_i)`$ consists of all the points at least as close to $`p_i`$ as to any other site

$$V(p_i)=\{x:\left|p_ix\right|\left|p_jx\right|,ji\}.$$

The Voronoi tessellation can be computed through the Delaunay triangulation, which is the triangulation that maximises the minimum internal angle among the triangles . In terms of the fluid dissipative particles, the Delaunay triangulation indicates the neighbourhood of each fluid particle in the tessellation and the Voronoi tessellation indicates its spatial extension (see Figure 2). By first triangulating the set of points $`P`$ with the Delaunay triangulation, the Voronoi tessellation can then be constructed by using the Delaunay-Voronoi duality , i.e. to each Delaunay vertex corresponds a Voronoi region and to each Delaunay edge corresponds a Voronoi boundary surface or line in two dimensions.

The physical model requires maintenance of the Voronoi tessellation for each time step following the dynamics in the Lagrangian frame of reference and specification of the boundary conditions for the tessellation. In our opinion, direct coding of the construction of the Voronoi tessellation is not desirable. This still is in its own right a research field in computational geometry in order to obtain a stable and efficient triangulation. It involves writing exact number types to check the validity of predicates on points (for instance the insphere property for four points) and combinatorial aspects to handle the data structure. Nevertheless, writing an inefficient and unstable code to construct the Voronoi tessellation is quite simple, but the result will likely be orders of magnitude slower and unstable in most practical situations. Whilst there are several public libraries  that are able to compute the Voronoi-Delaunay tessellation for a given set of points. However, the algorithmic complexity of the construction is at best $`N\mathrm{log}N`$, where $`N`$ is the number of points, and thus it is crucial to use an efficient library to match that complexity. We were unable to find any library able to dynamically change the tessellation and impose periodic boundary conditions.

Previously, in our research group, the Triangle library  was employed for the construction of the Voronoi tessellation . Although the Triangle library does not support three dimensional Voronoi tessellations, it is the fastest library known to us in two dimensions. *Triangle* is based on the divide-and-conquer algorithm , which allows one to compute the tessellation given the complete list of points. It furnishes as well an incremental algorithm, which provides the functions to construct the tessellation point by point.

With *Triangle* the Voronoi tessellation was maintained over time by deleting the set of points $`P`$ from the tessellation, determining the new positions of the points and then computing the tessellation for the new points, for every time-step. Given the algorithmic complexity of the construction of Voronoi tessellations, this is clearly not an efficient way to maintain the tessellation over time; it is desirable to be able to locally restore the tessellation at each time-step, thus avoiding deleting and then re-computing the entire tessellation.

Subsequent awareness of the Computational Geometry Algorithms Library (CGAL)  gave us a homogeneous environment for two and three dimensional tessellations. CGAL incorporates a wide-range of function modifiers, including the insertion or removal of specific points in a tessellation, thus offering a highly flexible code for the construction and maintenance of a Voronoi tessellation over time. Consequently, we are currently employing the CGAL library for the development of the Voronoi DPD code in two and three spatial dimensions.

## 4 Boundary conditions

In this section, we describe possible implementations of periodic and wall boundary conditions in two and three dimensions based on the insertion-removal of image cells (points) around the simulation box. The insertion of polymers and colloidal particles is described using boundary surfaces. We present these algorithms in two dimensions because our three dimensional implementation is a straightforward application of the same techniques. The library takes care of the extra complexity of handling the three dimensional tessellation. Differences are pointed out in Section 4.4.

### 4.1 Periodic boundary conditions

Boundary conditions can be implemented by inserting ‘image points’ around the boundaries of the space domain. We call image points the points of the tessellation which lie outside the simulation box and represent images of centres of fluid particles inside the simulation box. The simplest way to achieve a periodic tessellation is to copy the simulation region to the eight surrounding boundary regions (see Figure 3).

However, this requires that the tessellation has to be computed for a total of $`9N`$ points if $`N`$ is the number of points to be triangulated in the simulation region. It is obviously not necessary to triangulate this number of image points to maintain the periodic boundary conditions and in light of the algorithmic complexity of the construction of the tessellation, it is crucial to minimise the number of image points. Clearly, it is only the points that lie ‘near’ a boundary face that must be copied to the boundary regions. To isolate these ‘boundary points’ we use the following geometrical criterion: a point is in the neighbourhood of the boundary if the presence of another point outside the domain can affect its Voronoi region. This criterion can be implemented with a simple algorithm. This is given in Algorithm 1.

In the first step of the algorithm, $`N`$ points are triangulated in the simulation region. In step $`2`$, the following geometrical properties of the Delaunay-Voronoi tessellation  are employed:

* Each node $`p_i`$ of the Delaunay triangulation corresponds to a Voronoi region $`V(p_i)`$.
* If $`v`$ is a Voronoi vertex at the junction of $`V(p_1),V(p_2),V(p_3)`$ then $`v`$ is the centre of the circle $`C(v)`$ determined by the points $`p_1,p_2`$ and $`p_3`$.
* $`C(v)`$ is the circumcircle for the Delaunay triangle $`(p_1,p_2,p_3)`$ containing no other points $`p_i`$.

It follows from the Delaunay-Voronoi properties above that for a Delaunay triangle $`(p_1,p_2,p_3)`$, if the circumcircle $`C(v)`$ does not intersect a boundary face of the domain, then no image point can affect the Voronoi regions $`V(p_1)`$, $`V(p_2)`$ and $`V(p_3)`$ and hence, according to our geometrical criterion, $`p_1,p_2`$ and $`p_3`$ do not have to be copied to the boundary regions. However, if $`C(v)`$ does intersect a boundary face, then an image point can affect the Voronoi regions $`V(p_1)`$, $`V(p_2)`$ and $`V(p_3)`$. Unfortunately, it is not possible to know if any boundary points will produce an image point within $`C(v)`$ until the tessellation for these boundary points is constructed. Nevertheless, a good upper-bound can be given on the number of points that must be copied by considering the worst situation that for any Delaunay triangle $`(p_1,p_2,p_3)`$ whose circumcircle $`C(v)`$ interests a boundary face, an image point will always be located within $`C(v)`$. In this case $`p_1,p_2`$ and $`p_3`$ must be copied to the boundary region of the opposite face; if $`C(v)`$ intersects the corner (i.e. two boundary faces), $`p_1,p_2`$ and $`p_3`$ must be copied to the boundary regions of the two opposite faces and also to the boundary region of the opposite corner (see Figure 3).

The points which belong to the convex hull of the triangulation have their Voronoi centre $`v`$ positioned at infinity and thus it is not possible to compute $`C(v)`$ in these cases. The solution we applied is to insert the corner points of the periodic boundary regions in order to create an external convex hull. After the image points are calculated the corners points are removed. This is not needed when the periodic boundary condition is already set and we want to maintain it, because the image points will form the convex hull.

The periodic boundary tessellation produced by Algorithm 1 is shown in Figure 4. Thus Algorithm 1 avoids copying all the points and, in assuming the worst case scenario that an image point is always located within $`C(v)`$ for any boundary point, this algorithm reduces the number of image points to the minimum upper-bound possible.

### 4.2 Fixed (solid) wall boundary conditions

A possible way to implement fixed (solid) wall boundary conditions for the Voronoi tessellation is to introduce image points into the boundary regions which correspond to the reflection of the boundary points in their boundary face(s) (see Figure 5). The interaction between a boundary particle and the ‘wall’ is then handled by the interaction between the particle and its ‘mirror-image’ particle. We have successfully implemented fixed wall boundary conditions for the Voronoi tessellation with Algorithm 2.

The algorithm follows the same principle as Algorithm 1; we first insert the corner points of the simulation domain, and then triangulate our $`N`$ points, therefore enabling us to compute the circumcircle of every Delaunay triangle in the simulation domain. Next in step 3, we isolate all of the boundary points. However, the criterion for a boundary point for fixed wall boundary conditions is more complicated than that for periodic boundary conditions. For any circumcircle $`C(V)`$ which intersects a boundary face, two additional requirements must be satisfied before $`p_1,p_2`$ or $`p_3`$ can be isolated as boundary points.

Firstly, for a given Delaunay triangle $`(p_1,p_2,p_3)`$, if $`C(V)`$ intersects a boundary face, it is possible for $`(p_1,p_2,p_3)`$ to lie on the wrong side of the boundary. For example, consider the fixed wall boundary conditions for the ‘C’ shaped domain in Figure 5. Suppose we are isolating the particles that lie in the neighbourhood of the first ‘shorter’ horizontal boundary face from the bottom in Figure 5. It is possible for the circumcircle of a Delaunay triangle $`(p_1,p_2,p_3)`$ lying above this boundary face to intersect this boundary and yet clearly $`(p_1,p_2,p_3)`$ cannot physically lie in the neighbourhood of this boundary. Thus, if $`C(V)`$ intersects a boundary face, we must then check that $`p_1,p_2`$ and $`p_3`$ lie on the correct side of the boundary; if we are traversing the boundary faces in a clockwise direction then $`p_1,p_2`$ and $`p_3`$ are on the correct side of the boundary face if they lie to the right-hand side of the boundary. Similarly, an anti-clockwise traversal of the boundary faces requires that $`p_1,p_2`$ and $`p_3`$ lie to the left-hand side of the boundary.

If these first two requirements are met, then we must check one final geometrical property. We refer again to the first ‘shorter’ horizontal boundary face up from the bottom of the ‘C’ shaped domain in Figure 5. Suppose for a Delaunay triangle $`(p_1,p_2,p_3)`$, the circumcircle $`C(V)`$ intersects this boundary face and $`p_1,p_2`$ and $`p_3`$ lie on the correct side of the boundary. Consider the supporting line of this boundary face (that is, the infinite line which is parallel to the boundary face segment and passes through it) and the perpendicular projection of the Delaunay vertices $`p_1,p_2`$ and $`p_3`$ onto this supporting line. For the majority of cases, these three projections will also lie on the boundary face segment in which case $`p_1,p_2`$ and $`p_3`$ can be reflected in the boundary face and are thus boundary points. However, for Delaunay triangles whose circumcircle intersects the left-end of this boundary face, it is possible that some or even all of these projections may not lie on the boundary face segment, that is, they do not have a reflection in this boundary face and hence are not boundary points. Thus, if $`C(V)`$ intersects a boundary face and $`p_1,p_2`$ and $`p_3`$ lie on the correct side of the boundary, we must then determine if the perpendicular projection of each Delaunay vertex onto the supporting line is also on the boundary face segment. For those that are, they are boundary points.

Having isolated the boundary points with the above criterion, we then insert the reflection of these points in their boundary face(s). We note that if $`C(V)`$ intersects a corner (i.e. two boundary faces) then the two additional geometrical requirements detailed above must be applied separately to each boundary. Finally in step 4 of Algorithm 2, we remove the corner points from the simulation.

This algorithm can be applied to any two dimensional polygon that can be described by a set of vectors. We show in Figure 5 the fixed wall boundary tessellation for a ‘C’ shaped domain, obtained using Algorithm 2. We note that for corners where the angle between the two boundary faces in the fluid is greater than $`180^0`$, the actual corner is not produced exactly. This is because near the corner, the insertion of an image point above one of the boundary faces will affect the Voronoi region of the other boundary face (see Figure 5). However, for realistic simulations involving tens of thousands of points, the ‘approximate’ corner obtained is perfectly sufficient.

### 4.3 Boundary surfaces

One final aim of this research is to add specificity to the model by inserting polymers and colloidal particles in the fluctuating hydrodynamic solvent based on the Voronoi dissipative particle dynamics model. Then, the rheological properties of the resulting system can be computed via direct simulations of shear flows.

Polymers and colloidal particles are inserted in the Voronoi tessellation using surface boundaries formed by chains of Voronoi centres. This is represented in Figures 6 and 7. The boundary surface are fixed giving a couple of Voronoi centres sufficiently close to form a Voronoi surface, which represents the boundary. How close they need to be is computed from the typical scale of the system or the smallest Voronoi fluid particle.

The physical meaning of the two Voronoi centres depends on their position. The Voronoi cell internal to the polymer or the colloid is a coarse-grained representation of a part of the polymer. The Voronoi centre in the fluid is a fluid particle which does not slip on the colloidal/polymer surface.

This description of surface boundary could be used as well for wall boundaries, with the disadvantage that an higher number of cells is needed to construct the boundary compared to the approach presented in Section 4.2, but with the advantage of an easy construction and maintenance of the boundary during the dynamics. In fact, once the two cells forming the boundary surface are kept at a fixed distance, then the maintenance of the tessellation does not involve any particular extra work.

### 4.4 Efficiency of image point boundary conditions

The efficiency of the construction of periodic/solid boundary conditions using image points depends on the number of images that we need to insert. Intuitively, this is a limitation when we tessellate a small number of points, because in this case the surface:volume ratio is high. We quantitatively measured the number of images needed to impose boundary conditions in two and three spatial dimensions.

The results are reported in Table 1. First, we note that while in two dimensions we have 8 neighbour regions around the box, in three dimensions there are 26. The effect of this geometrical change is evident. The number of shell points represents the number of points that have at least one image. The number of fluid particles which are in the core, and therefore not affected by boundary conditions, is of course $`N=N_{fp}N_{shell}`$, where $`N_{fp}`$ is the number of fluid particles. Each point in the shell could be copied to more than one neighbour region. The number of images is the number of points that are actually inserted outside the simulation box. The number of extra images gives the number of point images that are not connected to any point inside the simulation domain, i.e. their presence does not affect the periodic boundary at all. The small number of extra images compared to the number of fluid particles shows that the criterion used to select image points is very effective.

Obviously, it is important that the number of images, which represent a significant computational overhead, is small compared to the number of points in the box. This is easily obtained in two dimensions where for 1000 points we have 190 images (less than 20%). However, in three dimensions, to have the same percentage we need more than 50000 points. This severely limits the performance of the current code with full periodic boundary conditions in three dimensions.

## 5 Maintaining the Voronoi tessellation dynamically

The dynamical maintenance of the tessellation is based on the reconnection algorithm (see ). The reconnection algorithm locally restores the Delaunay property on a valid triangulation. It is very effective, but we need to guarantee that after all the points are moved the tessellation does not present invalid triangles (tetrahedra). It is not simple to assure this condition. Usually, by reducing the time step it is possible to control the validity of the triangulation. However, in our opinion this method is not really stable, and brings with it the additional problem that the reduction of the time step effectively slows down the code. In three dimensions, we prefer to locally restore the tessellation moving one point by one. We are thus able to check that the triangulation is still valid and then to use asynchronous dynamics to update different fluid particles with different time steps. This last option is particulary interesting for a multiscale model like this. In Appendix A, we detail the code for the maintenance of the tessellation in two and three dimensions realised with the CGAL library .

### 5.1 Local reconnection algorithm in two dimensions

Given the algorithmic complexity of the construction of a Voronoi tessellation for $`N`$ points $`(N\mathrm{log}N)`$, it is crucial that we are able to efficiently maintain the Voronoi tessellation over time. We can make significant improvements by utilising the flexibility of the CGAL library. We are able to restore the Delaunay Triangulation (and thus the Voronoi tessellation) locally with the following algorithm (an implementation for the CGAL library is shown in Appendix A), which is an extension of the algorithm used by CGAL for restoring the Delaunay triangulation when individual points are inserted or removed.

We now consider the individual steps of the algorithm. In the first step we ensure that no point can move outside its Voronoi region by restricting the size of the time-step. This is achieved by iterating over the edges of the Delaunay triangles and computing the time-step for each edge such that neither point can move more than a quarter the distance towards the other one by setting $`\delta t`$ to the minimum of all these values. After integrating the stochastic Langevin equations of motion, we then update the position of each point in step $`2`$. It is clear that some triangles may now be invalid Delaunay triangles and so we restore the triangulation to a Delaunay triangulation in step $`3`$. We achieve this as follows: every triangle in the simulation domain has three neighbouring triangles, that is, three neighbouring vertices (by definition, the triangles that belong to the convex hull have the infinite vertex as one of their three neighbouring vertices). A triangle is a valid Delaunay triangle if it satisfies the empty circle property, that is, if none of its neighbouring vertices lies inside its circumcircle. Therefore, we iterate over the finite triangles in the triangulation and check the empty circle property for every triangle; for any triangle $`(p_1,p_2,p_3)`$ which does not satisfy the empty circle property due to a neighbouring vertex, say $`q`$, located within its circumcircle, we simply triangulate the quadrilateral formed by $`p_1,p_2,p_3`$ and $`q`$ the other way round. We note that Algorithm 3 must be implemented such that any new triangles formed are subsequently visited by the iteration process.

Clearly, there is far less work involved in restoring the Delaunay triangulation than with re-computing the whole tessellation for every time-step. Indeed, the above restore Delaunay algorithm is a very efficient way of maintaining a two-dimensional Voronoi tessellation over time.

### 5.2 Reconnection algorithm in three dimensions

In three dimensions, instead of performing a sequence of geometrical flips of the facets in order to locally restore the Delaunay property, we check the validity of the insphere property for all the facets and eventually, if the test fails, we move back the point to the preceding position, remove it and insert in the new position. Algorithm 4 shows this procedure.

For a standard system of fluid particles (see section 6), the number of topological events is usually small, around 10%, therefore the expensive remove-insert procedure is performed only on a very few points. To improve performance the remove-insert task can be used only in the case of invalid triangulation, flipping the facets to restore the Delaunay property for the others. We note that Algorithm 4 spends a large amount of its computational time without modifying anything in the triangulation if it does not not adjust the position of the fluid particle. Because of this, the algorithm is very suited for a data-parallel parallelisation. In fact, a processor acting on a fluid particle does not invalidate the memory of any other processor, while the few topological changes can be stored in a queue and performed sequentially by each processor.

### 5.3 Asynchronous time dynamics

The physical model does not present any restrictions regarding its temporal update. In principle, but not actually implemented in the code, we can add multiscaling in time to multiscaling in space. To each fluid particle is associated a length scale given by its volume, while its time scale can be fixed by its velocity. Therefore, we can set the time scale for each fluid particle and move it on this basis. Sequential temporal updating has already been studied in particle based methods and solved with a hierarchy of time scales of power of two . In practice, given the time scale for fluid particle $`i`$, $`dt_i`$, and the largest time scale $`dt_s`$, then the time steps are chosen such that

$$dt_i=\frac{dt_s}{2_i^n},$$

where $`n_i`$ refers to the time bin of particle $`i`$. Synchronisation is maintained at the end of each large time step $`dt_s`$. Then a new large time scale is computed, and so on. The Voronoi local maintenance algorithm can take advantage of this time advancement given its locality.

## 6 Simulation set-up

In previous dissipative particle models, *a priori* length and time scales were fixed by the cut-off radius . For this model, the only limiting condition is that the number of molecules $`N_k`$ in each DP must be large enough such that $`1/N_k1`$. For the purpose of the present paper we have chosen to work with an average number of molecules per dissipative particle larger than 500 hundred. There are of course no upper bounds to the number of molecules we may assign inside a DP, but when this number is sufficiently large, the fluctuations disappear and the model becomes a standard Lagrangian hydrodynamics code. The input parameters are the size of the simulation box $`L`$, the number of dissipative particles $`N_{DP}`$, the molar mass $`mol`$, the temperature $`T`$ and the density $`\rho `$.

In our three-dimensional simulations, we set up the parameters appropriate for argon. A set of numerical values for water and argon are expressed in cgs units in Table 2. The equation of state contains the contribution of the intermolecular potential in the pairwise approximation. The size of the simulation box is set up depending on the number of DPs and the number of molecules per DP and the tessellation initialised with the Voronoi tessellation composed of regular hexagons or a random configurations of points (three dimensions) in periodic boundary conditions.

We ran a set of equilibrium simulations to benchmark the code in order to establish the computational limits of the model and the Voronoi tessellation on our machine. Because of the thermodynamic consistency of the model, the temperature obtained from the internal energy of the fluid particles is equal to the kinetic temperature, i.e. the temperature of the set of particles, computing the total kinetic energy with zero mean and converting the internal energy to the temperature. We use the set of equations (1). In three dimensions, the equilibrium simulations are performed starting from an isothermal-isobaric configuration with fixed internal energy corresponding to $`300K`$ for argon but with zero momentum and equal pressures for all fluid particles. The number of fluid particles is 1000, while the number of molecules per fluid particle is set by varying the size of the system. In these simulations, we set an average of about 500 molecules per DP and a three dimensional simulation box with side $`2650\times 10^8cm`$. The simulations showed that the DPD temperature of the system increases until it reaches 300K beyond which value the dissipative term causes the extra kinetic energy to be dissipated, operating like a thermostat. This gives strong support for the fluctuation-dissipation relations in three dimensions as well, confirming for the three-dimensional model what has been tested already in two dimensions . The temperature is one of the first quantities to reach equilibrium. For momentum, mass density, pressure, volume, etc. the equilibrium distributions can be analytically computed and the predictions confirmed via a more extensive series of numerical simulations which have been reported elsewhere .

## 7 Parallel implementation and performance in two dimensions

Whereas a parallel code can be written in two and three spatial dimensions following exactly the same strategy, we have so far implemented and tested the parallel performances of the model only in two dimensions. In fact, the parallel code uses the same Algorithm 1 to isolate boundary points as the sequential code. However, whereas in the sequential code the image points are the points in the boundary region of the opposite face, in the case of the parallel implementation they are the boundary points on an adjacent processor. The use of image points gives a straightforward parallel implementation with a message passing paradigm that in this case has been accomplished using the MPI (Message Passing Interface)  standard library. The big advantage of the image points solution is that this parallelisation is independent of the library used. The library is used as a software engine to compute the tessellation only.

The model contains three phases of communication:

* to construct the tessellation using the sequential code for periodic boundary conditions, but communicating the ghost region to the adjacent processor;
* to move the particles to the corresponding processor when the particles fall outside the processor’s domain due to the position update;
* to preserve the force symmetry between two particles in different domains due to the stochastic force term in the Langevin equations.

The performance analysis of the code is reported for a Cray T3E (256 processors Alpha 600MHz) in Figure 8 and Table 3. All the benchmarks are run with the same configuration of $`128000`$ fluid particles (Npart) for $`50`$ iterations.

The simulation is a relaxation towards equilibrium of the isothermal model  starting from a random initial configuration. The timings have been taken measuring the CPU-clock time for the routines: move updates the position and communicates the particle outside the boundary, construct constructs the tessellation, integration integrates the Langevin equations and communicates the stochastic force, comm is the total time of communication without the time for buffering and total measures the total time without Input-Output. All the times measured are inclusive of communication, which is reported separately for comparison; the time for one processor is the time of the sequential version without any overhead due to buffering.

The important point which we wish to stress here is that although the construction and maintenance of the Voronoi tessellation represents a significant computational overhead due to the algorithmic complexity of $`NlogN`$, with this parallel implementation the speed-up obtainable is very promising. This is an intrinsic feature of the model. Each fluid particle interacts only with its nearest neighbour, thereby reducing the amount of communication. Previous standard dissipative particle methods  lacked this strong locality feature, interacting with many more neighbouring particles and severely limiting the parallel performance .

Looking at the speed-up index, it is worth noting that the super-linearity shown was expected because the algorithmic complexity of the construction is not linear, but $`Nlog(N)`$. The construction of the tessellation for $`N/2`$ takes less than half of the time spent for $`N`$ dissipative particles.

The present parallel implementation may suffer from load imbalance depending on the nature of the multiscaling in a given application. If the length scale is set by the polymer molecule and the polymer density is uniform in the simulation region then the parallel implementation should be well balanced. A different situation may arise in the case of multiscaling due to the flow structure, when some processors may have a higher computational load than others. In this case the parallel performance will drop to the slowest processor’s performance. A dynamic repartitioning of the simulation region among processors should reduce this effect, but the repartitioning routine would necessitate further communication, thereby producing an additional overhead.

## 8 Computational limits

In the preceding sections, we described several algorithms that we have used to implement boundary conditions and dynamic reconstruction of the tessellation. We showed as well the good parallel performance of the model due to the locality of its interactions. Now we want to estimate the computational limits of this model in view of future applications to complex fluids flows. The CGAL library version 2.3 is used for these runs.

In two dimensions, a sequential code can easily handle 50,000 fluid particles. Considering the good parallel performance, it is possible to run very large scale simulations. The situation in three dimensions is quite different. The poor efficiency of periodic boundary conditions limits the overall performance of the code. In Table 4, we report the computational time in seconds for a set of runs with increasing number of fluid particles, for 10 iterations. The simulation is a simple equilibration of the fluid particle system. We also report the average number of topological events (reconnections) that are performed for each iteration, the time to move and maintain all the fluid particles, the Eulerian time corresponding to the time needed to solve the Langevin equations once the tessellation and the volumes have been computed, and finally the time to impose boundary conditions.

Table 4 suggests that while the maintenance of the tessellation is reasonable compared to the time needed to integrate the Langevin equations, the boundary conditions do not allow to increase the size of the simulation to more than few thousands fluid particles. For the maintenance task, a non-local reconnection algorithm, which restores the Delaunay property after the entire set of points are moved, would be at least 3 times faster in the move task. However, the stability would be affected and the possibility of asynchronous time update is also lost. Instead, our code proved to be very stable. Further optimisations are of course possible using local flips and remove-insert procedures only in pathological configurations.

## 9 Conclusions

Within the framework of fluctuating hydrodynamics codes for complex fluids, we have described our present algorithmic approach to multiscale dissipative particle dynamics. The algorithms to dynamically maintain the Voronoi tessellation, impose boundary conditions and develop parallel codes for the Voronoi dissipative particle dynamics model have been described in two and three spatial dimensions. The use of publicly available computational libraries  is an advantage compared to direct coding because the evolution of the library is automatically incorporated into the model. This is an important factor at this stage of research in computational geometry where a really stable and efficient library is still under development. In this context, CGAL offers a first interesting approach, delivering a usable and powerful tool for scientists.

The reconnection algorithm used to dynamically maintain the tessellation in two dimensions is stable as long as it is applied on a valid albeit not Delaunay triangulation. This is not easy to guarantee even reducing the time step. In three dimensions, our approach in Section 5.2 has a limited additional cost, but offers the advantage of great stability. The overall efficiency is good for running simulations in real applications.

Boundary conditions are evidently not yet efficiently addressed in three dimensions. Periodic boundaries take most of the time in the simulation. The image points approach for boundaries is appealing because it provides at the same time periodic, fixed wall boundary and a domain decomposition parallel implementation, but it does not perform efficiently due to the large number of image points and the complexity and cost of the point removal from the tessellation in three dimensions. However, it is possible to use periodic boundary in the $`x`$ direction and a wall surface boundary on the $`y`$ and $`z`$ to significantly reduce the amount of image points. This corresponds to flowing a fluid in a duct, where shear flows can be applied to compute the rheological properties of the fluid.

Alternative implementations of three dimensional periodic boundary conditions can be envisaged. Although, we have not actually implemented it, it is nevertheless interesting to consider here an alternative option. Imposing periodic boundary conditions on a bounded region of space can be thought as actually bending the Euclidian space of the region to form a torus. This is actually equivalent to changing the topology of the region which, in fact, might lead to problems of finite size effects (in case the triangulation becomes locally non-euclidian). We have implemented boundary conditions by imposing periodicity on the Euclidian space around a box of fixed size. This makes use of image points outside the box in order to obtain the periodic boundary. A second alternative is to use a three dimensional torus, without image points, but directly connecting the fluid particles from one side to the other. Of course, this necessitates some care when handling the coordinates. A possible implementation of this solution requires changing the insphere function which is used to construct the tessellation. This is related to the distance function, which could be implemented as $`d_{torus}(𝐏_1,𝐏_2):=norm((𝐏_1𝐏_2)mod(𝐋/2)),`$ where $`𝐋=\{L_x,L_y,L_z\}`$ represents the three linear dimensions of the box and the modulus operation $`mod`$ is applied to each component. This function must consider the toroidal topology, in order that points close to the boundary of the previous Euclidian box are still considered close as they are in the torus. This implementation would require substantial modifications of the CGAL library.

G.D.F. thanks Queen Mary, University of London, and Schlumberger Cambridge Research for funding his Ph.D studentship. We thank CGAL for help and support, and P. Español and M. Serrano for useful discussions and exchange of preprints.

## A Maintaining the Voronoi tessellation

In two dimensions, we are able to restore the Delaunay triangulation after every time-step with the algorithm listed below, which is an extension of the restore\_delaunay function in Delaunay\_triangulation\_2.h. We first compute the time-step such that no Delaunay vertex can move outside its Voronoi region. This helps to ensure that we have a valid triangulation, but it does not guarantee it. We integrate the stochastic Langevin equations of motion and then update the position of each vertex by calling the CGAL set\_point function. For each vertex v, we then circulate around the incident faces, checking that v is not located inside the circumcircle of each neighbour of v (external flip) and that no vertices adjacent to v lie inside the circumcircle of each face (internal flip). For any invalid Delaunay triangle, we simply triangulate the quadrilateral the other way by calling the CGAL flip function. Thus we are able to efficiently restore the two dimensional Voronoi tessellation locally.

```

void restore_delaunay(Vertex_handle vh) {
        int i;
        Face_handle f=vh->face(),next,start(f);
        do {
                i=f->index(vh);
                if(!is_infinite(f)) {
                   if(!internal_flip(f,cw(i))) external_flip(f,i);
                   if(f->neighbor(i)==start) start=f;
                }
                f=f->neighbor(cw(i));
           } while(f!=start);
}

void external_flip(Face_handle& f, int i) {
        Face_handle n=f->neighbor(i);
        if(ON_POSITIVE_SIDE !=
           side_of_oriented_circle(n,f->vertex(i)->point())) return;
        flip(f,i);
        external_flip(f,i);
        i=n->index(f->vertex(i));
        external_flip(n,i);
}

bool internal_flip(Face_handle& f, int i) {
        Face_handle n=f->neighbor(i);
        if(ON_POSITIVE_SIDE !=
           side_of_oriented_circle(n,f->vertex(i)->point())) return false;
        flip(f,i);
        return true;
}

```

In three dimensions, the in\_sphere property is checked for all the facets incidents to a given vertex vh. If a topological event is found then the point is removed and inserted again. The facets are inserted in a set in order to avoid double checking the same facet. An implementation is listed below.

```

typedef std::set<Facet,compare_facets> Facets_cont;

bool move_reconnection(Vertex_handle& vh,const Vector& dr) {
    Point old_point=vh->point();
    Point new_point=old_point;
    move(new_point,dr);
    vh->set_point(new_point);
    if (check_topological_event(vh))
        {vh->set_point(old_point);move(vh,dr);return false;}
    return true;
    }

bool check_topological_event(const Vertex_handle& vh) const {
    Facets_cont facets;
    composing_facets(vh,facets);
    Facets_cont::iterator fi=facets.begin(),fiend=facets.end();
    for (;fi!=fiend;fi++) if (facet_in_sphere(*fi)) return true;
    return false;
}

void composing_facets(const Vertex_handle& vh,Facets_cont& facets)
const {
    std::set<Cell_handle> cells;
    this->incident_cells(vh,cells);
    std::set<Cell_handle>::iterator ci=cells.begin(),cend=cells.end();
    for(;ci!=cend;ci++)
    {
    facets.insert(Facet(*ci,0));
    facets.insert(Facet(*ci,1));
    facets.insert(Facet(*ci,2));
    facets.insert(Facet(*ci,3));
    }
}

bool facet_in_sphere(const Facet& f) const {
 return (ON_BOUNDED_SIDE==
        side_of_sphere(f.first,
                      (*(f.first->mirror_vertex(f.second))).point()));
}

```


format: text

text:
# The Redshift Evolution of the 2-8keV X-ray Luminosity Function

## 1. Introduction

Hard X-ray surveys provide the most direct probe of supermassive black hole (SMBH) accretion activity, since high energy X-rays can penetrate extremely large column densities of gas and dust. Early analyses of deep Chandra X-ray surveys found that optical selections of active galactic nuclei (AGN) had missed large numbers of accreting SMBHs. Barger et al. (2001b) showed that the redshift history of SMBH growth is more strongly peaked to low redshifts than had been inferred from optically-selected samples.

With the high spatial resolution and energy sensitivity of deep Chandra observations, X-ray samples can be selected at both low and high redshifts in the same rest-frame hard energy band and their optical counterparts, and hence spectroscopic redshifts ($`65`$% can be identified), unambiguously determined. The hardest band for which the hard X-ray luminosity function (HXLF) can presently be determined over a wide range of redshifts is rest-frame $`28`$ keV. Although some Compton-thick sources may still be omitted, the accounting will be far more complete than that of any other available sample.

In this paper we use the 1 Ms exposure of the Chandra Deep Field North (CDF-N) and three $`100`$ Ks exposures of the A370, SSA13, and SSA22 fields, together with previous ROSAT and ASCA data, to construct the rest-frame $`28`$ keV luminosity function from $`z=0`$ to 4 and model the SMBH accretion history. We assume $`\mathrm{\Omega }_M=1/3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=2/3`$, and $`H_o=65`$ km s-1 Mpc-1. We use $`L_x`$ to denote the rest-frame $`28`$ keV luminosity.

## 2. X-ray Sample Selection

The faintest data are taken from the 1 Ms CDF-N catalogs of Brandt et al. (2001; hereafter B01) and Barger et al. (2002). We exclude known stars and the small number of extended X-ray sources. The flux limits are $`3\times 10^{17}`$ ergs cm-2 s-1 ($`0.52`$ keV) and $`2\times 10^{16}`$ ergs cm-2 s-1 ($`28`$ keV). Below a flux of $`1.6\times 10^{14}`$ ergs cm-2 s-1 ($`28`$ keV) and $`6\times 10^{15}`$ ergs cm-2 s-1 ($`0.52`$ keV) we include only sources lying within $`8^{}`$ of the aimpoint where the redshift identifications are most complete; above these fluxes we include all the sources. This restricted sample contains 199 (52%) sources in the $`28`$ keV band and 219 (57%) in the $`0.52`$ keV band, where the parentheses are the percentages with redshifts. We supplement the sample at fluxes above $`2\times 10^{15}`$ ergs cm-2 s-1 ($`28`$ keV) and $`3\times 10^{16}`$ ergs cm-2 s-1 ($`0.52`$ keV) with data from the SSA13 field (Barger et al. 2001a), the A370 field (Barger et al. 2001b), and the SSA22 field (A. J. Barger et al., in preparation). This adds a further 67 (72%) sources in the $`28`$ keV band and 88 (60%) in the $`0.52`$ keV band. We further supplement the sample with the larger area bright flux samples in the soft band from ROSAT (Lehmann et al. 2001) and the hard band from ASCA (Akiyama et al. 2000). These have highly complete spectroscopic identifications.

For the conversion of counts to flux in the CDF-N, we use the X-ray spectral photon indices $`\mathrm{\Gamma }`$ determined in B01 from the hard-band to soft-band counts ratios. For the conversions in the remaining Chandra fields, we assume a fixed $`\mathrm{\Gamma }=1.2`$, which is approximately the average $`\mathrm{\Gamma }`$ in the B01 sample. For the flux conversions in the ROSAT and ASCA samples, we use the softer $`\mathrm{\Gamma }`$ values assumed by Lehmann et al. (2001) and Akiyama et al. (2000), which are appropriate for higher flux sources.

We determined the solid angle covered by the sample at a given flux by comparing the observed numbers of objects as a function of flux with the measured number counts in the appropriate energy band. For $`28`$ keV we use the averaged counts given in Cowie et al. (2002), and for $`0.52`$ keV we use the power-law fits of Mushotzky et al. (2000), which agree with other determinations (e.g., B01; Rosati et al. 2002). The solid angle covered by the combined $`28`$ keV samples ranges from just under 0.01 deg2 at the faintest fluxes to 5.8 deg2 at the highest fluxes. At $`10^{14}`$ ergs cm-2 s-1 ($`28`$ keV) the solid angle is 0.11 deg2. The maximum solid angle covered by the combined $`0.52`$ keV samples is 0.42 deg2. At $`10^{15}`$ ergs cm-2 s-1 ($`0.52`$ keV) the solid angle is 0.12 deg2.

In Fig. 2a and b we show redshift versus flux for the soft and hard selected samples, respectively, including sources that have not been spectroscopically identified. The solid curves correspond to loci of constant $`L_x`$. Sources with broad-line optical/UV spectra are denoted by a second, larger symbol. Any source more luminous than $`10^{42}`$ ergs s-1 is very likely an AGN on energetic grounds (Zezas, Georgantopoulos, & Ward 1998; Moran et al. 1999), though many of the intermediate luminosity sources do not show obvious AGN signatures in the optical. It is clear from Fig. 2 that the highest X-ray luminosity sources are mostly broad-line AGN.

## 3. Evolution of the AGN X-ray Luminosity Function

In the present section we concentrate on two redshift intervals: $`z=0.11`$, as shown by the dashed lines in Fig. 2b, and $`z=24`$, as shown in Fig. 2a. With these redshift intervals we can compare samples chosen in approximately the same rest-frame hard X-ray energy range of $`312`$ keV (the observed $`28`$ kev band at the center of the low redshift interval) and $`28`$ keV (the observed $`0.52`$ keV band at the center of the high redshift interval). The remaining small differential $`K`$-corrections to rest-frame $`28`$ keV are made by assuming that the spectra can be approximated by $`\mathrm{\Gamma }=1.8`$ power-laws. We have chosen to use this spectral index since the differential K correction is mostly significant in correcting from higher to lower energies where the spectra may be better represented by the unabsorbed intrinsic power law (see Barger et al. 2002). However, the results are only weakly sensitive to this assumption. For a $`\mathrm{\Gamma }=1.2`$ power-law K-correction the HXLF at $`z=24`$ is essentially unchanged while the HXLFs at $`z=0.11`$ shifts to slightly lower luminosities (a factor of 1.27) while the shape remains approximately the same.

We have chosen not to correct for the intrinsic absorption since this too requires an assumption about the intrinsic spectrum and also about the form of the absorption and any scattering. Since we are comparing in the same rest frame at different redshifts this will not affect the HXLF comparison unless there is substantial evolution in the intrinsic absorption. Moreover in the $`28`$ keV rest frame band the average absorption corrections are not very large. If we assume an intrinsic spectral index of $`\mathrm{\Gamma }=1.8`$ and a simple photoelectric absorption the average absorption correction to $`L_x`$ is 1.14 in the $`0.11`$ redshift range and 1.53 in the $`z=24`$ range.

We probe to an $`L_x`$ below $`10^{43}`$ ergs s-1 in the high redshift interval and below $`10^{42}`$ ergs s-1 in the low redshift interval. Both samples are highly complete in redshift identifications at high luminosities. At $`L_x>10^{44}`$ ergs s-1 in the high redshift interval there are 22 spectroscopically identified sources, 21 of which are broad-line AGN. At most there are a further 18 spectroscopically unidentified sources that could lie at these luminosities in this redshift interval; however, a substantial fraction of these are likely at other redshifts. At $`L_x>10^{44}`$ ergs s-1 in the low redshift interval, 85% of the identified sources are broad-line AGN, and there are only three unidentified sources that could lie at these luminosities in this redshift interval. Thus, most of the high luminosity X-ray sources are broad-line AGN rather than obscured AGN.

The differential HXLF as a function of X-ray luminosity and redshift, $`d\mathrm{\Phi }(L_x,z)/d\mathrm{log}L_x`$, is the number of X-ray sources per unit comoving volume per unit logarithmic luminosity. We adopt the $`1/V_a`$ method to estimate the binned HXLF in different redshift shells as

$$\frac{d\mathrm{\Phi }}{d\mathrm{log}L_x}=\frac{\underset{i}{}1/V_a^i}{\mathrm{\Delta }\mathrm{log}L_x},$$
(1)
where the $`V_a^i`$ is the total comoving volume in the redshift bin $`\mathrm{\Delta }z`$ in which source $`i`$ with luminosity $`L_x`$ could lie and still be included in the sample and the summation is over all sources in the given redshift interval and luminosity bin. (Log denotes the base 10 logarithm.)

Our measured HXLFs for the two redshift intervals are shown in Fig. 3 with Poissonian $`1\sigma `$ uncertainties based on the number of galaxies in each luminosity bin. However, since incompleteness is a potentially larger source of error, we recomputed the HXLFs by assigning all the spectroscopically unidentified sources to the center of each redshift interval (i.e., we included all the unidentified sources in both redshift intervals). Dotted ($`z=0.11`$) and dashed ($`z=24`$) lines in Fig. 3 show these maximal limits. Because the spectroscopic identifications are much more complete at the high X-ray fluxes, the associated systematic uncertainties are larger at lower $`L_x`$. In principle photometric redshifts can be estimated for the sources, and at higher redshifts this procedure is effective (A. J. Barger et al., in preparation). However, over the $`z=04`$ range there could be considerable uncertainties in assigning photometric redshifts to sources with substantial AGN contributions, so we prefer our more conservative approach. Analyses of the HXLF using alternate techniques suggest that the lowest luminosity point in the determinations may be biased low (e.g. Page and Carrera 2000, Miyaji et al. 2001) but the effect is small compared to the statistical and systematic uncertainties.

At $`L_x>10^{44}`$ ergs s-1 Fig. 3 shows the familiar rise of the number counts with redshift. These counts are dominated by broad-line AGN. However, at lower $`L_x`$ the $`z=0.11`$ counts lie above the $`z=24`$ counts, or are comparable if we use the maximal correction for incompleteness. Since it is extremely unlikely that all the unidentified sources lie in the $`z=24`$ range (this would result in a large deficiency at $`z=12`$), the present-day number density of sources at $`L_x10^{43}`$ ergs s-1 is higher than it was in the past. This result cannot be significantly changed by highly obscured sources: the dot-dashed line shows the upper bound on the $`z=24`$ HXLF if we selected in observed-frame $`28`$ keV, which corresponds to rest-frame $`832`$ keV, and placed all the unidentified sources in this redshift interval. The resulting upper limit is a factor of $`2`$ higher than the upper limit obtained from the observed-frame $`0.52`$ keV selection.

All the X-ray sources are included in Fig. 3 without regard to optical spectroscopic classification. Previous HXLF analyses have generally focused on broad-line AGN where the counterparts can be determined even in low resolution X-ray data (e.g., Boyle et al. 1998). To make a direct comparison with these earlier studies, in Fig. 3 we plot our spectroscopically identified broad-line AGN HXLFs in the intervals $`z=0.11`$ and $`z=24`$ with those derived by La Franca et al. (2002; $`z=0.5`$) from BeppoSAX data and by Miyaji et al. (2001; $`z=24`$) from ROSAT data. Our HXLFs agree with these brighter HXLFs in the overlap region where there is also some overlap in the datasets used. At $`L_x10^{43}`$ ergs s-1 there is little apparent evolution in the broad-line AGN HXLF, while at higher luminosities the HXLF increases with redshift.

## 4. Cosmic evolution

In Fig. 4a we show the comoving number density evolution with redshift of $`L_x>10^{42}`$ ergs s-1 sources. At these luminosities the probability that the X-ray sources are powered by anything other than AGN is extremely low, so we are mapping the evolution of intermediate luminosity AGN. In the $`z=12`$ redshift bin we show the number densities from both observed-frame $`28`$ keV (open squares) and observed-frame $`0.52`$ keV (solid diamonds); the results are similar, although the soft band may be missing some obscured sources. The solid ($`28`$ keV) and dotted ($`0.52`$ keV) horizontal bars show the number densities obtained if all the spectroscopically unidentified sources are assigned redshifts at the center of each redshift bin. The bars are not consistent with one another because all the unidentified sources are included in all the bins.

The number density of intermediate luminosity AGN does not rise rapidly from $`z0`$ to $`z2.53`$, in contrast to that of higher luminosity AGN (e.g., Miyaji et al. 2000). In fact, the evolution with increasing redshift is flat or declining. Thus, AGN traced by X-ray luminosity show their own version of downsizing, just as star-forming galaxies do (Cowie et al. 1996). High X-ray luminosity sources, which strongly overlap with broad-line AGN (quasar) populations, peak at higher redshifts, while intermediate luminosity sources are still common now.

In Fig 4b we show the evolution with redshift of the $`28`$ keV comoving luminosity density production rate of $`L_x>10^{42}`$ ergs s-1 sources ($`\dot{\lambda }_x`$). Integration of the power-law fits to the HXLFs shown in Fig. 3 give a similar answer: $`1.9\times 10^{39}`$ ergs s-1 Mpc-3 in the $`z=0.11`$ redshift interval and $`2.2\times 10^{39}`$ ergs s-1 Mpc-3 in the $`z=24`$ interval. The largest uncertainty is the redshift distribution of the unidentified sources rather than the small effects of extrapolation outside our observed luminosity range. This incompleteness uncertainty is less than a factor of 2.9 even in the high redshift interval (Fig. 4b).

We parameterize the evolution of $`\dot{\lambda }_x`$ as

$$\dot{\lambda }_x=A\left(\frac{1+z}{1.5}\right)^\alpha ,$$
(2)
where $`\alpha `$ ranges from $`0.4`$ (all unidentified sources in the low redshift interval) to 1.1 (all unidentified sources in the high redshift interval). If only identified sources are included, $`\dot{\lambda }_x`$ is roughly flat until the lowest $`z=0.10.5`$ interval in Fig. 4 where it may drop slightly and $`A`$ $`2\times 10^{39}`$ ergs s-1 Mpc-3. Thus, the $`28`$ keV luminosity density production rate evolves slowly with redshift, as was first noted in Barger et al. (2001b).

The SMBH accretion rate density is now given by

$$\dot{\rho }_{SMBH}=B\frac{\dot{\lambda }_x(1ϵ)}{ϵc^2},$$
(3)
where $`B`$ is the correction from $`L_x`$ to the bolometric luminosity of the AGN, and $`ϵ`$ is the radiative efficiency of the accretion flow. Setting $`B40`$, which seems to be roughly valid for both unobscured and obscured AGN (Elvis et al. 1994; Barger et al. 2001b), and $`ϵ0.1`$, which is close to the maximum possible radiative efficiency,

$$\dot{\rho }_{SMBH}=C\left(\frac{1+z}{1.5}\right)^\alpha ,$$
(4)
where $`C1.3\times 10^5`$ $`\mathrm{M}_{\mathrm{}}`$ yr-1 Mpc-3. Even within the uncertainty in $`\alpha `$, SMBH growth is continuing strongly to the present time. Integrating Eq. 4 for $`\alpha =0`$ gives a present universal SMBH density of $`2\times 10^5`$ $`\mathrm{M}_{\mathrm{}}`$ Mpc-3, comparable to the $`24\times 10^5`$ $`\mathrm{M}_{\mathrm{}}`$ Mpc-3 measured from local galaxies (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Yu & Tremaine 2002). The late formation of the X-ray background reduces the required local SBMH mass density of previous estimates that assumed the bulk of the production was earlier (Fabian and Iwasawa 1999).

In summary, our results show that while the higher accretion mass flow rates that power the most luminous AGN peaked at higher redshift and are now much rarer, lower accretion mass flow rates are still common at the present time. These less luminous events dominate SMBH formation. The longer duration of the lower mass flow rate systems provide constraints on theoretical modelling of the SMBH formation (e.g., Haehnelt, Natarajan, & Rees 1998; Haiman & Menou 2000; Yu & Tremaine 2002) and on the evolution of the accretion disks feeding the sources, which may evolve more rapidly if they are initially of higher mass (Duschl 2002, private communication).

Support came from CXC grants DF1-2001X and GO2-3187B (LLC), NSF grants AST-0084816 (LLC) and AST-0084847 (AJB), the University of Wisconsin Research Committee with funds granted by the Wisconsin Alumni Research Foundation (AJB), the Alfred P. Sloan foundation (AJB), NSF CAREER award AST-9983783 (WNB), and NASA grant NAS 8-01128 (GPG, PI).