text:
# A new general purpose event horizon finder for 3D numerical spacetimes

## 1 Introduction

An event horizon (EH) is defined as a 2+1 surface in 3+1 space, inside of which no null geodesics can reach future null infinity, while outside at least some can. For that reason an EH is a global concept and can in principle only be found when the full history of the spacetime is known. The horizon itself is generated by outgoing null geodesics, that once they have joined onto the horizon will forever stay on it.

In recent years 3D numerical evolutions of binary black hole spacetimes have become stable enough that it is possible in some cases to follow the final merged black hole for a significant time . Therefore it makes sense to start looking for EHs in these numerical spacetimes. However, previously published EH finders for numerical spacetimes  have either been too slow, limited to single black hole spacetimes or have taken advantage of special symmetries in order to handle changes in topology and can not be used in the general 3D case without symmetries.

This paper is organized in the following way. In section 2 I will present the basic ideas and methods of EH finding. In section 3 I will discuss different ways of describing the surface and present the level set description used here. The numerical implementation is described in section 4, where also a cure for a very serious problem is presented. In section 5 I will present tests of the EH finder using analytical data with known location of the EH to test the accuracy of the code while in section 6 I will present tests using numerical data in order to test the robustness of the code. Finally in section 7 I will discuss prospects of future work and uses of the code.

## 2 Event horizon finding

In principle the EH can be found by integrating null geodesics forward in time. Outgoing null geodesics just outside or inside the EH will all diverge away from the EH. Those outside will escape to infinity, while those inside will end up at the singularity. However, as pointed out in  it is practically impossible to follow the generators of the horizon forward in time, since small numerical errors will cause the null geodesics to deviate exponentially away from the EH. This can, on the other hand, be taken advantage of, if the generators are integrated backwards in time, since then the EH will be an attractor. A method based on this was presented in .

For the backward null geodesic integration method, the presence of a small tangential velocity component can cause the null geodesics to deviate from the EH. Even though the null geodesic may return to the EH, it will be in a different position and it might even cross other null geodesics, thereby creating spurious caustics. The cure is to evolve the complete horizon surface as a whole. Since the only way a surface can move is in its normal direction, tangential drift is not an issue. This can be done by representing the surface by a function

$$f(t,x^i)=0,$$
(1)
and requiring this surface to satisfy the null condition

$$g^{\alpha \beta }_\alpha f_\beta f=0.$$
(2)
Expanding equation (2) out yields a quadratic equation for $`_tf`$, which can be solved, giving the following evolution equation for $`f`$

$$_tf=\frac{g^{ti}_if+\sqrt{(g^{ti}_if)^2g^{tt}g^{ij}_if_jf}}{g^{tt}}.$$
(3)
Here the root was chosen so as to describe outgoing null geodesics. Notice that, in contrast to the geodesic equation, this equation does not contain any derivatives of the metric.

Note that the exact location of the EH is unknown at the end of a numerical evolution. However the location of the apparent horizon (AH) can serve as a good initial guess for the EH (see  for details about finding AHs). It will be completely inside the EH, but will be very close if the numerical spacetime is almost stationary. A practical way of locating the EH is to start with the AH as an initial guess for the EH and to use the attracting property of the EH to approach it asymptotically. It is therefore important to note that the surface $`f=0`$ will never exactly coincide with the EH, but will approach it exponentially. In order to estimate when the surface $`f=0`$ is a good approximation, a different initial surface (chosen so as to be completely outside of the EH) can be evolved backwards in time. The EH will then always be located between these two surfaces and when the two surfaces agree to within a small fraction of a grid cell the location of the EH is known with sufficient accuracy.

## 3 Description of the surface

In  either of the following parametrization was adopted to describe the surface in the presented axisymmetric cases

$`f(t,r,\theta )`$ $`=`$ $`rs(t,\theta ),`$ (4)
$`f(t,z,\rho )`$ $`=`$ $`\rho s(t,z),`$ (5)
where equation (5) was used in the colliding black hole case, in order to allow the surface to cross itself thereby describing the locus of generators before they join onto the EH. This was possible due to the high level of symmetry, but would be difficult to generalize to the non symmetric case, since it requires prior knowledge of the location of the caustic points (in this case the symmetry axis).

The choice made in this work, is to avoid any parametrization of the surface that might run into trouble, but to keep the description in equation (1). That is, the horizon surface is described as the 0-level isosurface of the scalar function $`f`$, where $`f`$ is negative inside and positive outside the surface. There are several advantages to this choice. The main advantage is that changes of topology (such as when black holes merge) are handled naturally with no special symmetry requirements. What happens in these cases is simply that the number of regions with negative $`f`$-values changes. Another advantage is that the function $`f`$ is defined at the same grid points as the numerical metric so that no interpolation is necessary in order to evolve the surface. The last point can also be seen as a disadvantage, since when $`f`$ is defined on fixed grid points only, there are in general no grid points exactly on the $`f=0`$ isosurface. This means that in order to analyze the surface it is necessary to find it first. A second (and more serious) problem is the fact that during evolution of $`f`$ its gradients steepen. This can be seen in the following way: The level set function $`f`$ actually defines an infinity of surfaces with different iso-values. Thus for example the triplet of iso-values $`f=[0.5,0.0,0.5]`$ defines three outgoing null surfaces that will all move towards the EH, when evolved backwards in time. This means that these three surfaces will move towards each other and since they are defined by constant values of $`f`$, the gradients will steepen. This will cause numerical problems if it is not addressed. Another potential problem is illustrated in A, where it is shown that for at least one choice of the metric (Schwarzschild in Kerr-Schild coordinates), the evolution of $`f`$ is ill conditioned.

## 4 Numerics

The evolution of equation (3) is performed using the Method of Lines (MoL) with either second order Runge-Kutta (RK2) or 3-step Iterative Crank Nicholson (ICN). The spatial derivatives of $`f`$ are calculated using one sided second order derivatives. The directions of the derivatives are determined by different schemes depending on the problem. If the shift is zero, the directions are determined using the values of $`f`$ itself according to the scheme in B. On the other hand if the shift is non-zero, the direction of the shift has to be taken into account in order to preserve stability, the main reason being that the evolution equation contains terms of the form $`g^{ti}_if`$. In this case the stencil direction is taken to be opposite the shift direction, since I am evolving backwards in time.

In order to be able to handle numerical data with excised regions, the code supports dynamic interior excision regions, that can be used even when the numerical data is defined everywhere. Since the EH normally does not occupy the whole numerical domain it is only necessary to output a rectangular box that contains the EH at all times. This can result in a very significant reduction of the amount of required disk space for the numerical metric data. Furthermore to save computation time, once the numerical data has been read in, the computationally active region can be further reduced so that only a smaller rectangular box surrounding the black holes is used. As the holes move across the grid, this smaller active region moves along with them by de-activating and activating grid points as necessary. At the interior boundaries and the boundaries of the rectangular box, one sided second order derivatives are used as well.

The initial guess for $`f`$ at $`t=T`$, where $`T`$ is the time at the end of the numerical evolution, is normally chosen as a sphere of radius $`r_0`$ centered at the point $`P=(x_0,y_0,z_0)`$

$$f(T,x,y,z)=\sqrt{(xx_0)^2+(yy_0)^2+(zz_0)^2}r_0,$$
(6)
which has the nice property that

$$|f|=1,$$
(7)
except at the origin where it is non differentiable. The parameters for the initial guess are chosen to be as good a guess for the EH as possible, while making sure that $`f=0`$ is either completely inside or outside of the EH. If the EH is judged (for example by looking at the shape of the AH) to be significantly nonspherical it is also possible to use an arbitrarily oriented and shaped ellipsoid as the initial guess. This has not been necessary so far since, with the currently used gauge conditions, the final horizon is almost always close to being a coordinate sphere.

As mentioned, there is a problem with gradients of $`f`$ steepening during the evolution. For that reason $`f`$ is re-initialized regularly during the evolution so that it again satisfies equation (7) approximately. To obtain this re-initialization the evolution equation

$$\frac{f}{\lambda }=\frac{f}{\sqrt{f^2+1}}\left(|f|1\right)$$
(8)
is evolved in an unphysical parameter $`\lambda `$ until a steady state has been achieved. The factor in front of the parenthesis in equation (8) consists of two parts. The numerator is present to make sure that no evolution takes place when $`f=0`$, since I do not want to move the surface during re-initialization. The denominator is there to make sure that the Courant condition for stability is limited by the constant value, 1, instead of the maximum value of $`f`$. The term in the parenthesis makes sure that the evolution stops when $`|f|=1`$. The evolution of equation (8) is also done with MoL, but since it is only the final steady state that are significant, it is sufficient to use a simple Euler scheme. However, a second order Runge-Kutta scheme has been implemented as well. The spatial derivatives used for the re-initialization are second order upwinded according to the scheme in B. Since the re-initialization equation does not depend on the shift, the shift direction is ignored.

If the re-initialization is performed often enough, I do not have to worry about potential stability problems due to the use of one sided differences of $`f`$ at the boundaries. This is the main reason, that it is possible to evolve only a small region around the $`f=0`$ surface. Note also, that without re-initialization all different $`f=const`$ surfaces would be evolved as null surfaces, but every time re-initialization is done all surfaces except $`f=0`$ will be changed, so they can not any more be considered to be null surfaces. In this way $`f=0`$ is picked out as the surface of interest.

The re-initialization equation (8) is similar to the one used in , the differences being introduced to make the surface move as little as possible, while still being reasonably fast.

This re-initialization scheme works well most of the time, but it has some problems when there is only a few points in either directions inside the surface. In those circumstances the level 0 isosurface can move significantly outwards. This typically occurs just before the topology of the surface changes. To avoid this, a scheme has been implemented to detect when this happens and avoid doing the re-initialization until after the topology has changed. The detection scheme essentially consists of searching for all places where 1) $`f`$ is negative (i.e. inside the surfaces) and 2) there is a local extremum of $`f`$ in all directions. Since the value of $`f`$ in such a grid point is approximately equal to the signed distance from that grid point to the level 0 isosurface, the absolute value of the maximum of $`f`$ in these points can be used to estimate the minimum width of the region with $`f<0`$. If this width is less than a given threshold, the re-initialization is not done.

In order to calculate the areas of the EHs, it is necessary to locate points on the surfaces. This task is complicated by the fact that in principle, the number of individual surfaces at any given time is not known. Therefore points on the surfaces present are located in several separate steps.

First the number of surfaces in the data is found. This is done using a “flooding” algorithm using an integer mask (surface\_mask) initialized to zero at the grid points of the level set function, $`f`$, and an integer counter (surface\_counter) initialized to zero. The algorithm proceeds as follows. Locate a grid point where $`f<0`$ and surface\_mask=0. If such a grid point does not exist there are no surfaces, so exit. Otherwise increment surface\_counter and mark this point surface\_mask=surface\_counter. Repeatedly find connected grid points inside the same surface, i.e. locate points with $`f<0`$ and at least one neighbour point with surface\_mask=surface\_counter and mark those likewise until no more such points can be found. In this way all points inside the first surface has been marked. Now check if another point with $`f<0`$ and surface\_mask=0 exists. This point must then be located inside another surface. Repeat the above marking procedure until all points with $`f<0`$ has been marked with a surface number.

Next an approximate coordinate centroid is found for each surface using the marked points. This is used as the center for a polar coordinate system used to explicitly parametrize the corresponding surface as $`r(\theta ,\varphi )`$. For each chosen direction (given by $`\theta `$ and $`\varphi `$) the radius of the surface is found using interpolated values of $`f`$ in a Newton iteration root finder111A Hermite interpolation polynomial is used, since not only the function but also the derivative has to be continuous in the Newton iteration..

With this explicit parametrization of the surface, the area, centroid and circumferences can be easily computed. Unfortunately it is not always possible to find this explicit parametrization. As shown later in Section 6.2 the EH can in some cases be so distorted that with any choice of center, $`f`$ would be multivalued along some angular directions. For this reason a more general surface finding and integration routine using a direct triangulation of the surfaces is under consideration.

I would like to emphasize that the use of interpolation is only necessary when finding points on the surfaces for analysis purposes and is not used at all in the evolution and re-initialization of the level set function.

The EH finding algorithm has been implemented as a thorn in Cactus  using the Einstein Toolkit  and is fully parallelized. Though it might happen that the computationally active region is contained completely on one processor, this is usually not a serious concern since the thorn uses a lot less memory than a comparably sized numerical evolution and can therefore be run on a lot less processors. Also it is not necessary to run the EH finder on the full computational grid of the spacetime evolution. It suffices to output just the part of the grid that completely contains the EH (or horizons) within it at all times, also cutting down on the required disk space for output.

The source code for the EH finder, will be released for public use in the near future.

## 5 Tests using analytical data

In this section I will present tests of the code using analytic spacetimes, where the location of the EH is known analytically at all times.

### 5.1 Non-rotating black hole in isotropic coordinates

The first test uses the metric for a non-rotating black hole in isotropic coordinates. In these coordinates the metric is

$$s^2=\left(\frac{1M/(2r)}{1+M/(2r)}\right)^2t^2+\left(1+\frac{M}{2r}\right)^4[r^2+r^2(\theta ^2+\mathrm{sin}^2\theta \varphi ^2)],$$

where $`M`$ is the mass of the black hole. In these coordinates the horizon is a stationary sphere with radius $`r_{EH}=M/2`$. The equation of motion for radial outgoing null geodesics is

$$\frac{r}{t}=\frac{1M/(2r)}{\left(1+M/(2r)\right)^3}.$$

Looking at null geodesics close to the EH $`r=M/2+ϵ`$ where $`ϵ<<M/2`$ and expanding it can be seen that

$$\frac{ϵ}{t}\frac{ϵ}{4M}.$$

So integrating forward in time the null geodesic diverges exponentially from the EH. However integrating backwards in time it should be expected that the null surface converges exponentially to the location of the EH with an e-folding time of $`4M`$.

I performed runs at three different resolutions $`\mathrm{\Delta }=0.05M,0.025M,0.0125M`$ for $`M=1`$ with an initial radius for the null surface equal to $`r_0=0.4M`$ (i.e. 2, 4 and 8 gridspacings away from the true EH). In figure 1 I show the maximum distance of the surface from the true location of the EH divided by the gridspacing as a function of time in a logarithmic plot.

Initially there is a clear exponential convergence until the finite difference solution is reached with a numerical e-folding time of $`4.13M`$, agreeing nicely with the analytic expectation. At higher resolution the finite difference solution is closer to the analytic solution. At a resolution of $`\mathrm{\Delta }=0.05`$ the maximum error at $`T=100M`$ is $`0.0097\mathrm{\Delta }`$. At the higher resolution of $`\mathrm{\Delta }=0.025`$ the corresponding error is $`0.0056\mathrm{\Delta }`$,and at the highest resolution of $`\mathrm{\Delta }=0.0125M`$ the error is $`0.0031\mathrm{\Delta }`$. Looking at the ratio of the error at successive resolutions, it can be seen that the code is approaching second order convergence. Looking at the absolute errors it is worth noting that the method is able to locate the EH to less than 1/100 of a grid spacing when the radius of the EH is only 10 grid spacings.

### 5.2 Non-rotating black hole in Kerr-Schild coordinates

The second test uses the metric for a non-rotating black hole in Kerr-Schild coordinates. In these coordinates the metric is

$$s^2=\left(1+\frac{2M}{r}\right)t^2+\frac{4M}{r}tr+\left(1+\frac{2M}{r}\right)r^2+r^2(\theta ^2+\mathrm{sin}^2\theta \varphi ^2).$$
(9)
Just like in the previous section it is easy to find the equation of motion for radial outgoing null geodesics and analyse the limit of null geodesics, close to the EH (in these coordinates $`r_{EH}=2M`$). The result is that the e-folding time is $`4M`$ for these coordinates as well.

Here I performed runs at three different resolutions $`\mathrm{\Delta }=0.2M,0.1M,0.05M`$ for $`M=1`$ with an initial radius for the null surface equal to $`r_0=1.8M`$ (i.e. 2, 4 and 8 gridspacings away from the true EH). In figure 2 I show the maximum distance of the surface from the true location of the EH divided by the gridspacing as a function of time in a logarithmic plot.

Again there is a clear exponential convergence until the finite difference solution is reached with a numerical e-folding time of $`4.01`$ in excellent agreement with the analytic expectation. Also in this case, the higher the resolution the closer the finite difference solution is to the analytic solution. The maximal error divided by the grid spacing at $`T=100M`$ is 0.01, 0.0057, 0.0031 respectively at the 3 different resolutions used. I.e. very similar to the results for Schwarzschild in isotropic coordinates.

### 5.3 Rotating black hole in Kerr-Schild coordinates

The third test uses the metric for a rotating black hole in Kerr-Schild coordinates. In these coordinates the EH is an ellipsoid with the following equation

$$d\frac{x^2+y^2}{r_+^2+a^2}+\frac{z^2}{r_+^2}1=0,$$
(10)
where $`r_+=M+\sqrt{M^2a^2}`$, $`M`$ is the mass and $`a=J/M`$ is the angular momentum per unit mass. For points not located on the EH the deviation, $`d`$, defined in equation (10) can be used as a measure of the error. However it will not be a direct measure of the distance from the surface to the EH. Usually it will be somewhat larger than the distance.

I performed runs at four different resolutions $`\mathrm{\Delta }=0.2M,0.1M,0.05M,0.025M`$ for $`M=1`$ and $`a=0.8`$ with an initial radius for the null surface equal to $`r_0=2.0`$ in all cases. In figure 3 I plot the deviation as defined in equation (10) divided by the gridspacing as a function of time in a logarithmic plot.

At all resolutions there is a nice exponential convergence until the finite difference solution is reached with a numerical e-folding time of $`5.18`$. The deviation, $`d`$, divided by the grid spacing at $`T=100M`$ goes down at higher resolution. For the four resolutions it is $`0.029,0.013,0.0062,0.0031`$ respectively, showing perfect second order convergence.

## 6 Tests using numerical data

In this section I will present results from tests based on metric data from highly dynamic spacetimes evolved numerically with Cactus. The exact location of the EH is not known in these numerical spacetimes, but I compare to other published results to the extent possible. At the same time it will be a test of the robustness of the code and how well it can handle changes of topology in multiple black hole spacetimes.

### 6.1 Misner data with $`\mu =2.2`$

In  initial data describing two black holes, initially at rest, was presented. The initial separation between the two black holes is described by a parameter $`\mu `$. These types of initial data have been extensively studied in the literature and especially the EH for these data was studied using numerical data from axisymmetric 2D codes and from 3D codes in ,  and .

The spacetime was evolved using the BSSN formulation of the Einstein equations developed in  and  and implemented in Cactus in . In order to be able to compare with the results in  I used maximal slicing. The run was done at fairly low resolution of $`\mathrm{\Delta }=0.128M`$ in order to be performed with a modest amount of output on a workstation. With the recent development of improved shift conditions in , it turned out to be advantageous to use the hyperbolic version of the gamma freezing shift. In this way the slice stretching was controlled, the expansion of the horizon in coordinate space was strongly reduced and it was possible to evolve for a very long time. In the present case, the evolution was halted at $`T=70M`$, since this is long enough to track the EH accurately for the first $`50M`$ of the evolution.

At $`T=70M`$ the AH was found in the numerical data and it turned out to be nearly spherical with a coordinate radius of $`r_{AH}=2.84M`$ (the ratio of the equatorial to the polar circumference was $`C_{eq}/C_{pol}=1.00074`$). Two different runs with the EH finder were therefore started at $`T=70M`$; one run where the initial surface described by $`f=0`$ was a sphere of radius $`r_I=2.5M`$ completely contained within the AH and another run with an initial sphere of radius $`r_O=r_I+0.5M=3.0M`$ which seemed to be a safe guess for a surface completely outside of the EH. This was confirmed by the fact that initially the inner surface was expanding while the outer surface was contracting when the surfaces were evolved backwards in time.

Both of these surfaces were evolved successfully backwards in time, through the change of topology all the way to the initial data slice. The change of topology occurred at around $`T2.75M`$.

In  it can be seen from figure 11 that the topology changes between $`T=2M`$ and $`T=3.3M`$. However, it turns out that a more precise transition time is $`T3.2M`$ , which is somewhat later than the transition time $`T2.75M`$ found above. Since in both cases maximal slicing was used, the difference must come from a difference in the initial slicing. In  the Čadež lapse profile, which is zero at the throats, was used initially, while in the current work the initial lapse was one. Therefore in this case the evolution proceeds somewhat faster until the lapse collapses, compared to the case with the lapse initially collapsed at the throats, resulting in a shorter time for the change of topology. It is not currently possible to use the Čadež initial lapse profile in Cactus, so a direct comparison is unfortunately not possible. However, the results are consistent.

In the left plot in figure 4 I show the total irreducible mass of the horizons for the Misner $`\mu =2.2`$ spacetime as a function of time for the two different choices for the initial surfaces mentioned previously. The total irreducible mass is defined as

$$M_{ir}=\sqrt{\frac{_{i=1}^nA_i}{16\pi }},$$

where $`n`$ is the number of horizons in the spacetime and $`A_i`$ is the area of the $`i`$’th horizon. As can be seen, even though $`M_{ir}`$ of the two surfaces are very different at $`T=70M`$, they are almost indistinguishable on this plot for $`T<50M`$. The fact that the curve crosses the dotted line and the subsequent rise in $`M_{ir}`$ from about $`T=20M`$ to $`T=50M`$ are due to the lack of accuracy of the underlying numerical data. The run was done at comparatively low resolution and there are also reflections from the outer boundary. The initial dip in $`M_{ir}`$ must be due to the low resolution (initially the horizons are approximately $`0.25M`$ in radius) but as the horizons expand in coordinates they become better resolved and the areas, and therefore $`M_{ir}`$, correspondingly more accurate.

The right plot in figure 4 shows the expected exponential convergence in the difference in $`M_{ir}`$ for the same two choices of the initial surfaces. Looking closer it can be seen that the convergence is in fact slightly better than exponential, but since the horizon mass is growing slightly with time, this is to be expected. The noise in the beginning, at $`T<10M`$, is a lot smaller than the truncation error in either the interpolation and the finite differencing used in obtaining the area.

### 6.2 “Head on” collision of 3 black holes

Since the Misner data is axisymmetric, the robustness of the code towards non-symmetric data was not really tested in the previous section. Therefore I performed a numerical run with a three black hole spacetime. To keep things simple I used Brill-Lindquist initial data  that contains no linear and angular momentum. However, I placed the three black hole punctures on a plane in coordinate space that was tilted with respect to the coordinate axis and with slightly different coordinate distances between the black holes. The masses were chosen to be equal $`M_1=M_2=M_3=0.5`$ giving a system with an ADM mass of $`M_{ADM}=1.5`$222For Brill-Lindquist data the ADM mass is just the sum of the individual black hole mass parameters, while the bare masses are different from the mass parameters.. The three position vectors were $`𝐫_1=(1.2,1.4,0.8)`$, $`𝐫_2=(1.6,0.6,0.8)`$ and $`𝐫_3=(0.4,0.8,1.6)`$ giving a system with the center of mass at the origin. The coordinate distances were thus $`d_{12}3.44`$, $`d_{13}3.35`$ and $`d_{23}3.13`$.

The spacetime was evolved similarly to the Misner spacetime in the previous section, though a hyperbolic K-freezing slicing was used instead of maximal slicing, simply because it is much faster. I used a fairly low resolution of $`\mathrm{\Delta }=0.2`$, since the purpose was not to test the accuracy but rather the robustness of the code. The run could then be performed on a workstation with a moderate amount of disk requirements. Using the hyperbolic gamma freezing shift it was possible to evolve until $`T=40=26.7M_{ADM}`$ at which time I stopped the evolution. The run could have continued a bit longer but the Hamiltonian constraint was starting to grow significantly333The growth of the Hamiltonian constraint was contained mostly within the apparent horizon and was rather small in the exterior. and I estimated that the evolution time was long enough to locate the EH accurately enough for the initial $`10M_{ADM}`$ of evolution.

At about $`T=10=6.7M_{ADM}`$ a common AH appeared and at $`T=40=26.7M_{ADM}`$ it had evolved into an approximately spherical shape with a coordinate radius of $`r_{AH}=3.6`$. The deviations from sphericity was less than 1%. Therefore it was safe to assume that a sphere of radius $`r_I=3.5`$ was completely contained within the EH at $`T=40=26.7M_{ADM}`$ and this was chosen as the inner initial guess for the EH, while the outer initial guess was chosen to be a sphere of radius $`r_O=3.7`$.

In figure 5 four frames showing the interesting part of the evolution of the EH can be seen.

Frame (a) shows 3 separated horizons at $`T=2.067M_{ADM}`$. The individual horizons are clearly distorted and it can be seen that they “feel” each other. Frame (b) shows that at $`T=2.1M_{ADM}`$ two of the horizons has merged while the third is still clearly well separated. This has changed at $`T=2.267M_{ADM}`$ shown in frame (c) where the final horizon has merged with the two other. In frame (d) the horizon is shown at $`T=2.667M_{ADM}`$, where it has started to increase in thickness. At later times the EH becomes more and more spherical.

The spatial and time resolution in the metric data is not high enough to find instances of toroidal EHs as predicted in , since these are generally short lived and have sharp features. However, looking at frames (b) and (c), it does not take much imagination to envision that a toroidal EH could appear for a short time, given higher spatial and time resolution.

The plots in figure 6 shows the total irreducible mass of the horizons, $`M_{ir}`$ in the spacetime as a function of time for the two different choices of initial surfaces at $`T=26.7M_{ADM}`$. In these plots I had to exclude a few points corresponding to times where the shape of the horizon is so distorted that describing it with $`r(\theta ,\varphi )`$ around a suitably chosen center is impossible (see e.g. frame (b) in figure 5). Again $`M_{ir}`$ goes beyond the limiting value set by the ADM mass of the system, but in this case the resolution is even lower than in the Misner case so this is not surprising. The convergence of the two surfaces are also here nicely exponential, though the surfaces do not get as close together as in the Misner case due to the shorter evolution time.

## 7 Discussion

The EH finder presented here is a very general and accurate method of locating and tracking EHs. It works both with analytic and numerical metric data and on spacetimes with single or multiple black holes and can handle changes of topology (and even though it has not been demonstrated here, it should be able to handle toroidal EHs). For analytic and stationary spacetimes the horizon can typically be located to within a small fraction of a cell size (on the order of $`1/10\mathrm{\Delta }`$$`1/100\mathrm{\Delta }`$) at typical resolutions and the code shows second order convergence with the grid resolution.

The currently implemented analysis tools are restricted to calculating the area, centroid and polar and equatorial circumferences444The circumferences only makes sense in spacetimes with sufficient symmetries. for all the surfaces present in the data as long as the surfaces can be represented in spherical coordinates as $`r(\theta ,\varphi )`$ around a point inside the surface. As shown in figure 5 this is not always possible, so a more general area integration routine is under consideration.

The robustness of the code was demonstrated by finding and tracking the EH of a three black hole spacetime with no symmetries through two changes of topology.

Analysis tools for calculating the Gaussian curvature of the horizon and for tracking a congruence of individual generators of the horizon surface have been shown in  to be useful for extracting physical information from highly dynamic spacetimes. With a congruence of horizon generators, it is possible to construct the membrane paradigm quantities  and thus to study the shear and expansion of the horizon in detail. It should be possible to track individual generators along with the horizon surface using the information contained in the level set function. I am currently investigating a way to implement this powerful analysis tool.

The prediction of the generic formation of toroidal EHs (at least for a short while) in non-symmetric spacetimes  is interesting. I believe that high resolution is required to see this numerically, but it might be within the realm of possibility with this code.

This EH finder will be applied to a large number of different black hole spacetimes that can be evolved long enough to reach an almost stationary state.

During the final phase of writing this paper, two other papers about EH finding methods  appeared on GR-QC. Tests of the code have been performed on the Hitachi SR8000-F1 at the Leibniz-Rechenzentrum (LRZ), the IBM SP RS/6000 at the National Energy Research Scientific Computing Center (NERSC), the IA-32 Linux cluser at the National Center for Supercomputing Applications (NCSA), various workstations at the Albert-Einstein-Institut (AEI), my laptop and on the new Linux cluster (Peyote) at the AEI. I would like to thank Ed Seidel, Sascha Husa, Denis Pollney, Fransisco Guzman, Ian Hawke, Jonathan Thornburg and Jeffrey Winicour for helpful discussions. On the practical side I would like to thank Ian Hawke for implementing the MoL thorn in Cactus making it easy to write an evolution code, Jonathan Thornburg for implementing the new interpolator in Cactus that can return both the interpolated function and its derivatives at the same time and Thomas Radke for helping me with the IO aspects of the code. I would also like to thank the Tapir group at Caltech for useful comments after my talk during my visit in February. This work was supported by the EU Programme ‘Improving the Human Research Potential and the Socio-Economic Knowledge Base’ (Research Training Network Contract HPRN-CT-2000-00137).

## Appendix A Level set evolution in Kerr-Schild coordinates

In this appendix I will show that, for a non-rotating black hole in Kerr-Schild coordinates, an analytic solution to the evolution equation for the level set function $`f`$ only exist for a limited time when evolving backwards in time. I will only consider the spherically symmetric case $`f=f(r)`$. With the metric given by equation (9) the relevant components of the contravariant metric is

$$g^{tt}=\left(1+\frac{2M}{r}\right),g^{tr}=\frac{2M}{r},g^{rr}=1\frac{2M}{r},$$

and the evolution equation for $`f`$ becomes

$$\frac{f}{t}=\frac{2Mr}{2M+r}\frac{f}{r}.$$
(11)
For $`r=2M`$ this equation has the trivial solution

$$f(t,r=2M)=c,$$

where $`c`$ is a constant depending only on the initial value of $`f`$ at $`r=2M`$. For $`r2M`$ this equation has the general solution

$$f(t,r)=h(u(t,r))=h(tr4M\mathrm{ln}|2Mr|),$$
(12)
where $`u(t,r)=tr4M\mathrm{ln}|2Mr|`$. I now limit myself to $`r<2M`$. At $`t=0`$ the following relation between $`u(0,r)=u_0(r)`$ and $`r`$ can be found

$$u(0,r)=u_0(r)=r4M\mathrm{ln}(2Mr).$$

The function $`u_0(r)`$ has the properties that $`lim_r\mathrm{}=\mathrm{}`$ and $`lim_{r2M^{}}=\mathrm{}`$ and it has exactly one minimum in between at $`r=2M`$ where the value is $`u_0(2M)=2M4M\mathrm{ln}(4M)`$. The range of $`u(0,r)=u_0(r)`$ is therefore the unbounded interval $`[2M4M\mathrm{ln}(4M),\mathrm{})`$, which again is the domain of the $`h(u)`$. The function, $`h(u)`$, does not change in time, so in order for the solution, $`f(t,r)`$ to exist, $`u(t,r)`$ must satisfy $`u(t,r)2M4M\mathrm{ln}(4M)`$. From this it can be seen that a solution exists only at $`r=0`$ for

$$t2M4M\mathrm{ln}20.77259M.$$

There is no limit to the existence of a solution when evolving forward in time. In the solution, the problem shows up as the development of infinite derivatives, $`f(t,r)/t`$ and $`f(t,r)/r`$, at $`r=0`$ in a finite time. Excising a region around $`r=0`$ does not help, as the point of infinite gradient moves outward and will reach the boundary of the excision region in a finite time. The time for the development of infinite gradients is independent of the initial shape of $`f(0,r)`$ and are purely a property of the metric. A non rotating black hole in isotropic coordinates do not have the same property. In that case it takes infinite time to develop infinite gradients.

## Appendix B Second order upwinded derivatives

Define $`a_l`$ and $`a_r`$ as the second order one sided derivatives of $`f`$ in the $`x`$-direction as

$`a_l`$ $`=`$ $`{\displaystyle \frac{3f_{i,j,k}4f_{i1,j,k}+f_{i2,j,k}}{2\mathrm{\Delta }}},`$
$`a_r`$ $`=`$ $`{\displaystyle \frac{3f_{i,j,k}+4f_{i+1,j,k}f_{i+2,j,k}}{2\mathrm{\Delta }}}.`$
Then define the negative and positive parts of $`a_l`$ and $`a_r`$ as

$`a_l^{}`$ $`={\displaystyle \frac{a_la_l}{2}},\text{ }a_l^+`$ $`={\displaystyle \frac{a_l+a_l}{2}},`$
$`a_r^{}`$ $`={\displaystyle \frac{a_ra_r}{2}},\text{ }a_r^+`$ $`={\displaystyle \frac{a_r+a_r}{2}}.`$
If $`f_i0`$ then

$$\frac{f}{x}\{\begin{array}{ccc}a_l^+\hfill & \text{if}& \hfill a_l^+>a_r^{}\\ a_r^{}\hfill & \text{otherwise}.& \end{array}$$

If $`f_i<0`$ then

$$\frac{f}{x}\{\begin{array}{ccc}a_l^{}\hfill & \text{if}& \hfill a_l^{}>a_r^+\\ a_r^+\hfill & \text{otherwise}.& \end{array}$$

Derivatives with respect to $`y`$ and $`z`$ are approximated by corresponding equations. This method of upwinding derivatives is a generalization of the first order method presented in  to second order derivatives. The first order method is shown to be a consistent monotone scheme for the equation under consideration in . The above second order generalization is not guaranteed to be a consistent monotone scheme for the equations considered here, but it works nicely in practice.

## References


format: text

text:
# Bubble dynamics in DNA

(March 19, 2024)

## Abstract

The formation of local denaturation zones (bubbles) in double-stranded DNA is an important example for conformational changes of biological macromolecules. We study the dynamics of bubble formation in terms of a Fokker-Planck equation for the probability density to find a bubble of size $`n`$ base pairs at time $`t`$, on the basis of the free energy in the Poland-Scheraga model. Characteristic bubble closing and opening times can be determined from the corresponding first passage time problem, and are sensitive to the specific parameters entering the model. A multistate unzipping model with constant rates recently applied to DNA breathing dynamics \[G. Altan-Bonnet et al, Phys. Rev. Lett. 90, 138101 (2003)\] emerges as a limiting case.

: J. Phys. A: Math. Gen.

Introduction. Under physiological conditions, the equilibrium structure of a DNA molecule is the double-stranded Watson-Crick helix. At the same time, in essentially all physiological processes involving DNA, for docking to the DNA, DNA binding proteins require access to the “inside” of the double-helix, and therefore the unzipping (denaturation) of a specific region of base pairs . Examples include the replication of DNA via DNA helicase and polymerase, and transcription to single-stranded DNA via RNA polymerase. Thus, double-stranded DNA has to open up locally to expose the otherwise satisfied bonds between complementary bases.

There are several mechanisms how such unzipping of double-stranded DNA (dsDNA) can be accomplished. Under physiological conditions, local unzipping occurs spontaneously due to fluctuations, the breathing of dsDNA, which opens up bubbles of a few tens of base pairs . These breathing fluctuations may be supported by accessory proteins which bind to transient single-stranded regions, thereby lowering the DNA base pair stability . Single molecule force spectroscopy opens the possibility to induce denaturation regions of controllable size, by pulling the DNA with optical tweezers . In this way, the destabilising activity of the ssDNA-binding T4 gene 32 protein has been probed, and a kinetic barrier for the single-strand binders identified . Denaturation bubbles can also be induced by under-winding the DNA double helix . A recent study of the dynamics of these twist-induced bubbles in a random DNA sequence shows that small bubbles (less than several tens of base pairs) are delocalised along the DNA, whereas larger bubbles become localised in AT-rich regions . Finally, upon heating, dsDNA exhibits denaturation bubbles of increasing size and number, and eventually the two strands separate altogether in a process called denaturation transition or melting . Depending on the specific sequence of the DNA molecule and the solvent conditions, the temperature $`T_m`$ at which one-half of the DNA has denatured typically ranges between $`50^{}`$C and $`100^{}`$C. Due to the thermal lability of typical natural proteins, thermal melting of DNA is less suited for the study of protein-DNA interactions than force-induced melting . On the other hand, the controlled melting of DNA by heating is an important step of the PCR method for amplifying DNA samples , with numerous applications in biotechnology .

The study of the bubble dynamics in the above processes is of interest in view of a better understanding of the interaction with single-stranded DNA binding proteins. This interaction involves an interplay between different time scales, e.g., the relaxation time of the bubbles and the time needed for the proteins to rearrange sterically in order to bind . Dynamic probes such as single-molecule force spectroscopy  and molecular beacon assays  may therefore shed light on the underlying biochemistry of such processes.

In a recent experiment by Altan-Bonnet et al, the dynamics of a single bubble in dsDNA was measured by fluorescence correlation spectroscopy. It was found that in the breathing domain of the DNA construct (a row of $`18`$ AT base pairs sealed by more stable GC base pairs) fluctuation bubbles of size $`2`$ to $`10`$ base pairs are formed below the melting temperature $`T_m`$ of the AT breathing domain. The relaxation dynamics follows a multistate relaxation kinetics involving a distribution of bubble sizes and successive opening and closing of base pairs. The characteristic relaxation time scales were estimated from the experiment to within the range of $`20`$ to $`100\mu `$s. Also in reference , a simple master equation of stepwise zipping-unzipping with constant rate coefficients was proposed to successfully describe the data for the autocorrelation function, showing that indeed the bubble dynamics is a multistate process. The latter was confirmed in a recent UV light absorption study of the denaturation of DNA oligomers .

In this work, we establish a general framework to study the bubble dynamics of dsDNA by means of a Fokker-Planck equation, based on the bubble free energy function. The latter allows one to include microscopic interactions in a straightforward fashion, such that our approach may serve as a testing ground for different models, as we show below. In particular, it turns out that the phenomenological rate equation approach, corresponding to a diffusion with constant drift in the space of bubble size $`n`$ used by Altan-Bonnet et al to fit their experimental data corresponds to a limiting case of our Fokker-Planck equation. However, the inclusion of additional microscopic interactions in such a rate equation approach is not straightforward . In what follows, we first establish the bubble free energy within the Poland-Scheraga (PS) model of DNA denaturation , and then derive the Fokker-Planck equation to describe the bubble dynamics both below and at the melting temperature of dsDNA.

Bubble free energy. In the PS model, energetic bonds in the double-stranded, helical regions of the DNA compete with the entropy gain from the far more flexible, single-stranded loops . The stability of the double helix originates mainly from stacking interactions between adjacent base pairs, aside from the Watson-Crick hydrogen bonds between bases. In addition, the positioning of bases for pairing out of a loop state gives rise to an entropic contribution. The Gibbs free energy $`G_{ij}=H_{ij}TS_{ij}`$ for the dissociation of two paired and stacked base pairs $`i`$ and $`j`$ has been measured, and is available in terms of the enthalpic and entropic contributions $`H_{ij}`$ and $`S_{ij}`$ . In the following we consider a homopolymer for simplicity, as suitable for the AT breathing domain in reference . For an AT-homopolymer ($`i=j=\mathrm{AT}`$), the Gibbs free energy per base pair in units of $`k_BT`$ yields as $`\gamma \beta G_{ii}/2=0.6`$ at $`37^{}`$C for standard salt conditions $`(0.0745\text{M-Na}^+)`$. Similarly, for a GC-homopolymer one finds the higher value of $`\gamma =1.46`$ at $`37^{}`$C. The condition $`\gamma =0`$ defines the melting temperature $`T_m`$ , thus $`T_m(\mathrm{AT})=66.8^{}`$C and $`T_m(\mathrm{GC})=102.5^{}`$C (we assume that $`G_{ij}`$ is linear in $`T`$, cf. reference ). Above $`T_m`$, $`\gamma `$ becomes negative. For given $`\gamma =\gamma (T)`$, the statistical weight for the dissociation of $`n`$ base pairs obtains as

$$W(n)=\mathrm{exp}(n\gamma ).$$
(1)
Additional contributions arise upon formation of a loop within dsDNA. Firstly, an initial energy barrier has to be overcome, which we denote as $`\gamma _1`$ in units of $`k_BT`$. From fitting melting curves to long DNA, $`\gamma _110`$ was obtained, so that the statistical weight for the initiation of a loop (cooperativity parameter), $`\sigma _0=\mathrm{exp}(\gamma _1)`$, is of order $`10^5`$ . As the energy to extend an existing loop by one base pair is smaller than $`k_BT`$, DNA melts as large cooperative domains. Below the melting temperature $`T_m`$, the bubbles become smaller, and long range interactions beyond nearest neighbours become more important. In this case, the probability of bubble formation is larger, and $`\gamma _1`$ ranges between 3 and 5, thus $`\sigma _00.05`$  (cf. section 5 in ). According to reference , the smallness of $`\sigma _0`$ inhibits the recombination of complementary DNA strands with mutations, making recognition more selective. Secondly, once a loop of $`n`$ base pairs has formed, there is a weight $`f(n)`$ of mainly entropic origin, to be detailed below. The additional weight of a loop of $`n`$ open base pairs is thus

$$\mathrm{\Omega }(n)=\sigma _0f(n).$$
(2)
For large bubbles one usually assumes the form $`f(n)=(n+1)^c`$ . Here, the value of the exponent $`c=1.76`$ corresponds to a self-avoiding, flexible ring , which is classically used in denaturation modelling within the PS approach. Recently, the PS model has been considered in view of the order of the denaturation phase transition . Reference  finds by finite size scaling analysis of measured melting curves of DNA oligomers that the transition is of second order. In reference , the value $`c=2.12`$ was suggested, compare the discussion in references . For smaller bubbles (in the range of 1 to a few tens of base pairs), the appropriate form of $`f(n)`$ is more involved. In particular, $`f(n)`$ will depend on the finite persistence length of ssDNA (about eight bases), on the specific base sequence, and possibly on interactions between dissolved but close-by base pairs (cf. section 2.1.3 in ). Therefore, the knowledge of $`f(n)`$ provides information on these microscopic interactions.

For simplicity, we here adopt the simple form $`f(n)=(n+1)^c`$ for all $`n>0`$, and consider the loop weight

$$\mathrm{\Omega }(n)=\sigma _0(n+1)^c.$$
(3)
We show that at the melting temperature the results for the relaxation times for the bubbles are different for the available values $`c=1.76`$ and $`c=2.12`$ quoted above. This shows that the specific form of $`f(n)`$ indeed modifies experimentally accessible features of the bubble dynamics .

In what follows, we focus on a single bubble in dsDNA, neglecting its interaction with other bubbles. Since due to $`\sigma _01`$ the mean distance between bubbles ($`1/\sigma _0`$ ) is large, this approximation is justified as long as the bubbles are not too large . It also corresponds to the situation studied in the recent experiment by Altan-Bonnet et al. According to the above, the total free energy $`(n)`$ of a single bubble with $`n>0`$ open base pairs follows in the form

$$\beta (n)=\mathrm{ln}\left[W(n)\mathrm{\Omega }(n)\right]=n\gamma (T)+\gamma _1+c\mathrm{ln}(n+1),$$
(4)
where the dependence on the temperature $`T`$ enters only via $`\gamma =\gamma (T)`$. We show the free energy (4) in figure 1 for $`c=1.76`$ and for the parameters of an AT-homopolymer, for physiological temperature $`T=37^{}`$C, at the melting temperature $`T_m=66^{}`$C, and at $`T=100^{}`$, compare the discussion below.

Bubble dynamics. In the generally accepted multistate unzipping model, the double strand opens by successively breaking Watson-Crick bonds, like opening a zipper . As $`\gamma `$ becomes small on increasing the temperature, thermal fluctuations become relevant and cause a random walk-like propagation of the zipper locations at both ends of the bubble-helix joints. The fluctuations of the bubble size can therefore be described in the continuum limit through a Fokker-Planck equation for the probability density function (PDF) $`P(n,t)`$ to find at time $`t`$ a bubble consisting of $`n`$ denatured base pairs, following a similar reasoning as pursued in the modelling of the dynamics of biopolymer translocation through a narrow membrane pore . To establish this Fokker-Planck equation, we combine the continuity equation (compare reference )

$$\frac{P(n,t)}{t}+\frac{j(n,t)}{n}=0$$
(5)
with the expression for the corresponding flux,

$$j(n,t)=D\left(\frac{P(n,t)}{n}+\frac{P(n,t)}{k_BT}\frac{}{n}\right)$$
(6)
where it is assumed that the potential exerting the drift is given by the bubble free energy (4). In equation (6), we incorporated an Einstein relation of the form $`D=k_BT\mu `$, where

the mobility $`\mu `$ has dimensions $`[\mu ]=\mathrm{sec}/(\mathrm{g}\mathrm{cm}^2)`$, and therefore $`[D]=\mathrm{sec}^1`$ represents an inverse time scale. By combination of equations (4), (5) and (6), we retrieve the Fokker-Planck equation for $`P(n,t)`$ 111Note that the operator $`\frac{}{n}`$ acts also on $`P(n,t)`$.

$$\frac{P(n,t)}{t}=D\left(\frac{}{n}\left\{\gamma +\frac{c}{n+1}\right\}+\frac{^2}{n^2}\right)P(n,t).$$
(7)
Thus, we arrived at a reduced 1D description of the bubble dynamics in a homopolymer, with the bubble size $`n`$ as the effective ‘reaction’ coordinate. For a heteropolymer, also the position of the bubble within the double helix, i.e., the index $`m`$ of the first open base pair, becomes relevant. In this case, the bubble free energy $``$ and thus the PDF depend on both $`m`$ and $`n`$, and on the specific base sequence; the corresponding generalisation of equation (7) is straightforward. In a random sequence, additional phenomena may occur, such as localisation of larger bubbles . Finally, to establish the Fokker-Planck equation (7), we assume that changes of the bubble size $`n`$ occur slower than other degrees of freedom of the PS free energy within the bubble region (e.g., Rouse-Zimm modes). This assumption seems reasonable due to the long bubble dynamics’ relaxation time scales of $`20`$ to $`100\mu `$s , and the good approximation of bubble independence .

By rescaling time according to $`tDt`$, the Fokker-Planck equation (7) can be made dimensionless, a representation we are going to use in the numerical evaluation below. The formulation in terms of a Fokker-Planck equation makes it possible to derive the characteristic times for bubble closing and opening in terms of a first passage time problem. That is, for bubble closing, the associated mean closing time follows as the mean first passage time to reach bubble size $`n=0`$ after starting from the initial bubble size $`n_0`$. We now determine these characteristic times for the three regimes defined by $`\gamma `$ with respect to temperature $`T`$.

(i) $`T<T_m`$. In this regime, the drift consists of two contributions, the constant drift $`D\gamma `$ and the loop closure component $`Dc/(n+1)`$, which decreases with $`n`$. For large $`n`$, we can therefore approximate the drift by the constant term $`D\gamma `$, and in this limit the Fokker-Planck equation (7) is equivalent to the continuum limit of the master equation used in reference  to describe the experimental bubble data. In particular, we can identify our two independent parameters $`D`$ and $`\gamma `$ with the rate constants $`k_+`$ and $`k_{}`$ to open and close a base pair introduced in , respectively: $`D(k_++k_{})/2`$ and $`\gamma 2(k_{}k_+)/(k_++k_{})`$. In this approximation, the correlation functions used successfully to fit the experimental results in  can be derived from the Fokker-Planck equation (7). Moreover, we can deduce the mean first passage time PDF for a bubble of initial size $`n_0`$ to close in the exact analytical form (compare )

$$f(0,t)=\frac{n_0}{\sqrt{4\pi Dt^3}}\mathrm{exp}\left\{\frac{(n_0D\gamma t)^2}{4Dt}\right\},$$
(8)
which decays exponentially for large $`n`$. In particular, from (8) the characteristic (mean) first passage time for bubble closing,

$$\tau =n_0/(D\gamma )$$
(9)
follows, which is linear in the initial bubble size $`n_0`$. In figure 2, we compare this analytical result for the value $`\gamma (37^{}\text{C})`$ with the characteristic closing times using the full drift term from equation (7), obtained from numerical integration. It can be seen that the qualitative behaviour for both cases with and without the $`Dc/(n+1)`$ term is very similar, but that in the presence of the loop closure component the characteristic times are reduced.

(ii.) $`T=T_m`$. At the melting temperature, the drift in equation (7) is solely given by the loop closure term $`Dc/(n+1)`$. In figure 2, we plot the characteristic bubble closing time obtained numerically. In comparison to the above case $`T<T_m`$, the faster than linear increase as function of initial bubble size $`n_0`$ is distinct. Keeping in mind that $`Dc/(n+1)`$ becomes very small for increasing $`n`$, this behaviour can be qualitatively understood from the approximation in terms of a drift-free diffusion in a box of size $`n_0`$, in which the initial condition $`P(n,0)=\delta (nn_0)`$ is located at a reflecting boundary, and at $`n=0`$ an absorbing boundary is placed. This problem can be solved analytically, with the result $`\tau =n_0^2/(2D)`$ for the characteristic escape time 222This result can, for instance, be obtained through the method of images, compare reference . in which the quadratic dependence on $`n_0`$ contrasts the linear behaviour in the result (9). Thus, at the melting temperature $`T_m`$, the tendency for a bubble to close becomes increasingly weaker for larger bubble sizes, and therefore much larger bubbles can be formed, in contrast to the case $`T<T_m`$. A further comparison to the value $`c=2.12`$ mentioned above demonstrates that a clear quantitative difference in the associated closing times exists, see figure 2. However, the qualitative behaviour remains unchanged. In principle, the study of the bubble dynamics can therefore be used to discern different models for the loop closure factor.

(iii.) $`T>T_m`$. Above the melting temperature, the drift is governed by the interplay between the loop closure component $`Dc/(n+1)`$ tending to close the bubble, and the bubble free energy $`D\gamma (T>T_m)<0`$, which causes a bias towards bubble opening. In figure 1, we show for the AT-homopolymer case how the overall drift potential after a small initial activation barrier becomes negative, and the dynamics is essentially governed by the $`\gamma `$-contribution. As a consequence, from the result (8), we find that the associated mean first passage time diverges, i.e., the bubble on average increases in size until the entire DNA is denatured, as expected from the PS model.

By symmetry, similar results hold for the bubble opening process. However, the existence of the bubble initiation energy $`\gamma _1`$ involves an additional Arrhenius factor, which is not included in the Fokker-Planck equation (7), and which reduces the opening probability, causing an increase of the associated opening time, cf. reference .

Conclusions. From the bubble free energy of a single, independent bubble in the Poland-Scheraga theory of DNA melting, we derived a Fokker-Planck equation for the PDF to find a bubble created by denaturation of $`n`$ base pairs at time $`t`$. This formulation allows for the calculation of the characteristic times scales for bubble closing and opening in terms of first passage time problems. Three different regimes were distinguished: (i) below the melting temperature, the characteristic bubble closing time increases linearly in bubble size, and the drift can be approximated by the constant value $`D\gamma `$. In this approximation, the Fokker-Planck equation matches the continuum version of the master equation employed previously in an experimental study of DNA breathing . (ii) At the melting temperature, an approximately quadratic growth of the bubble closing time as a function of bubble size is observed, which can be explained by noting that the $`1/(n+1)`$-dependence of the loop closure drift component can be neglected for larger $`n`$, leading to pure diffusion. In this approximation, the exact analytical results indeed lead to the quadratic dependence observed in the numerical results. (iii) Above the melting temperature, the characteristic closing time in our model diverges, consistent with the fact that on average the DNA follows the trend towards the thermodynamically favourable state of complete denaturation.

The expression for the Gibbs free energy used in our approach involves a purely entropic contribution for a single-stranded bubble. It was suggested in reference  that also in denaturation bubbles a residual stacking energy $`\epsilon _s`$ is present. In this case, the value of $`\gamma `$ would have to be corrected by this $`\epsilon _s`$.

From the Fokker-Planck equation (7), in which the microscopic interactions enter via the free energy (4), one can derive measurable quantities such as moments and dynamic correlation functions . In principle, the Fokker-Planck equation involves one free parameter, the time scale $`1/D`$, while the values for the other parameters are known. However, by fitting sufficiently accurate experimental data for DNA bubble dynamics at different temperatures, the values for additional parameters may be extracted by help of our general Fokker-Planck framework probing suitable free energy functions.

We are happy to acknowledge helpful discussions with Terence Hwa, Yacov Kantor, Mehran Kardar, Richard Karpel, Udo Seifert, and Mark Williams. We also thank Oleg Krichevsky for sending us a preprint of reference  prior to publication.


text:
# Untitled Document

Phase transitions, memory and frustration in a

Sznajd-like model with synchronous updating

Lorenzo Sabatelli1,2 and Peter Richmond1

1 Department of Physics, Trinity College Dublin, Ireland

2 HIM, IFSC, Dublin, Ireland

E-mail: .sabatell@tcd.ie

We introduce a consensus model inspired by the Sznajd Model. The updating is synchronous and memory plays a decisive role in making possible the reaching of total consensus. We study the transition between the state with no-consensus to the state with total consensus.

Key words: computer simulation, opinion dynamics, phase transitions

Basic features of Sznajd Model

It is well understood, within human societies, that it is generally easier to change someone’s opinion by acting within groups than by acting alone. For example, a single person stopping on the street and staring at the sky is usually ignored (or perhaps considered eccentric). However, if several people stare into the sky, they readily induce others to do the same.

Equally it is not unusual to come across door-to-door sale agents working in couples rather than on an individual basis. In democratic electoral campaigns, small groups (typically couples) of political party activists may visit, door-to-door, potential electors and seek to gain their votes. Trade union movements often try to coordinate their actions in order to strengthen their position against management than if everybody try to negotiate it alone.

The underlying principle, in all of the above-mentioned examples may be captured by the famous Abraham Lincoln’s injunction: ”United we stand, divided we fall”. This principle was developed into a computational model by Katarzyna Sznajd and her father.

The simplest (non-trivial) version of their model can be implemented on a two-dimension lattice of spins. Each site carries a spin, $`S`$ , that is either up or down. This represents one of the two possible opinions on any question. Two neighbouring parallel spins, i.e. two neighbouring people sharing the same opinion, convince their neighbours of this opinion. If they do not have the same opinion, then they do not influence their neighbours.

The system evolves from one time step to another through a random sequential updating mechanism.

The system always reaches an overall consensus, within a sufficiently long time and ends up with all spins up (or down) if the initial fraction of up opinions is larger (smaller) than $`1/2`$ . Furthermore, the smaller the size of the lattice, the smoother is the transition.

If the random sequential updating is replaced by a synchronous updating mechanism the possibility of reaching a consensus is reduced quite dramatically.

The updating is performed by going systematically through the lattice to find the first member of the pair, then choosing randomly the second member of the pair within the neighborhood of the first. Having in this way completed the assembly of couples, each agent then orients her/himself according to her/his neighbours at time step $`t`$ . Like-minded couples will induce their neighbours to turn to the same state (opinion). However a single agent may often belong, simultaneously, to the neighbourhood of more than one couple (of likeminded agents). In this case, if the couples have different opinions, she/he doesn’t know what to do (frustration) and ends up doing nothing, i.e. sticks with her/his previous opinion. Frustration may prevent the system from reaching total consensus.

Memory

A feature of the agents in the models discussed above is the complete absence of memory. The past plays no role.

In this note we assume that agents are endowed with memory. For the sake of simplicity, they are all thought to have the same memory span $`T`$ and updating mechanism is synchronous. Agents are keen to change opinion when in the neighbourhood of a like-minded couple, as in the models above. An agent resorts to her/his individual history when frustration occurs. In that case, the new state turns out to be the most frequent of her/his own $`T+1`$ most recent ($`T`$ accounts for the past, one for the present). Of course, this rule proves to be more efficient when $`T`$ is an even number.

Phase transitions

In a two-dimensional Sznajd model with asynchronous updating, the system always ends up at a fixed point: either all spins point up or they all point down. When synchronous updating is chosen things are different. The system converges to a fixed point (all spins up or all spins down), only if the asymmetry $`\mathrm{\Delta }p`$ (absolute difference) in the initial distribution of opinions is above a critical value $`\mathrm{\Delta }p_c`$ that depends not only on the size of the lattice  but also on the memory length $`T`$ \[fig.1\].

Fig.2 shows 1-$`\mathrm{\Delta }p_c`$ as a function of $`L`$ (the lattice side size). The slope decreases as $`T`$ increases. For small values of $`L`$ the three curves displayed seem to converge to the same point. This is easily understood because for small lattices one expects the impact of frustration on total consensus to be small. As a consequence, the role of the memory $`T`$ is not as important as it would be for large values of $`L`$ .

For large $`L`$ the curves monotonically decay to zero, this suggests that asymptotically (for infinitely large lattice) no consensus is possible.

Fig.3 shows $`\mathrm{\Delta }p_c`$ as a function of $`T`$ . The slope of the two curves appears to be independent of $`L`$ \- an aspect that requires further study.

Conclusions

We have studied a model based on the Sznajd consensus model with synchronous updating that includes memory. This feature plays a very important role and helps overcome frustration and achieve total consensus.

The memory length, $`T`$ , affects the phase transition from the no-consensus state to the total consensus state. The bigger $`T`$ the closer to zero is the critical point. For a given side size $`L`$ and $`T`$ enough large an always-consensus situation, as with the random sequentially updated Sznajd model, may be attained. Anyway, for infinitely large (but with finite memory) lattices, no consensus is possible, just as in the case without memory.

Acknowledgments

We would like to thank Dietrich Stauffer who provided continual advice, stimulating remarks and good humour during this study.

The authors acknowledge support from the EU via Marie Curie Industrial Fellowship MCFH-1999-00026

References

 K. Sznajd-Weron and J. Sznajd. Opinion Evolution in Closed Community. Int. J. Mod. Phys. C 11: 1157-1166 (2000).

 D. Stauffer. Monte Carlo Simulations of the Sznajd model, Journal of Artificial Societies and Social Simulation 5, No.1, paper 4 (2002) (jasss.soc.surrey.ac.uk).

 D. Stauffer. Frustration from Simultaneous Updating in Sznajd Consensus Model. cond-mat/0207598

Fig. 1.1

The phase transition from a no-consensus state to a total consensus state is driven by the parameter $`\mathrm{\Delta }`$p for $`L`$ =17 and $`T`$ =0, 2, 8. The transition point is shifted towards zero as the agent memory length $ T$ increases.

Fig. 1.2

The phase transition from a no-consensus state to a total consensus state is driven by the parameter $`\mathrm{\Delta }`$p for $`L`$=101 and $`T`$=0, 2, 8. The transition point is shifted towards zero as the agent memory length $`T`$ increases.

Fig. 1.3

The phase transition from a no-consensus state to a total consensus state is driven by the parameter $`\mathrm{\Delta }`$p for $`L`$=301 and $`T`$=0, 2, 8 The transition point is shifted towards zero as the agent memory length $`T`$ increases.

Fig. 1.4

The phase transition from a no-consensus state to a total consensus state is driven by the parameter $`\mathrm{\Delta }`$p for $`L`$=1001 and $`T`$=0, 2, 8. The transition point is shifted towards zero as the agent memory length $`T`$ increases.

Fig. 1.5

The phase transition from a no-consensus state to a total consensus state is driven by the parameter $`\mathrm{\Delta }`$p for $`L`$=50 and $`T`$=0, 2, 8, 100, 500, 1000. The transition point is shifted towards zero as the agent memory length $`T`$ increases.

Fig. 2

Variation with $`L`$ of 1-$`\mathrm{\Delta }`$p (one minus the absolute value of the difference in the initial probabilities for +1 and –1) for which in half of the cases a consensus was reached. That may be seen as the phase transition point from the state without consensus to the state with consensus. The estimated slope is –0.39 for $`T`$=0, -0.21 for $`T`$=2, - 0.11 for $`T`$=8.

Fig. 3

Variation with $`T`$ of the difference in the initial probabilities for which in half of the cases a consensus was reached. That may be seen as the phase transition point from the state without consensus to the state with consensus. The estimated slope is 0.46 for $`L`$=17 and 0.47 for $`L`$=50.