text:
# Adaptive Event Horizon Tracking and Critical Phenomena in Binary Black Hole Coalescence

## I Introduction

There is a two-fold motivation for this work. On one hand, as argued in , , , and , a robust event horizon solver provides an abundance of intuition and information in numerical analysis of the binary black hole coalescence problem. Alternatively, several open questions concerning the dynamics of black hole event horizons in the strong field domain remain unanswered in the literature. In this work, these motivations are pursued in development of a robust event horizon solver that is capable of establishing an intriguing analogy between black holes undergoing merger and fluid droplets undergoing bifurcation.

A theorem due to Penrose (cited in ) establishes that the event horizon of a black hole is generated by null geodesics with no future end point. Specifically,

1. Followed into the past, a null geodesic can only leave the event horizon at a special event denoted a caustic.
2. Through each non-caustic event of the event horizon there is a unique null geodesic.

As demonstrated below, these properties of the generators yield a rich structure for the dynamical evolution of the event horizon.

For example, it can be shown by way of the no hair theorem that the topology of sections of the event horizon for a stationary black holes is necessarily spherical. It was believed, at first, that in the dynamical strong field regime similar behavior persists. In a surprising result, Hughes et al., and others found numerically a momentarily toroidal section of the event horizon in the nonstationary case of the gravitational collapse of a ring of particles . Subsequent work by Shapiro, Teukolsky and Winicour  and others , , ,  has shown that this result is correct and is currently conjectured to be generic for asymmetric gravitational collapse . Briefly, in an exactly future asymptotically stationary spacetime the horizon generators can be shown to merge in the asymptotic past at a zero dimensional point, which is unstable to perturbation of this zero dimensional locus of generators. Thus, horizons that are dynamical in the asymptotic future yield higher dimensional loci — curves and surfaces — of merged generators. Consequently, there exists a higher genus horizon for some interval in the evolution. Numerical and exact studies have presented similar evidence in  leading to the conjecture that event horizons of black holes are unrestricted in their genus .

Related work on the phenomenology of black hole event horizons concerns the differentiability of the two dimensional membrane of the event horizon sections , , . Creases and caustics of the event horizon, where the horizon is no longer differentiable, serve an important role both in the topology of single black hole event horizons and in the merger of multiple black holes into a single horizon. Where the event horizon of a black hole undergoes a change in topology, the surface of section becomes strictly singular . Thus, where two or more sections of the event horizons merge to form a single horizon each surface of section is singular and no longer differentiable at the points of merger , . These special crease and caustic points of the horizon are currently unrestricted in their presence on the horizon. For example, it can be shown  that working only with the definition of a black hole event horizon in terms of null geodesic generators having no future end point, it is possible to construct event horizons that are nowhere differentiable. Physically this can be approximated by a “cloud of sand” falling into a large black hole: Each “grain” falling into the hole generates a caustic that terminates where the grain crosses the horizon.

Both of these considerations are problematic to numerical studies — particularly ones that seek a generic solution or work from the finite difference approximation. These interesting phenomenological aspects of the caustics and crease sets of black hole event horizons undergoing merger or other topological transitions are an open field of research.

Here we concentrate on the development of a generic numerical black hole event horizon tracking method. It is possible to generically solve this numerical problem, and the solution is presented in . The computational technique uses the finite difference approximation, thus assumes a degree of differentiability which will not be achieved at caustics, creases, etc. Nevertheless, we find that the approach provides very good approximation in those cases; e.g., gives surfaces very close (and convergently close) to analytical expectations. In particular, problems associated to differentiability of the horizon that arise due to creases and caustics of the horizon can be handled on a case by case basis. One such approach, advocated here, is the incorporation of adaptive mesh refinement (AMR). Incorporating AMR into a numerical method allows special points, such as crease sets or caustics, to be resolved to within machine epsilon — a resolution that is typically sufficient to determine the global behavior of the phenomena. But it remains that special attention is required if very fine scales are to be resolved.

A subject related both to the need for AMR in numerically tracking black hole event horizons and to the presence of crease or caustic points in continuum black hole event horizons, is the question of power law scaling of black hole event horizons undergoing topological transitions. To motivate the possibility of this effect, it is useful to consider the ‘membrane paradigm’ of black hole event horizons. In this approach the event horizon sections are considered in analogy to a two dimensional fluid, such as the surface of a liquid droplet and in fact, the dynamics of the sections of the horizon have been shown to obey evolution equations analogous to those of a two dimensional fluid . Accordingly, the event horizon of a black hole can be viewed as a distinct dynamical physical entity within a spacetime. Further, the membrane paradigm equations were demonstrated for numerically detected black hole event horizons in  and, as demonstrated in later sections, offer insight into the dynamics governing the numerical evolution of the critical membrane.

The formal analogy between black hole event horizons and the surface of a fluid droplet, spelled out by the membrane paradigm, suggests that the critical phenomena associated to fluid droplet bifurcation are also present in the binary black hole coalescence problem. Study of the fluid problem typically establishes power law scaling of the throat using both numerical and exact analysis of the Navier Stokes equation governing the dynamics; in almost all of those studies, use is made of AMR . The successes of the membrane paradigm and fluid studies of droplet bifurcation suggest that AMR applied to black hole event horizon solvers will produce analogous black hole results.

In section II we present a new AMR method for tracking black hole event horizons. We call our method the comoving front tracking method and notate it here as cmft. In section II the specifics of our method are shown in detail by building on the work of , which addressed the problem of numerically tracking black holes as a computational front tracking problem. The accuracy of an implementation of our method is shown in sections II. A. and B. which consider the case of Kerr - Newman black holes with and without coordinate deformations of the source. Section III applies our method to the Kastor - Traschen solution , describing the merger of two charged black holes in a spacetime with cosmological constant. In that section we study the minimal radius of the neck connecting the two black holes immediately following merger. We find, in analogy to the case of fluid droplets undergoing bifurcation, that the minimal radius of the neck undergoes power law scaling with the minimal throat radius proportional to time. In section IV, we summarize and discuss our conclusions.

## II Comoving Front Tracking

Adaptive mesh refinement is a numerical method developed generically for hyperbolic systems . The method was first applied in numerical relativity, with considerable success, by Mathew Choptuik, in study of the problem of black hole formation under the gravitational collapse of massless scalar fields . The method embodies Brandt’s rule of numerical analysis: Computational resources should be applied proportionally to the physical processes involved . As the method’s name implies, adaptive mesh refinement (AMR) is typically applied where solution or field variables require variable resolution over the computational mesh: Some sectors of the grid require higher resolution (where the field variables undergo interesting changes), while other sectors undergo less activity and so need not be computed to high accuracy. There are two distinct approaches to AMR. In one approach, there is one grid of coarsest resolution that is fixed but within which higher resolution grids are nested according to the behavior of the solution. In the alternate method of AMR, which is the one used here, there is one grid that moves with the solution. Such methods have been applied with great success to other fields of computational physics, particularly in problems involving the motions of fronts, such as those found across phase transitions.

To understand the cmft method for tracking black hole event horizons, it is useful to first consider the method as applied in a fixed mesh to the problem of tracking black hole event horizons in , although those studies did not consider the application of adaptive meshes. Since the event horizon is a null surface (away from caustics) one must study the null geodesic equations; we solve the null geodesic equations by solution of an eikonal equation. In the front tracking approach a clever coordinate system, adapted to the black hole’s surface of section $`\mathrm{\Gamma }`$, is chosen. Let $`\left\{\sigma _i\right\}`$ be such a coordinate system. Then the front $`\mathrm{\Gamma }`$ of a single level set of $`S`$, where $`S`$ solves the eikonal equation

$$g^{ab}_aS_bS=0,$$
(1)
can be tracked by elimination of one coordinate; vis,

$$S(x^i,t)=\sigma _1u(\sigma _1,\sigma _2,t)=0.$$
(2)
In terms of these coordinates, the eikonal equation written in the ADM variables

$$_tS=\beta ^i_iS\pm \sqrt{_jS\gamma ^{jk}_kS}$$
(3)
becomes

$$_tu=\beta ^1+\beta ^I_Iu\pm \alpha \sqrt{\gamma ^{11}+2\gamma ^{1I}_Iu+_Iu\gamma ^{IJ}_Ju}.$$
(4)
Here $`I=2,3`$ and $`\gamma ^{IJ}`$ is the two dimensional metric obtained by the choice of coordinates.

Several comments are relevant. First, this method considerably extends the non - linearity of the eikonal equation: The geometric variables $`\alpha ,\beta ^I,\gamma ^{IJ}`$ are each functions of the coordinates $`x^i`$ and consequently, the geometric variables are also functions of the grid function $`u`$; e.g., $`x^i=x^i\left(\sigma ^k\right)x^i(u,\sigma ^I)`$. That is, whereas the eikonal (3) was a hyperbolic equation of the general nonlinear form

$$_tS=F(t,x,_xS),$$
(5)
equation (4) is a hyperbolic equation of the general form

$$_tu=G(t,x,u,_xu).$$
(6)
As a direct consequence, while the nonlinearity of the gradient of $`u`$ is given in terms of the root in (4), the nonlinearity in terms of $`u`$ typically cannot be classified; particularly when the geometric variables are only provided numerically. For example, in numerically generated spacetimes, a case for almost all problems of interest, the geometric variables $`\alpha ,\beta ,\gamma ^{IJ}`$ are defined on a global grid $`𝒢`$ of $`N^3`$ points

$$x^i\left(k^i\right)=s^ik^i+x^i\left(0\right),$$
(7)
where the integer $`k^i`$ satisfies $`1k^iN`$ for each $`i`$. The front tracking method then generically requires interpolation methods since the two - dimensional mesh $`G^0`$ used for the finite difference approximation of the two - dimensional surface of sections of the horizon $`\mathrm{\Gamma }`$ will not in general coincide with the points of the global grid $`𝒢`$. Update of the grid function $`u`$ then requires knowledge of the geometric variables $`\alpha ,\beta ,\gamma ^{IJ}`$ at points that do not lie on $`𝒢`$, which can only be found by interpolation. A further complication of the front tracking method relates to the question of the coordinates $`\left\{\sigma _i\right\}`$. As described in  the method requires different coordinate systems depending on the physical processes involved. \[In fact, the original studies  found dependencies and sensitivities of results on the choice of coordinates, which is an unappealing feature of the method, particularly from the viewpoint of general relativity.\] For example, in the case of a single hole with static topology a spherical coordinate system is typically sufficient. However, for the case of head on binary black hole merger a cylindrical coordinate system was employed , while it is unclear what global coordinates should be chosen, or if there is even one set of coordinates which can be used for the problem of two black holes undergoing inspiral to merger.

The method of adaptive mesh front tracking developed here addresses several of the complications of the method described above. First, within our approach we fix the choice of coordinates as spherical coordinates although the cmft method is not contingent on this choice. Instead, spherical coordinates are chosen due to their boundary conditions, which are advantageous in implementation. Note that while spherical coordinates clearly cannot handle evolution into any topologies beyond the genus zero, $`S^2`$ topology of stationary black hole event horizons, the cmft method compensates for this by allowing the mesh to adaptively track black hole event horizons that are stationary in their asymptotic past and future; more importantly, through refinement, to naturally detect the onset of topology change, where (strictly) the surface becomes singular. In such circumstances, the surface can be monitored to within arbitrary proximity of the transition; and then continued past the transition by applying the code individually to the resulting black holes.

To make the method precise, at some fixed time level $`t=t_n`$ let $`c^i=x^i`$ be the average value of $`M^2`$ points distributed over a surface $`\mathrm{\Gamma }`$ having $`S^2`$ topology and satisfying a level set condition $`S\left(\mathrm{\Gamma }\right)=0`$. That is, in a two dimensional mesh $`G^n`$ of $`M^2`$ points $`x^i=x^i(I,J)`$ on $`\mathrm{\Gamma }`$ where the integers $`I,J`$ are $`1I,JM`$ let

$$c^i=x^i\frac{1}{M\left(M1\right)}\underset{I,J=1}{\overset{I,J=M,M1}{}}x^i(I,J).$$
(8)
Local coordinates on the surface $`\mathrm{\Gamma }`$ can then be written

$$xc^1=r\mathrm{cos}\left(\varphi \right)\mathrm{sin}\left(\theta \right)$$
(9)
$$yc^2=r\mathrm{sin}\left(\varphi \right)\mathrm{sin}\left(\theta \right)$$
(10)
$$zc^3=r\mathrm{cos}\left(\theta \right).$$
(11)
According to this choice $`r=u(\theta ,\varphi ,t_n)`$ can be updated in the two dimensional mesh $`G^n`$ according to a finite difference representation of equation (2) expressed in the spherical coordinates That is, the form of (2) is chosen here as

$$r=u(\theta ,\varphi ,t)=\sqrt{\left(xc^1\right)^2+\left(yc^2\right)^2+\left(zc^3\right)^2}.$$
(12)
The coordinates $`x,y,z`$ are defined both in terms of the global grid $`𝒢`$ and the $`M^2`$ mesh $`G^n`$ having center $`c^i`$. By comparison, the grid function $`u`$ is defined locally on the comoving mesh $`G^n`$. In the case that the geometric variables $`\alpha ,\beta ^i,\gamma ^{ij}`$ are provided numerically from solution on the global $`N^3`$ grid $`𝒢`$, interpolation of those variables onto the local $`M^2`$ surface $`\mathrm{\Gamma }`$ is required for update of $`u`$. This is a generic feature of the front tracking method. Further, in our method in one cycle of the construction, the radial coordinate $`u`$ will evolve in accordance with the update since the surface center $`c^i\left(t\right)`$ will also update: $`c^i\left(t\right)c^i\left(t+dt\right)`$. This update of the surface center corresponds to creation of a new $`M^2`$ mesh $`G^nG^{n+1}`$ and therefore we need a numerical change of variables $`u_{IJ}^{n+1}u_{I^{}J^{}}^{n+1}`$ between the two meshes. This change of variables will always, irrespective of the nature of the geometrical variables, require interpolation of the grid function $`u`$ from the one grid, $`G^n`$, to the other, $`G^{n+1}`$, since the points of those grids do not coincide in general. This procedure can become in practice a very detailed and intricate computational step since it amounts to passing data from one spherical mesh to another spherical mesh and the procedure can potentially fall prey to grid tangling effects. Most of the intricacy of the data passing $`G^nG^{n+1}`$ is restricted to the choice of spherical coordinates. For the purposes of this work ordinary second order interpolation proves sufficient both for grid passing and for interpolation of the geometric variables onto the surface. Finally, the iterated Crank Nicholson scheme is generically well suited as a finite difference approximation for the partial differential equation (4). Studies were conducted with other schemes, such as the method of lines used by other researchers in the front tracking problem, but superior performance was found with the iterated Crank Nicholson method. Discussion of this finite difference approximation can be found in .

In summary, a pseudo code expression for one complete update of the comoving geometry is then:

Pseudo-Code: A Complete Update Iteration

* Load $`\left\{x_{IJ}^n\right\}`$
* Build $`c^n`$ and $`u_{IJ}^n`$
* Interpolate $`\gamma `$ onto $`\mathrm{\Gamma }`$
* Update $`u_{IJ}^nu_{IJ}^{n+1}`$
* Update $`c^nc^{n+1}`$
* Pass data $`u_{IJ}^{n+1}u_{I^{}J^{}}^{n+1}`$
* Return $`\left\{x_{I^{}J^{}}^{n+1}\right\}`$

### A Kerr-Newmann Black Holes

In this section the accuracy of an implementation of the cmft method is considered in detail. For more detailed discussion of the signatures of black hole event horizons in the eikonal equation see  and . Briefly, since the event horizon of a black hole is a global structure of spacetime, its detection cannot be determined without the complete history of the Cauchy evolution of the spacetime. However, as shown below, the horizon is a critical outgoing null surface that neither expands to infinity nor collapses into the gravitational singularity. Numerical event horizon solvers have typically employed this property as the signature behavior of event horizon detection and tracking. In these approaches, the space of outgoing null surfaces propagated into the past is surveyed for evidence of the critical behavior of the horizon. Since outgoing null data followed into the past tend to approach the horizon to great precision these methods typically produce highly accurate approximations of the event horizon itself. For example, in figure (1) we show outgoing data propagated into the past in a background spacetime contatining a spherically symmetric black hole. In that figure three classes of data are considered: Null data initially interior to the black hole event horizon, null data initially exterior to the black hole event horizon and null data that is exactly on the black hole event horizon. In all cases, we see that outgoing null data propagated into the past approaches the black hole event horizon to high accuracy.

To establish the accuracy of our implementation, we consider first the case of stationary, spinning black holes, which is completely described by the Kerr - Newmann axisymmetric solutions of Einstein’s equation. According to Carter’s theorem , this family of solutions is the unique, asymptotically flat, stationary and axisymmetric black hole solutions of the vacuum equations. As such, it is sufficient to consider this class of solutions in an account of stationary black hole event horizons. It is convenient to the discussion to make use of the Kerr-Schild form for the metric :

$$g^{ab}=\eta ^{ab}2Hl^al^b.$$
(13)
Here $`\eta _{ab}=\mathrm{diag}(1,1,1,1)`$ is Minkowski’s metric, $`H`$ is a space time scalar, and $`l^a`$ is an ingoing null vector with respect to both the Minkowski and full metric. The Kerr solution is the two parameter family of solutions such that

$$H=\frac{Mr^3}{r^2+a^2z^2}$$
(14)
and

$$l^t=1,$$
(15)
$$l^x=\frac{rx+ay}{r^2+a^2},$$
(16)
$$l^y=\frac{ryax}{r^2+a^2},$$
(17)
$$l^z=\frac{z}{r},$$
(18)
$$r^2=\frac{1}{2}\left(\rho ^2a^2\right)+\sqrt{\frac{1}{4}\left(\rho ^2a^2\right)+a^2z^2}$$
(19)
where

$$\rho ^2=x^2+y^2+z^2.$$
(20)
Two points are worth mentioning. Firstly, the parameter $`M`$ corresponds to the gravitational mass of the source, while the parameter $`a`$ corresponds to the source’s spin: It can be shown that observers at asymptotic infinity measure the angular momentum of the source to be $`J=aM`$. Secondly, the cartesian coordinates are chosen such that the $`z`$ axis is aligned with the direction of spin and so is an axis of symmetry.

The Kerr family of solutions is interesting from the perspective of tracking black hole event horizons since the solutions actually contain two event horizons; with one nested interior to the other. Furthermore, the solutions have a ring curvature singularity located at $`r=z=0`$. In the cmft method, provided the event horizon is approached from the exterior, neither of these features requires special attention.

It is analytically known that the outermost horizon of a spinning black hole is located at the surface

$$r_+=M+\sqrt{M^2a^2}$$
(21)
with an area

$$A=4\pi \left(r_{+}^{}{}_{}{}^{2}+a^2\right).$$
(22)
We will now validate our cmft codes’ ability to determine these values.

We start with nonspinning black holes. With $`a=0`$, figure (1) shows the signature of the black hole event horizon for outgoing null data propagated backwards in time. Figure (2) similarly shows this effect for outgoing data and includes variation in the mass with the cases of $`M=1/2,1,2`$ considered.

In all cases interior and exterior data converge exponentially to the event horizon as predicted by the exact solution. The $`e`$ \- folding time of this exponential behavior (which can be shown in perturbation theory to satisfy $`\gamma =\left(4M\right)^1`$ for ingoing Eddington - Finkelstein coordinates) is shown numerically in figure (3). The slope of each line is $`\left(4M\right)^1M`$ in agreement with the perturbation prediction.

We turn now to consideration of spinning black holes. The ellipsoidal geometry of sections of spinning black hole event horizons are demonstrated in figures (4), (5), (6), (7), which show the intersection of the event horizon section both with a $`\varphi =\pi `$ slice and with a $`\theta =\pi /2`$ slice. These figures are the asymptotic limit cross sections of spheres tracked backwards until the area achieves stationarity. They are thus approximations of the black hole event horizons. The figures are drawn for $`a/M=0.25,0.5,0.75,`$ and $`0.9`$. As can be seen there is a slight asymmetry (an error) of the resulting surfaces. The accuracy of these results is shown in figure (8), which contains the L2 norm of the numerical truncation error $`e_r`$ of the $`r_+`$ function calculated over the surface of the final state. With $`\widehat{r}_+`$ denoting the finite difference approximation of the continuum funtion $`r_+`$, the truncation error is defined to be

$$e_r=\widehat{r}_+r_+.$$
(23)
In figure(8) time in units of $`M`$ is measured increasing into the past.

Figure (9) shows several values of the percent error in the area for a fixed unit mass black hole with spin parameter $`0a0.9`$. These figures were generated using a numerical resolution of $`m^2`$ with $`m=50`$ on the detected event horizon surface section<sup>*</sup><sup>*</sup>*Throughout the studies in this work an area calculation algorithm first proposed in  is used. In this method the metric $`h`$ of the space time is determined on the surface $`\mathrm{\Gamma }`$ and the area is calculated using numerical integration of $`A=𝑑\theta 𝑑\varphi \sqrt{h}`$. Also, note carefully that lower case is a number, giving the resolution of the surface discretization. M is a mass.. The truncation error of the numerical integration of the area scales at least as well as $`𝒪\left(h^2\right)`$ as required for second order convergence. This scaling of the truncation error is shown in figure (10), which shows the percentage error of the area of a nonspinning black hole for resolutions of $`m^2`$ with $`m=25,50,100`$.

In terms of the consistency of these results for the truncation error of $`r_+`$ and the percent errors of the calculated areas, note that there are at least two sources of systematic error in the numerical calculation of the area of any surface. One source of error is the accuracy with which each point of the surface is known, while the second source of error is the finite resolution of the discrete version of the surface. Let $`A`$ denote the continuum area and $`\widehat{A}`$ be the discrete version of the surface calculated using a resolution of $`m^2`$. Then by the finite difference approximation

$$A=\widehat{A}+𝒪\left(1/m^2\right)$$
(24)
for second order convergence. However, if the points of the surface are systematically inaccurate $`\widehat{A}`$ will not converge to $`A`$ in the continuum limit; instead $`\widehat{A}`$ will converge to the bias of the implementation. This result suggests that the bias in $`\widehat{A}`$ is on the order of $`0.1`$ percent.

The $`e`$-folding time for relaxation of outgoing null data onto the event horizon, or equivalently for formation of the solution singularity in the eikonal, is shown in figure (11) with various values of the spin $`a/M=0.25,`$ $`0.5,`$ $`0.75,`$ $`0.9`$. Note that only for the cases of rapidly spinning black holes $`a/M0.9`$ does the $`e`$-folding time differ appreciably from the spin zero case of $`\gamma =1/4M`$. These results suggest using the ‘rule of thumb’ relaxation time of $`\gamma 1/4M`$ for any spin with $`a/M0.9`$.

### B Non-linear Coordinate Deformations

The method of comoving adaptive meshes makes use of an intricate numerical procedure: The passing and interpolation of data from one spherical grid to another. To ensure that the computational implementation of the method is both accurate and robust, consider the case of the following coordinate transformation:

$$tt^{}=t,$$
(25)
$$xx^{}=x+b\mathrm{cos}\omega t,$$
(26)
$$yy^{}=y+b\mathrm{sin}\omega t,$$
(27)
$$zz^{}=z$$
(28)
A nonspinning black hole located at $`x=y=z=0`$, is then seen to rotate with angular frequency $`\omega `$ at a radius of $`r=b`$ about the point $`x^{}=y^{}=z^{}=0`$. The physics of this coordinate transformation is trivial and the area of the source event horizon remains $`A=16\pi M^2`$. However, from the perspective of the numerical implementation written in the coordinates $`(t^{},x^{},y^{},z^{})`$, the black hole appears highly dynamical and moving with an angular velocity. As such, the coordinate transformation $`(t,x,y,z)(t^{},x^{},y^{},z^{})`$ applied to the Kerr - Newmann family of analytic sources is a stringent test of an implementation of the method of comoving adaptive meshes.

Consider first the case of a linear shift of coordinates: $`\omega =0,b0`$. Figure (12) shows the percent error in the area of the detected horizon versus time increasing under propagation into the past, while figure (13) shows the L2 norm of the error in the function $`r_+`$ for $`b=3/2,1,1/2,0`$. Both of these figures are generated using surfaces of $`m^2`$ points with $`m=25`$ and spherical initial data of radius $`r=2`$ centered at the origin. The source is that of a nonspinning unit mass black hole. The solution ultimately settles down to the correct shift zero result, which suggests that this implementation of the method is stable with respect to translation of the source.

Additional complexity is obtained by fixing $`b=1`$ and varying $`\omega =2\pi /T`$. Figure (14) shows the percentage error in the area of the detected horizon section, while figure (15) shows the L2 norm of the error in $`r_+`$. In these figures the surface used was that of $`m^2`$ points with $`m=25`$, and spherical initial data of radius $`r=3`$ located at the origin. Also here, the source considered in both figures is a unit mass black hole with spin parameter $`a=m/2`$. According to these results, this implementation of the method of comoving adaptive meshes converges to the bias of the surface resolution as $`\omega 0`$. Further, for time scales on the order of the relaxation time $`4M`$, the error of the implementation is fairly significant although it remains below $`10\%`$. These result suggests that the particular implementation is sensitive to coupling of the time scale of the event horizon’s dynamics. Further accuracy can be obtained by considering a more stringent implementation.

## III Symmetric Binary Black Hole Coalescence

The cmft method cannot continuously monitor a topological transition in the event horizon of a black hole. However, with high resolution or refinement the method can detect the onset of a topological transition and come arbitrarily close to the transition itself.

In this context we now consider the Kastor - Traschen analytic solution of the Einstein - Maxwell $`Q=M`$ equations with cosmological constant $`\mathrm{\Lambda }=3H^2`$  . The solution is simply

$$ds^2=U^2dt^2+U^2\left(dx^2+dy^2+dz^2\right)$$
(29)
where

$$U=Ht+\frac{M_1}{r_1}+\frac{M_2}{r_2}.$$
(30)
$`H=\pm \sqrt{\mathrm{\Lambda }/3}`$ and we set the coordinates and choose the sign so that $`H=1`$. These are the choices made in ref . Both refs.  and  show that with these choices the black hole merger occurs in the future (i.e. as $`tt`$ increases into the future). The sign convention differs between refs.  and ; we follow .

Here also, $`M_1`$ and $`M_2`$ denote the masses of the two holes located at $`a_1`$ and $`a_2`$ with

$$r_{1,2}=\sqrt{x^2+y^2+(za_{1,2})^2}.$$
(31)
We consider the equal mass case ($`M_1=M_2=0.1`$), symmetrically oriented ($`a_1=a_2=a=0.1`$) on the $`z`$-axis. These are the parameters chosen in a calculation by Siino. The choice of relatively small $`M`$ keeps the final (postmerger) single black hole horizon inside the cosmological (deSitter) horizon. Ref.  demonstrates that for small negative $`t`$ there is one component of the event horizon that encompasses the origin; at somewhat earlier times there are two separated components of the event horizon. Figure 16 shows the cmft front tracker result (running backward from a late time t but still with $`t<0`$ ), with an “initial” guess equal to the late time single apparent horizon. Refs  and  estimate the late time single horizon at the apparent horizon:

$$r_{bh}=R_{bh}/(Ht)=(12MH\alpha )/(2tH^2),$$
(32)
where

$$\alpha =\sqrt{14MH},$$
(33)
and this $`M=M_1+M_2`$. ($`R_{bh}`$ is a constant, though $`r_{bh}`$ is not.)

The coalescence occurs at a small negative time estimated in  as: $`t^{}R_{bh}/(2aH)<0`$, which for our parameters is $`t^{}=0.763932`$. Following coalescence the merged black hole settles down to the final state of a charged black hole in a spacetime with cosmological constant. Figure (16) shows several frames of a movie indicating bifurcation of the black hole event horizon when tracked backwards in time from the outermost, guess surface at the latest time. Note that the solution will eventually break down when the grid function $`u`$ becomes multiple valued. Note also that near the bifurcation the solution exhibits multiple time scales. For example, away from the throat the solution appears quasi-stationary, while at the throat the solution is clearly dynamic.

As we investigate the behavior of the horizon found by cmft near breakdown, we can look for evidence of power law scaling of the topological transition. In analogy to the topological transitions found in the bifurcation (“pinch off”) of fluid droplets, such power law scaling would be expected in the neck of the event horizon at the point and time of merger, or “pinch on”   (the opposite of “pinch off”). Using cylindrical coordinates, let $`\rho `$ denote the radius of the throat connecting the two black holes. By the axial symmetry of the solution, $`\rho =\rho (z,\varphi ,t)=\rho (z,t)`$. At the “pinch on” time, $`t=t^{}`$, $`\rho (0,t^{})=0`$. If the solution exhibits power law scaling then

$$\rho =C\left(tt^{}\right)^\gamma $$
(34)
for times after the “pinch on”.

Figure (17) shows an order of magnitude in the evolution of $`\rho `$ just after the topological transition, using the cmft code. For $`0.2tt^{}0.6`$ the radius of the throat appears to vary exponentially, while for $`tt^{}0.0`$ the radius appears to vary linearly. Figure (18) shows that behavior of $`\rho `$ is indeed purely exponential in $`t`$ far from the topological transition. Further, figure (19) (also found using the cmft code) shows evidence near the transition for power law (in fact, linear) scaling of the throat radius. Note that without arbitrarily high resolution or nested adaptive mesh refinement any scaling that is present in the solution cannot be resolved to machine epsilon. Therefore, the results presented here using the cmft code can at best establish qualitative evidence for such scaling. Figure (19) does indeed present such evidence that the solution exhibits power law scaling in analogy to the bifurcation of fluid droplets. Further, according to the results of the cmft code, it is here conjectured that the Kastor - Traschen solution scales like $`\gamma 1`$. These results indicate that the cmft method is capable of probing into the fine structure of black hole event horizons undergoing a topological transition. Using the results from the cmft code, with $`100^2`$ points on the tracked surface, we obtain the computational estimate $`t^{}=1.454.`$

This result is somewhat different from the analytical estimates of Ref. . However, we can verify this result because we can carry out a much more accurate direct (ordinary differential equation) integration to find the merge time. We use the equatorial symmetry of this equal mass case merger, for which the $`\rho t`$ behavior is simply understood. In the equal mass case the throat is circular in cross section, and this cross section lies in the plane of symmetry of the problem. Suppose we introduce a cylindrical coordinate system for the metric (29). The event horizon is defined by outgoing null rays, and the throat thus is generated by null rays moving in the $`\rho `$ direction (the throat moves at the speed of light):

$$\frac{d\rho }{dt}=U^2$$
(35)
where

$$U=Ht+2M/\sqrt{\rho ^2+a^2}.$$
(36)
Figure (20) shows a log vs log plot of $`\rho `$ vs $`tt^{}`$, again begun from the late time single apparent horizon. It indicates a critical exponent of exactly $`\gamma =1`$. This plot was generated using an adaptive, Runge Kutta , numerical integration of the ordinary differential equation for $`\rho \left(t\right)`$ to drive the difference $`t^{}t`$ to below machine epsilon. The $`t^{}`$ found in this simulation is $`t^{}=1.45845468374`$, in agreement within accuracy to the cmft result.

Straightforward physics explains this answer. The instant of merger $`t^{}`$ is not a special value for the Kastor - Traschen metric function $`U`$. Hence, near this event, equation (35) for $`\rho `$ is simply

$$\frac{d\rho }{dt}=U^2\left(t^{}\right)=\mathrm{const}0$$
(37)
plus higher order terms. Thus, in particular, the geometrically significant quantity, the circumference of the throat, $`C(t)U(t^{})2\pi \rho (t)`$ grows linearly with time, $`\gamma =1`$.

## IV Conclusions

We have argued here that adaptive mesh refinement applied to the front tracking approach to tracking black hole event horizons is a potentially fruitful method; and moreover proves necessary for black hole processes in which the event horizon undergoes either a change of topology or develops creases or caustics. Towards that end we have developed an adaptive mesh technique that makes use of comoving meshes and demonstrated the accuracy of our approach. Applying our method to the symmetric collision data of the Kastor Traschen analytic solutions, we have found the surprising result of a power law in the minimal radius of the throat connecting the black holes following merger. We have shown power law scaling for the throat diameter at merger of black holes in a special situation: $`\mathrm{\Lambda }=3(0)`$, $`M_1=M_2=0.1`$, $`Q=M`$, and axisymmetic. This effect is in analogy to that of fluid droplet problems. Since it is generically true that the black hole event horizon satisfies dynamical evolution equations that are exactly the form of a two dimensional fluid, we conjecture and give some qualitative arguments below concerning the specializations in the Kastor Traschen to argue that the power law scaling found in the Kastor Traschen solution must be generic for the symmetric collision of two black holes. Extending the analogy further, we also expect a similar power law for the asymmetric problem, although generalized in that case to accomodate the asymmetry of the problem as we now argue.

It is clear that $`\mathrm{\Lambda }0`$ is irrelevant to this phenomena. The scaling behavior applies to an arbitrarily thin neck that first connects merging black holes. The curvature of the neck scales as $`1/(\rho \times \mathrm{separation})\left(\rho M\right)^1`$. For small $`\rho M`$ this dominates the fixed contribution.

The assumption of extremal black holes ($`Q=M`$) is a strong one. In Newtonian terms, it implies no acceleration of the motion between the holes (the electrical force is of equal magnitude and opposite the gravitational force). A first intuition may be that this has something to do with the exponent $`\gamma `$ in the scaling of the throat radius ($`\gamma =1`$ both numerically and analytically). However, the generic axisymmetric merger of black holes will deal with the scaling in a very brief interval at the beginning of merger; the “relative velocities” of the holes will be essentially constant in any case. We thus predict that extremality ($`Q=M`$) is not essential to the $`\gamma =1`$ power law behavior found here. For similar reasons we do not expect the spin of holes (zero in the case of the Kastor - Traschen solution) to affect the results found here. \[The special (non-generic) case where non-extremal equal mass hole data are set with the holes just at merger will be different, because we expect the holes to “accelerate from zero” as the throat grows.\]

Still in the context of axisymmetry, we can investigate whether the assumption of equal mass is relevant to the value of $`\gamma `$. In the limit of $`M>>m`$ we essentially consider a null topological tube (the small hole $`m`$) merging with a null plane $`M`$. We are investigating this configuration for merger. We again expect power law scaling; the question is to determine the exponent. If in this case $`\gamma 1`$, then $`\gamma `$ will in general be a function of the ratio of the masses: $`\gamma =\gamma (m/M)`$; $`\gamma (1)=1`$. We no longer have the symmetry plane to simplify analysis; the throat cannot be followed by following the individual $`\rho `$ direction null trajectories, but Lindblom  has given an argument, based on the smoothness of the throat after its formation, that gives $`\gamma =1`$ for the growth of the circumference in this case also.

Finally, we consider relaxing axisymmetry. In analysis and simulations based on perturbation of late time horizons evolved backwards, toroidal behavior is sometimes seen in the non-axisymmetric case. The horizons grow more then one “point” as they approach merger; two points that touch simultaneously in very quick succession will create an evanescent toroidal horizon. To date, there has been no study directly relating parameters of the toroid to parameters of the colliding holes (we have seen only suggestions of toroidal behavior in any of our computations (c.f. the “offset points” of the event horizon in , Figure 11) , perhaps because of resolution limitations). However, we would expect each merging “point” to behave like a seperate merger and thus to undergo power law scaling, but not necessarily with $`\gamma =1`$. In fact, analysis by Winicour  in the generic case shows that the toroid closes superluminally; thus prevents sending signals “through the hole”. This requires that the “points” grow superluminally, which requires $`\gamma <1`$. The crucial difference from the case explained by Lindblom, is that the toroidal structure is not described by a smooth surface. In particular, there is an inner “crease” where geodesics are joining the torus, allowing it to close superluminally.

## V Ackmowledgment

We thank L. Lindblom for very insightful and helpful comments. This work was supported by NSF rant PHY0102204. Computations were done at the University of Texas Advanced Computing Center.


format: text

text:
# The urban economy as a scale-free network

## I Introduction

The Zipf rank-size law for city sizes is one of the most widely known power laws in scienceZipf\_1949 . It is also but one out of many similar power laws from systems in biology, economy and society. We continue this research by presenting empirical evidence that land values are scale-free. The data we use is based on a database delivered by Sweden Statistics that covers estimations of the market value of all land in Sweden (2.9 million data points).

Although power laws are common they are not easily reconstructed from realistic underlaying dynamics. Their ubiquity suggests that they could be caused by some general systemic property common to a range of systems. Recent research suggests that many empirically observed power laws may be due to fundamental properties of these systems viewed as networks of interacting nodesBarabasi\_Albert\_1999 ; Albert\_Barabasi\_2002a ; Dorogovtsev\_Mendes\_2002a . We investigate the mechanisms causing land values to follow these statistics and present a network model that reproduces the empirical results. The model is based on basic definitions of city formation in urban economics theoryOSullivan\_2002 .

According to urban economics theory, the formation of modern cities is primarily caused by the advantages of trade between specialized producers. The exchange of goods and services between localized and largely immobile activities in trade economies makes a network representation natural: the nodes are units of land and the edges represent the exchange of goods and services between them.

It has been shown that the node degrees of a certain class of growing networks are power law distributedBarabasi\_Albert\_1999 ; Dorogovtsev\_Mendes\_Samukhin\_2000 . This class of networks is important because their growth mechanisms can be mapped to the microscopic dynamics of several real-world systems. We demonstrate how trade in an urban system can be represented as a scale free network and that, as a consequence, land values can be expected to follow the same distribution. We also verify that the model retains these properties when spatial constraints are taken into account. We do this by using a spatial network model to reproduce the empirically observed distribution of land values.

The network approach solves a fundamental issue in the problem of modelling urban systems by representing the system at the level of the underlaying market structure. This allows us to produce prices in units of currency rather than undefined and subjective fitness measures. It thereby opens up doors for several extensions of the scope of the model and provides a natural interface for integrating it with other models.

When we refer to an urban system we do not necessarily refer to individual cities but rather to systems of specialized trading activities. To clarify further, our use of the term *activity* refers to trade gains in units of *currency per unit area and unit time*. Activities can be resolved to any resolution down to individual transactions.

## II The urban economy as a network

In section II.1 we define a non-spatial model of urban economic growth and in section II.2 we extend it to a spatial model where growth is mapped to a 2D surface. In section II.3 we motivate the model ontology and the basic assumptions on which we have based the model. In particular we discuss the connection between node degree and land value and how the urban system meets the criteria for being a scale-free network.

### II.1 Formulation of the non-spatial model

A geographic area on which an urban economic system can grow is represented by an enumerated set of nodes, $`\{1,2,\mathrm{},N\}`$ corresponding to non-overlapping land areas. Trade of goods and service between activities in the nodes is represented by undirected edges. Since activities within the same node can trade with each other, an edge can connect a node with itself. The amount of activity of a site $`x_i`$ is defined as the degree of the corresponding node. See figure 1.

The network is initialized by connecting $`n_0`$ nodes so that each has a degree of $`x_0`$. The $`Nn_0`$ undeveloped nodes have no trade interactions and thus no activity.

At each time step we update the network as follows:

1. With probability $`q_1`$ we add $`m`$ edges between sites that already are developed: the first edge end-point is selected uniformly among developed nodes. The probability of a node $`i`$ to be selected is

$$\mathrm{\Pi }_i^u=\delta _i^{(D)}\frac{1}{n_t^{(D)}},$$
(1)
where $`n_t^{(D)}`$ is the number of developed nodes at time $`t`$ and $`\delta _i^{(D)}=1`$ if node $`i`$ is developed and $`\delta _i^{(D)}=0`$ otherwise.

The second end-point is selected preferentially among developed nodes. Preferential selection corresponds to the uniform selection of an activity in the system and the subsequent location of its node. It was defined by Barabasi et al Barabasi\_Albert\_1999  as

$$\mathrm{\Pi }_i^p=\frac{x_i}{_jx_j},$$
(2)
where $`\mathrm{\Pi }_i^p`$ is the probability of node $`i`$ to be attached to a new edge and $`x_i`$ is the degree of node $`i`$.
2. With probability $`q_2`$ we add edges between $`m`$ pairs of nodes that are both selected preferentially according to Eq. (2).
3. With probability $`q_3`$ we add $`m`$ units of initial activity on land that is previously undeveloped: the first end-point is selected without bias similarly to Eq. (1) but among undeveloped nodes. However, since nodes have no properties beside their degree in the non-spatial model, any undeveloped node can be added. The second end-point is selected preferentially according to Eq. (2).

We will refer to the growth classes as type-1, type-2 and type-3 growth and to their relative rates as $`q_1`$, $`q_2`$ and $`q_3`$ with $`q_1+q_2+q_3=1`$ throughout the paper.

#### II.1.1 A continuum formulation

Instead of assuming any particular set of activities and interactions we can use a continuum formulation where we consider each edge end-point to be an average of a large set of urban activities. It follows that an average activity must be assumed to interact equally much with all other average activities. Rather than counting explicit interactions to determine the activity level, we study the evolution of the expected node degrees. The time evolution of activity on a developed site $`i`$ follows the equation

$$x_i(t+1)=x_i(t)+q_1\frac{m}{n_t^{(D)}}+m(q_1+2q_2+q_3)\frac{x_i(t)}{_jx_j(t)},$$
(3)
which is solved by the continuous-time method introduced by Barabási et al Albert\_Barabasi\_Jeong\_1999 . After sufficiently long time the degree distribution approaches the form

$$P[x_i=x](x+A)^\gamma ,$$
(4)
with

$$A=\frac{2mq_1}{q_3(1+q_2)}$$
(5)
and

$$\gamma =1+\frac{2}{1+q_2}.$$
(6)
According to Eq. (4) the node degree distribution will be power law distributed for the non-spatial model.

### II.2 Formulation of the spatial model

The non-spatial model does not include any information about inhomogeneous relations in the network. An important departure from such a simple model that can be readily observed in real systems is of course the fact that different pairs of sites have different distances between them. Because of transportation cost optimization this has the potential to affect the edge distribution of the network model.

We will now extend the basic model to incorporate spatial interactions. This allows us to verify that the power law distribution of activities that is predicted by the non-spatial model is retained when space is introduced. Furthermore, it allows us to better map model output to empirical measurements. Scale-free spatial network models that have been studied recently are not directly applicable to urban economics since they require an *a priori* distribution of node locationsBarthelemy\_2002 ; Brunet\_Sokolov\_2002 ; Manna\_Sen\_2002 .

#### II.2.1 Network evolution

In the continuum non-spatial model we were able to treat network evolution as a purely local phenomenon since interactions in the network are homogenous, see Eq. (3). In the spatial model we have to again explicitly separate the dynamics of the edge end-points. This is because the spatial context of the first end-point of a new edge modifies the probabilities with which other sites will become the second end-point (transportation costs).

The same types of growth that are used in the non-spatial model are also used in the spatial model. But, as stated, the two end-points here need to be dealt with separately. The first site (primary increase) is selected either uniformly or preferentially depending on growth type. This will be discussed shortly. The second end-points of the trade relation edges (secondary increase) are always preferential.

#### II.2.2 Secondary increase

The spatial model tells us how a given primary increase $`p_j`$ in activity at site $`j`$ causes secondary increases $`s_i`$ in activity at all other sites $`i`$ where

$$s_i=p_j\frac{D_{ij}x_i}{_kD_{kj}x_k},$$
(7)
This is analogous to Eq. (2) but in the spatial model, transportation costs will bias the choice of trade partners. $`D_{ij}`$ is a matrix representing the interaction strength, which is assumed to decay with the transportation costs. For the results in this paper we have used $`D_{ij}=(1+cd(i,j))^\alpha `$, where $`d(i,j)`$ is the Euclidean distance between sites $`i`$ and $`j`$. The non-negative parameters $`c`$ and $`\alpha `$ controls the impact of spatial. Another function that can apply is $`D_{ij}=max(0,1d(i,j))`$ if the transportation characteristics of the activity is known, which can be the case in a model where activity types are modelled separately rather than as an average. Exponential $`D_{ij}=e^{d(i,j)}`$ can also be interesting to the extent that shielding is important, i.e. an activity tends to trade exclusively with the nearest supplier.

#### II.2.3 Primary increase

When trade takes place there is a mutual benefit in efficiency that is often used for further increasing the activity in the city. This feedback process makes it possible for the amount of activity in the urban system to increase considerably faster than the population. It is useful to think of Eq. (7) as a black box system to which a driving force, primary increase, is applied. It should also be noted that the primary effects used for the model we present here are by no means neither exhaustive nor final: most earlier models of urban growth could be introduced as primary effects in our framework.

#### II.2.4 Primary effects in type-1, type-2 and type-3 growth

Network evolution in the spatial model is similar to the non-spatial model. Secondary increases are always preferential following Eq. (7) and primary increases now reflect the spatial distribution:

Type-1 growth: the primary uniform increase is identical to the non-spatial case following Eq. (1).

Type-2 growth: the primary preferential increase is identical to the non-spatial case following Eq. (2).

Type-3 growth is split into two related processes where one corresponds to growth in the perimeter of clusters and the other corresponds to growth in connection to infrastructure in the rural areas between clusters. Such infrastructure is not explicitly represented in our model and instead we use a parameter $`ϵ`$ to tune the amount of ambient infrastructure and thus the rate with which seemingly isolated clusters will appear.

Type-3a growth: with a probability of $`q_3(1ϵ)`$ we set the activity of a perimeter node to $`m`$. Perimeter nodes are nodes that are not developed but borders to at least one developed cell. The site of the new node is selected randomly and with uniform probability among the perimeter sites on the grid

$$\pi _i^a=\delta _i^{(P)}\frac{1}{n_t^{(P)}},$$
(8)
where $`\pi _i^a`$ is the probability with which node $`i`$ is selected to undergo type-3a growth, $`\delta _i^{(P)}=1`$ if the node $`i`$ is on the perimeter, $`\delta _i^{(P)}=0`$ if the node is not on the perimeter and $`n_t^{(P)}`$ is the number of perimeter nodes at time $`t`$.

Type-3b growth: With a probability of $`q_3ϵ`$ we set the activity of an external node to $`m`$. An external node is a node that is undeveloped and that has no developed neighbors. The site of the new node is selected randomly and with uniform probability among external sites on the grid

$$\pi _i^b=\delta _i^{(E)}\frac{1}{n_t^{(E)}},$$
(9)
where $`\pi _i^b`$ is the probability with which node $`i`$ is selected to undergo type-3b growth, $`\delta _i^{(E)}=1`$ if the node $`i`$ is external, $`\delta _i^{(E)}=0`$ if the node is not external and $`n_t^{(E)}`$ is the number of external nodes at time $`t`$.

#### II.2.5 A continuum formulation of the spatial model

A continuum formulation for the evolution of developed nodes in the spatial model can, as in section II.1.1, be constructed by studying the time evolution of expected node degrees.

$$x_i(t+1)=x_i(t)+E[p_i(t)]+E[s_i(t)],$$
(10)
where $`E[s_i]`$ is the expected secondary increase which can be calculated by a weighted summation over the expected primary increments,

$$E[s_i]=\underset{j}{}E[p_j]\frac{D_{ij}x_i}{_kD_{kj}x_k}.$$
(11)
For the special case of our simple model for the primary effects we have

$$E[p_j]=\delta _j^{(D)}(p_j^{(1)}+p_j^{(2)})+\delta _j^{(P)}p_j^{(3a)}+\delta _j^{(E)}p_j^{(3b)},$$
(12)
with

$`p_j^{(1)}=q_1{\displaystyle \frac{m}{n_t^{(D)}}},`$ (13)
$`p_j^{(2)}=q_2{\displaystyle \frac{mx_j}{_kx_k}},`$ (14)
$`p_j^{(3a)}=q_3(1ϵ){\displaystyle \frac{m}{n_t^{(P)}}},`$ (15)
$`p_j^{(3b)}=q_3ϵ{\displaystyle \frac{m}{n_t^{(E)}}}.`$ (16)
The total growth for a developed node can in this case be separated into a uniform and a preferential part as

$$x_i(t+1)=x_i(t)+\zeta _i(t)+\eta _i(t)x_i(t),$$
(17)
with

$$\zeta _i=q_1\frac{m}{n_t^{(D)}},$$
(18)
and

$$\eta _i=\frac{q_2m}{_kx_k}+\underset{j}{}E[p_j]\frac{D_{ij}}{_kD_{kj}x_k}.$$
(19)
Noting that $`_i\eta _ix_i=m(q_1+2q_2+q_3)`$, equation 17 can be rewritten as

$$x_i(t+1)=x_i(t)+\zeta _i(t)+m(q_1+2q_2+q_3)\frac{\eta _i(t)x_i(t)}{_j\eta _j(t)x_j(t)},$$
(20)
which, in a comparison with equation (3), reveals that the only difference between the the spatial and the non-spatial model is the site and time-dependent parameter $`\eta _i`$. This is similar to the concept of node fitness, as presented in Bianconi\_Barabasi\_2001a ; Ergun\_Rodgers\_2001 , which can affect the node degree distribution. However, our simulation results indicate that $`\eta _i`$ falls within a sufficiently narrow interval for the power law to be essentially preserved (Fig. 2). This is also supported by calculations of $`\eta _i`$ for both simulated and empirical data.

### II.3 The network model in an urban economics context

#### II.3.1 The connection between node degree and land value

An approximate linear relationship between node degree in the model and land value in the real system is crucial for the interpretation of our results. The motivation follows from i) market pricing of goods and services and ii) the connection between trade benefits and land value.

i) Market pricing of commodities provides an adaptive measure that allows us to compare the activities that generate them. Hence, on average, an edge contributes identically to the value of both nodes to which it connects. This contribution is exactly our definition of activity, which implies that the degree of a node is proportional to its benefits due to trade.

ii) This connection consists of two proportionalities. For a node $`i`$ we have

$$v_ir_ix_i,$$
(21)
where $`v_i`$ is the value of the corresponding land area, $`r_i`$ is the bid-rentAlonso\_1960 ; Alonso\_1964 , and $`x_i`$ is the total trade benefits as outlined above. Capitalizing periodic rent income from the site $`i`$ gives land value $`v_i=\frac{r_i}{i}`$ where $`i`$ is the interest rateOSullivan\_2002 . The second proportionality is a weak form of the leftover principle from urban economics, which states that, in a competitive land market, rent equals the amount of money leftover after all expenses (except rent) are paid. This amount of money equals the sum of all trade benefits at the site. For our results it is sufficient that, on average, a certain proportion of each new unit of trade benefit goes to the landowner.

Together, i) and ii) suggest an approximately linear relationship between node degrees and land values.

#### II.3.2 Types of growth

Urban activity can increase in essentially two ways: either a new activity is related to, or it is unrelated to, an existing activity at the site. In the former case (preferential growth) this could be a new employee hired as a response to increased demand, in the latter case (uniform) it could be the establishment of a new firm. Preferential growth corresponds to a per-unit activity rate. Uniform growth corresponds to establishment among lots on a competitive land market where, for the average land use, we can not expect any lot to be more profitable than any other.

#### II.3.3 Reasons for treating perimeter growth separately from internal growth

The jagged perimeter of urban areas exposes large areas of undeveloped land to urban infrastructure, thus making it attractive for urban land use. Because perimeter land is in ample supply and currently has a low revenue, even land uses with a very low trade gain can be competitive. Many low-activity land uses in the outskirts of the urban area can likely just barely out-bid agriculture and would not be viable in competition with other urban land uses. Among high-activity land uses some have very specific demands on land improvements and can therefore not benefit from buying existing buildings inside the urban area. This creates a special case for perimeter land. Note that just like for type-1 growth, competition prevents prediction of where the next growth event will take place among the perimeter nodes, and no preferential growth is possible since there is no previous activity that can expand.

#### II.3.4 Node fitness and growth

If we only regard trading activities the only difference between two sites with identical activity is the value of their spatial context. Therefore, in the spatial model, secondary growth is not homogenous, see Eq. (7). Node fitness (Eq. 17) can be viewed as a local interest rate that predicts the growth rate of activity investments made in that site. The result of this is non-trivial growth predictions since the expected amount of new local development does not become a simple fraction of current development.

Most notably, the model predicts the emergence of urban sub-centers. This is realized by examining the expected secondary effects in the spatial model (Eq. 11). Apart from being proportional to the amount of present activity $`E[s_i]`$ is also subject to site competition and will increase for nodes that have high activity in relation to their own neighborhoods. For each possible primary effect in $`j`$, the node development $`x_i`$ is weighted with the fraction between trade intensity between site $`j`$ and $`i`$ and the sum of the all other trading options for the primary effect under consideration. Thus, nearby high-intensity nodes will not necessarily benefit a small neighbor.

## III Results

### III.1 Land values are power law distributed

From empirical data, land values in Sweden are demonstrated to follow a power-law distribution for the higher range of land values (see figure 2). The sharp transition that can be observed around $`60`$kSEK suggests that two truly different mechanisms generate the prices below and above this point. This is in agreement with the observation that the pricing mechanism we suggest would apply only to trading activity.

The data we have used is based on the land value component of market value estimations of about 2.9 million units of real estate in Sweden. The data was originally compiled by the Swedish National Land Survey and coded by Sweden Statistics to geographical coordinates. The data points we use are aggregated land values into $`100m\times 100m`$ squares.

Our results are supported by a recent study by Kaizoji that shows scale-free behavior of land prices in Japan, with an exponent $`\gamma `$ ranging from $`2.53`$ to $`2.76`$Kaizoji\_2003 .

### III.2 Power law distributed prices are predicted by the non-spatial model

We have developed a simple network model of the urban economy based on the Barabási-Albert model by mapping fundamental assumptions from urban economics to the ontology of the network model (see section II.1). In our derivation, we have determined the exponent of the node degree distribution, and thus, per our definition, also the predicted land values, to follow equation (6). To the extent that our interpretation of the underlaying dynamics is correct this demonstrates why urban land values can be expected to follow a power law and how the exponent may depend on parameters.

### III.3 The spatial model retains power law statistics

As discussed in section II.2.5, the impact of spatial constraints closely resembles that of node fitnessBianconi\_Barabasi\_2001a ; Ergun\_Rodgers\_2001 . This could potentially result in the distribution for the spatial network model becoming a sum of many power laws with different exponents.

In figure 2 we plot results from simulations showing that node degrees in the spatial model follow a power law distribution. The model parameters have been tuned (see sec. II.2.5) to reproduce the distribution of the observed land prices.

## IV Conclusion

We present empirical evidence that land prices follow a power law distribution for urban land uses, and present a generic model, based on the underlaying trade network, that reproduces this behavior. The model is a version of the Barabási-Albert scale free network that is also extended to incorporate spatial constraints.

The applicability of the network paradigm to urban growth suggests that scale invariance in urban systems is caused by: i) growth and ii) preferentiality in how new trade connections are formed between areas of land. Growth in this context refers to the continual development of new land. Preferentiality is a consequence of point-to-point interactions between activities occupying the land areas. Note that many other observables, such as population and urban land use intensity, might be highly correlated with land value.

The spatial model that we present can have more general applicability beyond urban economics. Other spatially growing networks are communication networks, transportation networks, electricity and utility networks. These can be expected to follow a similar type of growth since they are intimately connected to urban activity.

The network architecture is generic and allows for addition of any amount of detail. Also, being based on trade relations, the model produces output in units of currency. Because of this, such network models can provide a bridge between a microscopic dynamics that can be found empirically and emergent economic properties.

Further possible directions for research on urban economic networks include interpretation of other theoretical network results in terms of urban dynamics, finding empirical parameters for scenario predictions and model validation.

## V Acknowledgements

We would like to thank Martin Nilsson and Steen Rasmussen for valuable input and discussions.


text:
# Spin dependent quantum transport effects in Cu nanowires

## References