text:
# On the formation of an eccentric disk via disruption of a bulge core near a massive black hole

## 1 Introduction

Ever since the double-peaked nucleus of M31 was clearly resolved by Hubble Space Telescope (HST) imaging (Lauer et al., 1993), the morphology of this system has been a challenge to explain. The most successful kinematic model is the eccentric stellar disk proposed by Tremaine (1995) where stars in apsidally aligned elliptical orbits about a massive black hole result in an over-density at their apoapses, thus accounting for the brighter peak P1. The black hole itself resides near the fainter, more centrally located peak denoted P2. Though this model has been successful at matching the observed velocity and luminosity distribution (Kormendy & Bender 1999; Statler et al. 1999; Bacon et al. 2001; Sambhus & Sridhar 2002), eccentric disk formation remains a mystery.

Touma (2002) showed that a small fraction of counter rotating stars, possibly originating from an disrupted globular cluster on a retrograde orbit, could cause a more massive pre-existing stellar disk to develop a lopsided or $`m=1`$ instability. The self-consistent kinematic modeling of Sambhus & Sridhar (2002), which requires a small percentage of counter rotating stars, supports this scenario. Bacon et al. (2001) proposed that a collision from a passing molecular cloud or globular cluster could knock a pre-existing stellar disk off-center, resulting in a long lived, precessing mode. The N-body simulations of Jacobs & Sellwood (1999); Taga (2002) have illustrated that lopsided stellar disks could be long-lived.

The scenarios discussed by Sambhus & Sridhar (2002); Touma (2002); Jacobs & Sellwood (1999); Taga (2002); Bacon et al. (2001), begin with an initially axisymmetric stellar disk, which then becomes lopsided either because of a violent event (such as a collision with a globular cluster) or the growth of an instability. However Bekki (2000a) proposed that M31’s eccentric disk could have been a result of a single disruption event. He succeeded at producing an eccentric stellar disk by disrupting a globular cluster near a massive black hole in an N-body simulation, though an SPH simulation of the disruption of a massive gaseous cloud did not yield an eccentric disk (Bekki, 2000b). The massive black hole in in M31 is moderate with a mass $`37\times 10^7M_{}`$ (Kormendy & Bender, 1999; Bacon et al., 2001). The eccentric disk itself is nearly as massive as the black hole, $`3\times 10^7M_{}`$ (Peng, 2002; Sambhus & Sridhar, 2002; Bacon et al., 2001), presenting a problem for the globular cluster disruption scenario proposed by Bekki (2000a) which assumed that the globular cluster was typical of Galactic globular clusters and of order a million solar masses.

Here we reconsider Bekki’s proposal, that a single disruption event could have produced the eccentric disk in M31. The disruption event is likely to be complex so we focus on what final disks are capable of supporting lopsided a slow modes, drawing from the recent work of Tremaine (2001) who developed a formalism for the purpose of predicting the precession rates and eccentricities of discrete modes for low mass, nearly Keplerian disks.

We place limits on the density and core radius of progenitor clusters for eccentric disks near massive black holes. The resulting diagram is useful toward predicting what types of galaxies are most likely to harbor eccentric disks, for predicting the probability that eccentric disk formation events occur and for guiding the initial conditions of N-body simulations which may determine if and how they form.

## 2  Tidal disruption of a cluster

The observed correlation between mass of a black hole and the bulge dispersion (Gebhardt et al., 2000; Ferrarese & Merritt, 2001) allows us to relate the black hole to the bulge of its host galaxy.

$$\mathrm{log}(M_{bh}/M_{})=(8.13\pm 0.06)+(4.02\pm 0.32)\mathrm{log}\left(\frac{\sigma _{}}{200\mathrm{k}\mathrm{m}\mathrm{s}^1}\right)$$
(1)
where $`\sigma _{}`$ is the stellar bulge dispersion, and we adopt constants in the relation given by Tremaine et al. (2002). For a discussion on the systematics of measuring the velocity dispersion within $`r_e`$, the effective or half light radius, see Tremaine et al. (2002).

The transition radius (sometimes called the sphere of influence) where the gravity from the black hole takes over from the bulge is

$$r_t=\frac{GM_{bh}}{2\sigma _{}^2}=0.5\mathrm{pc}\left(\frac{M_{bh}}{10^7M_{}}\right)\left(\frac{200\mathrm{k}\mathrm{m}/\mathrm{s}}{\sigma _{}}\right)^2$$
(2)
Using the above relation between the bulge and the black hole mass,

$$r_t1.8\mathrm{pc}\left(\frac{M_{bh}}{10^7M_{}}\right)^{0.5}.$$
(3)
We now consider disruption of a stellar cluster or molecular cloud by a massive black hole. The tidal force from the black hole and bulge is

$$F_{tidal}=\left[\frac{2GM_{bh}}{r^3}+\frac{2\sigma _{}^2}{r^2}\right]s$$
(4)
at distance $`r`$ from the galaxy nucleus and distance $`s`$ from the cluster nucleus. We set this equal to the gravitational force from the cluster on itself to determine how the cluster disrupts.

We describe the nucleus of a cluster with two parameters, the core radius, $`r_0`$, and the density within this radius, $`\rho _0`$. The cores of isothermal spheres and King models are characterized with with these two parameters (e.g., Binney & Tremaine 1987). Within the core radius of the cluster, we assume that the cluster has nearly constant density, so its gravitational potential is proportional to $`s`$. We set the tidal force from the galaxy equal to that from the cluster and determine the distance $`r_d`$, (the disruption radius) when the two are equal. The distance from the galaxy nucleus where the entire core of the cluster disrupts is then set only by the cluster core density, $`\rho _0`$.

$$\frac{GM_{bh}}{r^3}+\frac{\sigma _{}^2}{r^2}G\rho _0$$
(5)
Inside the transition radius the disruption radius

$$r_d\left(\frac{M_{bh}}{\rho _0}\right)^{1/3}=4.6\mathrm{pc}\left(\frac{M_{bh}}{10^7M_{}}\right)^{1/3}\left(\frac{10^5M_{}\mathrm{pc}^3}{\rho _0}\right)^{1/3}.$$
(6)
Outside the transition radius

$$r_d\left(\frac{\sigma _{}^2}{G\rho _0}\right)^{1/2}=4.7\mathrm{pc}\left(\frac{\sigma _{}}{100\mathrm{k}\mathrm{m}/\mathrm{s}}\right)\left(\frac{10^5M_{}\mathrm{pc}^3}{\rho _0}\right)^{1/2}.$$
(7)
We expect that it is difficult to make a slowly precessing eccentric disk in a region where the black hole does not dominate the gravitational field. Outside the sphere of influence, the potential is strongly non-Keplerian. The requirement that the cluster be disrupted (survive to) within the transition radius is $`r_d<r_t`$. Using the previous equation, the definition for $`r_t`$, and the relationship between the bulge stellar dispersion and the black hole mass this constraint gives the requirement

$$\rho _07.4\times 10^5M_{}\mathrm{pc}^3\left(\frac{M_{bh}}{10^7M_{}}\right)^{0.5}.$$
(8)
From this we see that the less massive the black hole, the more dense a cluster would be required to form an eccentric disk. This condition is our first constraint on progenitor clusters for eccentric disks.

At radii larger than the disruption radius, the outer parts of the stellar cluster or galaxy core would be stripped. We expect that the radial distribution of the resulting disk will depend upon the concentration of the cluster or radial profile of the infalling galaxy bulge or core. Since the densest part of the cluster is disrupted at $`r_d`$, we expect that this part is most likely to be involved in any large eccentricity amplitude variation. So we can estimate the mass of the active disk formed at disruption using the core radius of the cluster, $`r_0`$, resulting in an active disk of mass equal to the mass within the core of the cluster, $`M_c\rho _0r_0^3`$. The width of this active disk would be of the order of the core radius of the cluster.

Tremaine (2001); Lee & Goodman (1999) assume that the mass of the disk is less than that of the black hole. If we consider the active part of the disk, this is even true in M31 (the center of mass is closest the black hole). It’s difficult to imagine a stable situation when the lopsided part of the disk is as massive as the black hole so we arrive at another constraint on our progenitor cluster,

$$M_c<M_{bh}.$$
(9)
It is unlikely that a thin disk will result following disruption unless the disk Mach number, $`(r)>1`$ (Tremaine, 2001),

$$(r_d)\frac{\mathrm{\Omega }(r_d)r_d}{\sigma _c}>1$$
(10)
where $`\mathrm{\Omega }(r_d)`$ is the angular rotation rate at the disruption radius $`r_d`$, and $`\sigma _c`$ is the central cluster dispersion.

Following disruption of the cluster, the energy spread of the particles is set by the cluster core radius, $`r_0`$, so we expect that the resulting width of the disk $`r_0`$. The angular momentum spread at disruption is also set by the size of the cluster core, $`\frac{\mathrm{\Delta }J}{J_0}\frac{r_0}{r_d}`$ where $`J_0`$ is the angular momentum of the bulk motion of the cluster at disruption. Consequently following disruption we expect a disk radial velocity dispersion at least

$$\sigma _c\sqrt{\frac{GM_{bh}}{r_d}}\frac{r_0}{r_d}.$$
(11)
This scaling is similar to that outlined by Johnston et al. (2001) for the disruption of dwarf galaxies. Alternately we can use the relation between central dispersion, core density and radius for a King model ($`\sigma _c^2=4\pi G\rho _or_0^2/9`$). Using either estimate for $`\sigma _c`$ and inserting into equation(10), assume that we are within the transition radius so that $`\mathrm{\Omega }=\sqrt{GM_{bh}/r^3}`$ we estimate

$$(r_d)\frac{M_{bh}^{2/3}}{\rho _0^{2/3}r_0^2}1$$
(12)
which is equivalent to requiring $`M_{bh}>M_d`$.

In summary we can constrain properties of an incoming cluster which subsequently forms an eccentric disk:

1. $`r_d<r_t`$. The disruption radius should be less than the transition radius.
2. $`r_0<r_t`$. The cluster core radius should be less than the transition radius.
3. $`M_c<M_{bh}`$. The cluster core mass should be less than the black hole mass. This constraint is equivalent to requiring $`r_0<r_d`$ or $`>1`$.

We plot these three constraints as lines in Figure 1 for three different black hole masses, and as a function of cluster core radius, $`r_0`$, and density, $`\rho _0`$. The allowed region of parameter space for a disrupted cluster to form an eccentric disk is the upper left hand corner, corresponding to core radii above about a pc, and core densities above about $`10^5M_{}\mathrm{pc}^3`$.

Also shown in Figure 1 are measured core radii and central densities (assuming a mass to light ratio of $`M/L_V=5`$ in solar units) for the Milky Way globular clusters based on the quantities compiled by Harris (1996).111We used the updated catalog available at http://physun.physics.mcmaster.ca/ harris/mwgc.dat

We also include on Figure 1, the massive globular cluster G1 (based on core properties estimated by Meylan et al. 2001) and the cores of M31, M32 and M33, based on density profiles (or limits in the case of M32) estimated by Lauer et al. (1998). Limits on the density and core radius are placed on the figure 1 for the two densest galaxies (other than M31, M32, and M33) studied by Faber et al. (1997). Higher angular resolution observations would be needed to find out if more distant galaxies can have cores with the high stellar densities of M31, M32 and M33.

We see from Figure 1 that Galactic globular clusters are capable of forming eccentric disks if they are extremely dense and have compact cores, and they are more likely to do so in galaxy nuclei containing larger black holes. However, to provide a plausible progenitor for an eccentric disk, a cluster must not only have a dense and compact core but be massive enough to account for the mass of the entire disk. Galactic globular clusters are not massive enough to account for the disk in M31.

### 2.1 Bekki’s simulations

We examine the properties of the simulated disruptions done by Bekki (2000a, b) to see if they are consistent with our proposed progenitor cluster limits. In an N-body simulation, Bekki (2000a) succeeded in producing a stellar eccentric disk by disrupting a globular cluster of mass 0.1 times the black hole mass, for $`M_{bh}=10^7M_{}`$. The bulge he adopted was described by an NFW profile with a scale radius of $`100`$pc. For this adopted bulge profile we calculate that the mass enclosed within 10pc was less than a million $`M_{}`$ so that the entire disruption took place within what we have called the transition radius, where the gravity from the black hole takes over from the bulge (equation 2). His simulated cluster disrupted at a radius of $`r_d3`$pc from the black hole, and since the cluster core was much smaller than this, $`r_00.4`$pc, the stellar system was capable of forming an eccentric disk. We conclude that the success of his simulation in producing the eccentric disk is consistent with the crude limits we have placed on progenitor cluster properties.

The SPH simulation by Bekki (2000b), however, produced a circular disk, not an eccentric one. The gas cloud in this simulation was of similar mass to that of the stellar cluster, (Bekki, 2000a) but had a significantly larger scale length of 10pc. The cluster disrupted at $`r_d10`$pc which is approximately the same as the core radius of the gas cloud. We suspect that a gas cloud with a smaller core radius would have been capable of producing an eccentric disk. Again, we find that the failure of this simulation in producing the eccentric disk is consistent with the crude limits we have placed on progenitor cluster properties.

### 2.2 M31’s eccentric disk and possible progenitor cluster or galaxy core properties

The two peaks observed in the nucleus of M31 are separated by 1.87 pc (Lauer et al., 1993) and the black hole mass is $`3\times 10^7M_{}`$ (Kormendy & Bender, 1999), so we expect a transition radius of $`r_t3.1`$ pc which is, as we expect, outside the location of the eccentric part of the disk. The width of the eccentric part of the disk is about 0.2” or $`0.7`$pc so we would estimate that if the disk formed from a disrupted cluster that its core radius would have been smaller. The total mass in the most eccentric part of the disk is about $`10^6M_{}`$ (Bacon et al., 2001) so we would require a progenitor cluster of $`r_00.7`$pc and $`\rho _0`$ a few times $`10^6M_{}\mathrm{pc}^3`$. For $`\rho _04\times 10^6M_{}\mathrm{pc}^3`$ the disruption radius $`r_d`$ (estimated by equation 6) is consistent with the location of the eccentric disk itself.

To provide a plausible progenitor, a cluster must not only have a dense and compact core but be massive enough to account for the mass of the entire disk. The extended, nearly circular part of the M31 stellar disk dominates its mass and has been estimated to be $`3\times 10^7M_{}`$ (Peng, 2002; Sambhus & Sridhar, 2002; Bacon et al., 2001). This is massive compared to Galactic globular clusters, presenting a problem for the globular cluster disruption scenario proposed by Bekki (2000a). However extremely massive extra-galactic clusters have been identified. For comparison, we show in Figure 1 the location of the massive globular cluster (or perhaps core of a dwarf galaxy) which is 40kpc from the nucleus of M31. This cluster has a mass of $`717\times 10^6M_{}`$ (Meylan et al., 2001), a core radius $`r_00.14\mathrm{"}=0.52`$pc, and a core density $`\rho _04.7\times 10^5M_{}\mathrm{pc}^3`$ (Meylan et al., 2001). G1 is massive enough that it could have accounted for the mass in the extended (and almost circular) part of M31’s eccentric disk, though its tidal radius is large, $`200`$pc, which implies that much of the cluster would have disrupted well outside the location of the eccentric disk in M31.

The estimated progenitor core properties required to form the M31 eccentric disk are surprising if we compare them to Milky Way globular clusters and nearly matched by the exotic and extremely massive extra-galactic cluster G1. However, they are similar to what has been measured in galaxy nearby cores, such as M33, M32 and M31 itself, which we have also placed on Figure 1. We see from this figure that galaxy cores can have the high densities and compactness to be eccentric disk progenitor candidates. Moreover, they are usually more massive than globular clusters and so can more easily account for the mass in M31’s eccentric disk.

### 2.3 NGC 4486’s eccentric disk and its progenitor

The other well known example with a double nucleus is NGC 4486B, with a separation of 0.15” or 12 pc between the two isophotal peaks Lauer et al. (1996) and a black hole mass of $`M_{bh}10^9M_{}`$ (Magorrian et al., 1998). The transition radius in this galaxy would be $`r_t17`$pc; so again we find that the eccentric stellar disk lies within $`r_t`$. From equation(6) and using 12 pc for the disruption radius we estimate that the progenitor core density would be $`\rho _06\times 10^5M_{}\mathrm{pc}^3`$. The width of the peaks comprising the disk (and hence the progenitor core radius) is estimated to be 2–4 pc. However, NGC 4486A is much more distant than M31, and consequently the eccentric disk mass is probably much higher than that of M31, probably exceeding $`10^8M_{}`$ and this is too massive for a globular cluster. However as in the case of M31’s estimated progenitor, the central densities, core radius and total mass are reasonable for a galaxy core.

### 2.4 Comparison to the properties of galaxy cores

In equation(2) we have assumed an isothermal profile for the stellar density profile of the background galaxy. However, as shown by multiple studies (e.g., Faber et al. 1997; Lauer et al. 1998) galaxy cores at small radii are seldom fit with a surface brightness profile proportional to $`R^1`$ (where $`R`$ is the observed angular radius from the nucleus on the sky) corresponding to an isothermal density profile $`r^2`$. Images from the HST have been part of a major effort to classify the nuclear stellar profiles in early-type galaxies, resulting in the classification of light profiles into two categories, galaxies with shallow inner cusps, denoted ‘core-type’ profiles and galaxies with ‘power-law’ light profiles (Lauer et al., 1995; Faber et al., 1997). Power-law type galaxy cores tend to have the steepest nuclear surface brightness profiles $`\mu (R)R^\gamma `$ with $`\gamma `$ nearly equal to 1 in only the the most extreme cases. Most of the galaxies studied by Faber et al. (1997) had shallower profiles. This implies that equation(3) somewhat underestimates the transition radius in most galaxies, particularly in the galaxies classified as core-type in which $`\gamma `$ is close to zero. We can regard our estimated radius (equation 3) as a lower limit on the transition radius for all but the galaxies with the steepest nuclear profiles.

We now consider the problem of two merging galaxies, both with more complex stellar density profiles and both with massive black holes. We can approximately describe the primary as having a transition radius given by equation(3) (which is actually a lower limit as mentioned above). If the secondary has a power-law form for its density profile then it will not completely disrupt at a particular radius, like the non-singular isothermal sphere or King models.

One component of the Nuker surface brightness profile (Faber et al., 1997) is the break radius, $`r_b`$, where the slope of the surface brightness profile changes. For profiles denoted core-type, within $`r_b`$, the surface brightness (and so density) rises much less steeply with decreasing radius than outside it. Because the surface brightness for core-type galaxies only increases slightly within $`r_b`$, we can associate the break radius and density estimated at this radius, with the core radius $`r_0`$ and density $`\rho _0`$ which we have used in the previous sections. This lets us determine if the bulge or stellar component of the secondary galaxy will survive tidal truncation to within the transition radius of the primary. A secondary with break radius exceeding its own transition radius $`r_b>r_t`$ (likely to be true for luminous galaxies) will almost completely disrupt when the tidal truncation radius reaches its break radius. The black hole of the secondary will then spiral in to the nucleus of the primary, leaving the stellar component of the secondary behind.

However for bulges and fainter ellipticals, we expect $`r_b<r_t`$. For example, M32’s profile is steeply rising well within its transition radius $`r_t0.5`$pc and increases in density to its last measured point (Lauer et al., 1998). In this case the core would be continuously disrupted, though we expect a change in the surface brightness profile of the resulting disk set by the location at which the tidal truncation radius of the secondary was equal to its own transition radius. The continuous stripping of the secondary in the presence of the secondary’s black hole may prevent the formation of a lopsided disk.

### 2.5 Tidal disruption of a power-law galaxy

The radius of tidal disruption can be estimated by comparing the mean density of the object to that of the object disrupting it (e.g., Sridhar & Tremaine 1992). At a distance $`r`$ from the nucleus of the primary galaxy, the mean density resulting from the black hole is $`\overline{\rho }(r)=M_{bh,1}/(\frac{4\pi r^3}{3})`$ where $`M_{bh,1}`$ is the mass of the primary’s black hole. We have assumed that $`r<r_{t,1}`$ where $`r_{t,1}`$ is the transition radius of the primary, so the contribution from the primary’s bulge is negligible.

For a secondary galaxy with a power-law law density profile (see Binney & Tremaine 1987 section 2.1e)

$$\rho (s)=\rho _a\left(\frac{a}{s}\right)^\alpha $$
(13)
the mass integrated out to radius $`s`$ is

$$M(s)=\frac{4\pi \rho _aa^\alpha }{(3\alpha )}r^{3\alpha }$$
(14)
so the mean density within this radius

$$\overline{\rho }(s)=\frac{3\rho _a}{(3\alpha )}\left(\frac{a}{s}\right)^\alpha .$$
(15)
for radii $`s>r_{t,2}`$ outside the transition radius of the secondary. We set the mean density of the primary within $`r`$ (due to the primary’s massive black hole) equal to the that of secondary (within stripping radius $`s`$), $`\overline{\rho }(s)=\overline{\rho }(r)`$, and find that

$$r^3\frac{M_{bh}(3\alpha )}{4\pi \rho _a}\left(\frac{s}{a}\right)^\alpha .$$
(16)
When the secondary is at a distance $`r`$ from the nucleus of the primary, it will be stripped out to a distance $`s`$ from its own nucleus, where $`s`$ is related to $`r`$ by the previous equation.

The stripping radius, $`s`$, is the same size as the distance from the primary’s nucleus, $`r`$, at a radius

$$r_e=\left[\frac{M_{bh,1}(3\alpha )}{4\pi \rho _aa^3}\right]^{\frac{1}{3\alpha }}a$$
(17)
where $`M_{bh,1}`$ is the mass of the primary’s black hole. When $`s`$ estimated from equation(16) exceeds $`r`$ then the secondary engulfs the nucleus of the primary.

For a disrupted secondary galaxy core to produce an eccentric disk, we expect that the radius at which the secondary engulfs the nucleus of the primary, $`r_e`$, should be smaller than the transition radius of the primary ($`r_{t,1}`$). In this case, the core of the secondary will survive intact within the transition of the secondary. A reasonable condition for formation of an eccentric disk via the disruption of a galaxy with a power-law density profile should be

$$r_e<r_{t,1}.$$
(18)
Using equations(17,3) this is equivalent to

$$\rho _a>M_{bh,1}^{\frac{\alpha 1}{2}}a^\alpha \kappa ^{\alpha 3}\frac{(3\alpha )}{4\pi }$$
(19)
where $`\kappa =5.7\times 10^4\mathrm{pc}M_{}^{1/2}`$. For an isothermal profile, the limit on the density at $`a=1`$pc

$$\rho _{a=1\mathrm{p}\mathrm{c}}>5\times 10^6M_{}\mathrm{pc}^3\left(\frac{M_{bh,1}}{10^7M_{}}\right)^{1/2}\mathrm{for}\alpha =2.$$
(20)
### 2.6 Existing simulations of galaxy mergers

With the exception of Merritt & Cruz (2002); Holley-Bockelmann & Richstone (2000), few simulations have been carried out with realistic stellar density profiles and central black holes. Merritt & Cruz (2002) presented N-body simulations of two merging galaxies, the primary with a shallow central density profile and the secondary with a steeper profile at the nucleus. They did simulations with and without massive nuclear black holes to illustrate the difference in remnant profiles following the merger. The density profiles adopted for their simulations were generated from the family introduced by Dehnen (1993).

To explore the role of the black holes on the disruption of the secondary in their simulations we estimate the transition radii for each galaxy. At small radii, $`r<a`$, both galaxies were power-law in form

$$\rho (r)=\frac{(3\alpha )M}{4\pi a^3}\left(\frac{a}{r}\right)^\alpha $$
(21)
where $`M`$ was the mass of the galaxy and $`a`$ was its scale length. Merritt & Cruz (2002) used $`\alpha _1=1`$ for the primary, and $`\alpha _2=2`$ for the secondary (subscripts here refer to the primary and secondary galaxy respectively). We estimate the transition radius for each galaxy (sphere of influence) by calculating the radius at which the integrated stellar mass within equals that of the central black hole. For the primary galaxy

$$r_{t,1}a_1\left(\frac{M_{bh,1}}{M_1}\right)^{1/2}=0.035a_1$$
(22)
where we have integrated equation(21) and used the ratio $`M_{bh}/M=0.0012`$ (that adopted by Merritt & Cruz 2002). For the secondary galaxy we estimate

$$r_{t,2}a_2\left(\frac{M_{bh,1}}{M_1}\right)=0.0012a_2.$$
(23)
Because the transition radius of the secondary was so small, its black hole should not have affected the simulations (they were probably not well resolved on that scale).

Both simulations produced similar remnants at radii larger than 0.1 times the primary’s scale length. The change in the surface brightness profiles at this radius was probably due to the break radius (scale length at which the density profile changes slope) of the secondary.

For the simulations containing both massive black holes, a change in radial density profile is seen in the remnant at approximately the primary’s transition radius (see their figure 2d.) Using their profiles, we estimate that $`r_e0.01`$ times the scale length of the secondary, which is consistent with the survival of a fraction of the secondary’s core to well within the transition radius of the primary. The constraints developed in the previous section were not violated ($`r_e<r_{t,1}`$), so this type of simulation might have been capable of forming an eccentric disk in the remnant. If the simulation did not form one, then simultaneous disruption of the densest part of the cluster might be required to form an eccentric disk. This would occur when the core has a well defined core radius (and nearly constant density within), and may not occur when the secondary has a steeply increasing density profile all the way to its nucleus.

In the simulations lacking the massive black holes (see their figure 2h), the secondary did not disrupt within its break or scale radius. This follows because the mean density of the secondary significantly exceeds the mean density of the primary within their scale radii, consequently the primary could not have disrupted the secondary within the secondary’s scale radius.

The merger simulations of Holley-Bockelmann & Richstone (2000) were done with a quite massive primary black hole $`M_{bh}>10^9M_{}`$ and lacked a secondary black hole. From their initial density profiles we estimate that the size of the secondary at its disruption radius in their 10:1 simulation was only about a third the disruption radius itself. This implies that a disk could be formed from the secondary, following disruption. Holley-Bockelmann & Richstone (2000) reported that their simulations produced a spinning remnant disk, however it was was thick, consistent with our crude estimate of the secondary’s largish size during disruption. Much of the secondary survived to within the transition radius of the pimary’s black hole (about 30 pc). It is possible that a denser secondary (better satisfying the conditions in eqns 18–20) could result in a thinner stellar disk which might then be more likely to be lopsided. Alternatively, a secondary with a shallow core (that would would disrupt all at once, as simulated by Bekki 2000a) might be required for eccentric disk formation.

## 3 Discussion and Summary

In this paper we have considered the possibility that a single disruption event could result in the formation of an eccentric stellar disk such as are found in M31 and possibly NGC 4486B. We explore the scenario of a galaxy core which is stripped as it spirals in toward a nuclear massive black hole until it reaches a critical radius, denoted here as $`r_d`$, the disruption radius. In the simplest model of an isothermal sphere, the disruption radius is set by the cluster core or King radius and central density. We suggest that an eccentric disk cannot form unless the core of the galaxy or stellar cluster is disrupted within the sphere of influence of the massive black hole, referred to here as the transition radius, $`r_t`$, and the core radius is smaller than the radius at which it disrupts. To succeed in making an eccentric disk, we find that the progenitor cluster core must be denser and more compact for less massive black holes.

For the eccentric disk in M31, the progenitor cluster core density must be greater than $`10^5M_{}\mathrm{pc}^3`$ and core radius smaller than a pc. From the radius, width and mass of the most eccentric region of the disk we estimate that the density of the core was a few times $`10^6M_{}\mathrm{pc}^3`$ and the core radius was of order a half pc.

Massive extragalactic globular clusters such as G1 might be dense, compact and massive enough to be an eccentric disk progenitor, however they are probably too diffuse to account for the large disk mass in M31 within a few parsecs of its black hole.

Dense bulges such as found in M31 itself, and in lower luminosity ellipticals galaxies such as M32, are dense and compact enough that they could be similar to progenitors for M31’s and NGC 4486B’s eccentric disks. Nuclear star clusters, such as found in M33, may also be dense and compact enough to be progenitors. Though the mass of M33’s cluster (a few million $`M_{}`$) is too small to be a progenitor for M31’s eccentric disk, the nuclear stellar clusters observed by (Böker et al., 2002) in late-type galaxies have half light radii of order 5pc, range in their estimated masses from $`10^610^8M_{}`$, and so could be progenitors for eccentric disks following disruption by a massive black hole. Because of the large mass of bulges and nuclear star clusters at small radii, it would be easier to account for the high masses of the two known eccentric disks by disrupting them, than possible by disrupting a globular cluster. While the massive black hole in M32 might prevent a progenitor similar to M32 from forming an eccentric disk, nuclear clusters such as found in M33 can lack massive black holes and so might provide better progenitor candidates. To date the high angular resolution of HST has resolved extremely high stellar densities or order $`10^6M_{}\mathrm{pc}^3`$ in only the nearest galaxy bulges. Until higher angular resolution observations are available, we will not know if such high density galaxy cores are common.

The two galaxies with double nuclei, M31 and NGC 4486B, exhibit only moderate color variations in their nuclei (Lauer et al., 1996, 1998), a situation that could be a natural consequence of a scenario that involves the merging of galaxy bulges (proposed here), and is more difficult to explain with a scenario that forms a younger disk in situ (an ingredient of the formation scenarios discussed by Bacon et al. 2001; Touma 2002; Taga 2002; Jacobs & Sellwood 1999).

The scenario proposed here is based on a simple tidal disruption argument and can be most quickly tested with N-body simulations such as have been carried out by Merritt & Cruz (2002); Bekki (2000a); Holley-Bockelmann & Richstone (2000). The simulations of Bekki (2000a) established that the disruption of a cluster could result in the formation of an eccentric disk, and Merritt & Cruz (2002); Holley-Bockelmann & Richstone (2000) have carried out simulations based on realistic galaxy profiles and with massive black holes. Holley-Bockelmann & Richstone (2000) established that a spinning nuclear stellar disk can be a remnant. It remains to be seen whether the merger of two galaxy cores can result in the creation of an eccentric stellar disk.

We suspect that the merger of a primary galaxy containing a very massive black hole with a secondary with a ‘core-type’ or shallow central surface brightness profile, and significantly lower mass black hole, would be most likely to form an eccentric disk with N-body simulations. For the core-type galaxies, the break radius and density at this radius provide an equivalent for the King or core radius and central density used in Figure 1 and so can be used to estimate the likelihood that the core can disrupt to form an eccentric disk.

If the merger of galaxy bulges can result in the formation of an eccentric disk, then their formation would be a natural consequence of hierarchical galaxy formation, and can also be used to probe the properties of the parents of galaxies which contain them.

This work was initiated by discussions with Joel Green and Rob Gutermuth in the class Astronomy 552 at the University of Rochester during the fall of 2002. We thank Chien Peng, Ari Laor and Eric Emsellem for helpful discussions and correspondence.


format: text

text:
# Ferromagnetic resonant tunneling diodes as spin polarimeters and polarizers

## Abstract

A method for measuring the degree of spin polarization of magnetic materials based on *spin*-dependent resonant tunneling is proposed. The device we consider is a ballistic double-barrier resonant structure consisting of a ferromagnetic layer embedded between two insulating barriers. A simple procedure, based on a detailed analysis of the differential conductance, allows to accurately determine the polarization of the ferromagnet. The spin-filtering character of such a system is furthermore addressed. We show that a $`100\%`$ spin selectivity can be achieved under appropriate conditions. This approach is believed to be well suited for the investigation of diluted magnetic semiconductor heterostructures.

The operation of spintronic devices Wolf  requires the availability of efficient techniques to inject and detect spin-polarized currents. In fact the search for new materials with a high degree of spin-polarization represents today an important technological challenge in spintronics. A key parameter in these devices is the spin polarization of the electric current defined as

$$𝒫=(I_{}I_{})/(I_{}+I_{}),$$
(1)
where $`I_\sigma `$ is the contribution to the current due to $`\sigma `$-spin carriers Mazin . Structures containing ferromagnetic materials are obvious candidates as sources of spin-polarized currents. These can be metallic yuasa  or consist of diluted magnetic semiconductors (DMSs) Ohno . Two different electric transport methods have been devised for determining this quantity both based on contacting the ferromagnets under investigation to a superconductor (S). In the first one, $`𝒫`$ can be estimated by realizing a F/S tunnel junction and Zeeman-splitting the superconducting density of states (DOS) through the application of an external magnetic field TM1 . More recently the point-contact Andreev reflection (PCAR) method was introduced. Here the suppression of Andreev reflection (AR) at F/S ballistic junctions is exploited Ji . This method proves most effective for highly-transmissive F/S contacts where the AR probability is high parker , consequently its application in DMS materials is hindered by the presence of an unavoidable Schottky barrier at the interface.

In this letter we propose an alternative route to determine the polarization of a ferromagnet (either metallic or DMS) based on a double-barrier resonant structure (see Fig. 1(a)). The operation principle of this device resides in the fact that, when the resonator is ferromagnetic, the resonant levels of the two spin species occur at different energies whose difference is a function of the exchange field in F (and therefore $`𝒫`$). The advantage of this method stems from the fact that this energy difference is little dependent on the quality of the barriers (provided their transmissivity is low enough as required for resonant tunneling to occur) and that no direct coupling to a superconductor is necessary, as imposed in PCAR. In addition the same device can be operated as an extremely efficient spin filter Ralph ; rocca .

The device consists of a ballistic ferromagnetic layer of thickness $`L`$ sandwiched between two non-magnetic barriers connected to electrodes (in Fig. 1(b) the potential profile of the structure is sketched). In the following we refer to this structure as to the ferromagnetic resonant tunneling diode (FRTD). As long as phase-coherence in preserved, current through the FRTD can flow via the resonant levels of the quantum well defined by the barriers. Such resonances show up as peaks in the differential conductance spectrum. Since the resonator is ferromagnetic, the spin degeneracy is lifted and the resonant levels relative to the two electron spin species occur at different energies (see Fig. 1(c)). By noting that this energy difference depends on the polarization $`𝒫`$ of the ferromagnet, it is possible to determine $`𝒫`$ from the analysis of the resonance energy position.

In order to clearly explain how the method works, we consider a single-band model as a generic example, and we describe the ferromagnet (or DMS) with an effective exchange field $`h_{exc}`$ (Stoner model) koenig . If the resonator is decoupled from the contacts, the quantum well levels

$$E_\sigma ^{n_\sigma }=(k_{n_\sigma })\sigma h_{exc},$$
(2)
(characterized by the spin-dependent index $`n_\sigma `$) are found from elementary quantum mechanics. In Eq. (2), $`k_{n_\sigma }`$ is the quantized longitudinal quasiparticle wavevector, $`(k)`$ is a generic dispersion relation, and $`\sigma =\pm 1`$ is the spin index, where the sign $`+`$($``$) corresponds to $``$($``$) spin channel. Note that $`(k)`$ contains a transverse contribution to the total quasiparticle energy. For fixed $`h_{exc}`$ and $`L`$ the requirement that $`E_\sigma ^{n_\sigma }>0`$ (i.e., resonant states above the Fermi energy $`E_F`$) determines the first resonance index $`n_\sigma `$. Due to presence of the exchange field, $`n_\sigma `$ can be in general quite different for the two spin species. By connecting the resonator to the electrodes through weakly-transmitting barriers such bound states acquire a finite width at the energies $`E_\sigma ^{n_\sigma }`$. As a result they manifest themselves as peaks in the differential conductance spectrum with a broadening dependent on barrier strength and temperature.

As in the PCAR method, our aim is to determine $`h_{exc}`$, and therefore $`𝒫`$, assuming that the band structure parameters contained in $`(k)`$ are known. Since the exchange field enters the dispersion relation additively, one might think that it can be determined from the energy difference between a spin-up and a spin-down resonance relative to the same index $`n_\sigma `$. This however is possible only if $`h_{exc}`$ is small with respect to such energy difference, a situation that is typically not fulfilled in ferromagnets and DMS. Nevertheless, $`h_{exc}`$ can still be determined through an analysis of the low-temperature differential conductance spectrum. One can employ the following simple procedure: i) one identifies two successive spin-up resonances and measures their energy difference from which the index $`\overline{n}_{}`$ of the first chosen spin-up resonance is determined solving Eq. (2); ii) the same operation is performed on spin-down resonances in order to obtain $`\overline{n}_{}`$; iii) finally one measures the energy difference $`\delta E_{}`$ between one spin-up resonance of i) and one spin-down resonance of ii) so that the exchange field is easily obtained from

$$h_{exc}=\frac{1}{2}\left[(k_{\overline{n}_{}})(k_{\overline{n}_{}})\delta E_{}\right].$$
(3)
It is worth stressing that for linear or constant dispersion curves the proposed method fails. In these situations indeed the energy difference of points i) and ii) would be independent of the resonance index, thus preventing its determination. To identify the spin character of the resonances we envision the following two possibilities. The first one is to make use of ferromagnetic leads spin , so that the amplitude of the resonances will be spin-dependent (for the same reason for which, for example, the Sharvin conductance of a ferromagnetic point contact of area $`𝒜`$ and Fermi wavevector $`k_\sigma `$, $`G_{Sharvin}=(e^2/h)(𝒜k_\sigma ^2/4\pi )`$, is different for the two spin species). The second possibility relies on the application of an external magnetic field which will produce a Zeeman shift of the resonances in opposite energy direction for the two spin species g .

So far we presented an ideal situation. In order to check the actual feasibility of the procedure outlined above we have calculated, within a transfer-matrix approach, the differential conductance spectrum of a three-dimensional FRTD with rectangular barriers, assuming translational invariance in the plane of the junctions and large spin-flip length in the ferromagnet with respect to the resonator width spinflip . For definiteness, we assumed a parabolic dispersion relation FRT  and typical parameters relative to a metallic ferromagnet calcolo . The resonator width was chosen in order to show a few resonances per spin within a bias voltage range of one tenth of $`E_F`$. Barrier transmissivity ($`10^2`$) was taken such that the resonances width is small with respect to their relative separation, and a finite temperature was also taken into account in order to make the simulation more realistic. Having computed the differential conductance spectrum, as if obtained directly from a measurement, we have identified the peaks positions and spin character. By applying our procedure we have determined $`h_{exc}`$ with a negligible error with respect to the nominal value. This proves the effectiveness of this method. Furthermore, it is remarkable that the method works properly even if the dispersion curve for bias voltages in the range considered is not far from linearity. We wish to remark that all the above simulations are just an exemplification of how the method works and that some experiments might need to be analyzed employing, for example, multi-band models. In this last case the method is still valid even though the analysis becomes more cumbersome.

Finally we discuss the spin-filtering action of the proposed structure. For bias voltages such that only the *first* resonance contributes to transport, i. e., for voltages which correspond to the first step in the current-voltage characteristic Ralph , the system behaves as an ideal spin filter, providing fully-polarized ($`100\%`$) currents even in the absence of a *half-metallic* (i.e., with $`𝒫=1`$) resonator. The spin polarization of the current can be chosen either by varying the diode geometry (the resonator width $`L`$) or by selecting an appropriate material (with a given $`h_{exc}`$). In both situations, the application of an external magnetic field offers an additional way of tuning the spin-filtering character of the structure. To elucidate this property we show in Fig. 2(a) a contour plot of the zero-temperature differential conductance as a function of resonator width and bias voltage. In this plot we assumed that the FRTD electrodes are of the same magnetic material as the resonator. The figure clearly shows that the spin character of the first resonance (indicated by $`,`$ white arrows) can be changed upon variation of $`L`$. For example by choosing $`L=2.2`$ nm the first resonance occurs at $`eV=0.08E_F`$ and is spin-up in character, while choosing $`L=2.6`$ nm a spin-down resonance occurs at $`eV=0.027E_F`$ (see red dashed-lines). A similar dependence is observed in Fig. 2(b) where the differential conductance is plotted versus the exchange field for fixed resonator width ($`L=3`$ nm).

We conclude by elaborating on the actual realization of the FRTD as a spin polarizer. Figure 1(a) shows a FRTD where the magnetization direction of the two ferromagnetic electrodes (F) is pinned through their contact with antiferromagnetic layers (AF), and the magnetization of the resonator can be changed from parallel to antiparallel direction upon the application of weak magnetic fields. The device could be engineered through the sequence of metallic ( for example NixFe1-x alloy) or III-V DMSs (typically Ga1-xMnxAs) stacked layers, and taking advantage of the $`x`$-dependence of their exchange field in order to tune the spin-filter.

In summary, we have proposed a general method which exploits a FRTD for determining the degree of polarization of a ferromagnet. The spin-filtering character of such structure was additionally addressed and ideal spin selectivity proven. The proposed procedure seems to be more appropriate than PCAR for the case of DMSs, where the formation of a transmissive contact with the superconductor is expected to be difficult due to the presence of an unavoidable Schottky barrier at the interface with the metal. On the contrary, the technology for the fabrication of magnetic resonant tunneling diodes, both metallic Ralph ; RT  and consisting of DMSs RT2 , is nowadays already available.

The authors gratefully acknowledge S. Sanvito for useful discussions. This work was supported by INFM under the PAIS project TIN and European Community through grant (RTN2-2001-00440).


text:
# Scaling in Numerical Simulations of Domain Walls

## I Introduction

The idea that the symmetries in nature are not respected by the vacuum plays a crucial role in the unification of forces. Moreover, any broken symmetries that can be identified today were likely to have been restored at high temperatures Kirzhnits:1972iw ; Dolan:1974qd ; Kirzhnits:1976ts . These two facts suggest that the Universe underwent phase transitions in its early history. It was realised that phase transitions may leave the Universe filled with topological defects Zeldovich:1974uw ; Kibble:1976sj , which if massive enough could be observed through their density fluctuations Kibble:1976sj ; Zeldovich:1980gh ; Vilenkin:1983jv  (see also VilShe94 ; Hindmarsh:1995re  for reviews).

Probably the most interesting and important property of defect networks is that they seem to exhibit dynamic scaling. This means that they quickly lose memory of their initial conditions and evolve towards configurations which can be characterised by a single length scale (or perhaps a few Austin:1995dz ) $`\xi `$. This length scale is thought to increases with time with a universal exponent. In a relativistic field theory one can argue from dimensional analysis that $`\xi (t)tf(Mt)`$, where $`M`$ is the mass scale of the defect. The large-scale dynamics of defects are independent of $`M`$, and so the dynamical scaling exponent can be naively estimated as 1.

This behaviour has been checked by numerically simulating classical field theories for domain walls Press:1989yh ; Coulson:1996nv ; Larsson:1997sp , gauge strings Vincent:1998cx ; Moore:1998gp ; Moore:2001px , global strings Ryden:1989vj ; Yamaguchi:2000dy , global monopoles Bennett:1993fy ; Yamaguchi:2001rf ; Yamaguchi:2001xn  and textures Bennett:1993fy ; Durrer:1994da ; Pen:1994nx . All the simulations are consistent with the linear scaling law over the range of the simulations, but do allow other behaviours: in particular, Press, Ryden and Spergel suggested that the results for domain walls would be better fitted by $`\xi t/\mathrm{ln}(t)`$.

Since the original naive scaling arguments were put forward, a more quantitative approach to the dynamics of domain wall networks has been proposed by one of the authors Hindmarsh:1996xv ; HinCop97 ; Hin02 , which predicts not only the linear scaling law for domain walls, but also the amplitude of the relation. Clearly, a logarithmic correction would be a problem for the approach. The main purpose of this paper is to give a more accurate numerical determination of the scaling law (in 2 and 3 dimensions) and to determine the amplitude, to check the accuracy of the theoretical predictions in Hindmarsh:1996xv ; Hin02 .

Our results for the scaling exponents and amplitude can be found in Tables 1, 2. They are consistent with the linear scaling law for $`\xi `$, but the numerically determined amplitudes are higher than the theoretical predictions. A detailed comparison can be found in Hin02 .

## II Dynamics of domain walls

Domain walls occur in field theories whose manifold of minimum energy states is topologically disconnected (see e.g. VilShe94 ). The canonical example in relativistic field theory is a theory of a single scalar field $`\varphi (x)`$, with action

$$S=d^4x\sqrt{g}\left(\frac{1}{2}_\mu \varphi ^\mu \varphi V(\varphi )\right),$$
(1)
where the potential $`V`$ is a function with more than one minimum, which we take to be the renormalisable form

$$V(\varphi )=\frac{\lambda }{4}(\varphi ^2\mu ^2)^2.$$
(2)
In the cosmological context, the metric $`g_{\mu \nu }`$ is taken to have the Friedmann-Robertson-Walker form

$$g_{\mu \nu }=R^2(t)\eta _{\mu \nu },$$

where $`\eta _{\mu \nu }=\mathrm{diag}(1,1,1,1)`$ is the Minkowksi space metric and $`t`$ is conformal time. The field can be conformally rescaled $`R(t)\varphi \varphi `$, giving an Euler-Lagrange equation

$$\ddot{\varphi }+2\frac{\dot{R}}{R}\dot{\varphi }^2\varphi +\lambda \varphi \left(\varphi ^2\mu ^2R^2\right)=0.$$
(3)
In the broken phase, ($`\mu ^2>0`$) there are domain wall solutions in which the field changes vacuum over a distance of order $`M^1`$, passing through zero. The solution for an infinite static planar wall is well known VilShe94  and with an appropriate choice of coordinates can be written

$$\varphi =\mu \mathrm{tanh}(Mz).$$
(4)
At high temperature the mass parameter receives thermal corrections from the fluctuations in fields to which $`\varphi `$ is coupled Kirzhnits:1972iw ; Kirzhnits:1976ts :

$$\mu ^2(T)=\mu _0^2cT^2,$$

where $`c`$ is a model-dependent numerical factor. For the pure scalar field theory $`c=1/6`$. Thus as the Universe cools the field undergoes a phase transition from $`\varphi =0`$ to $`\varphi =\pm \mu (T)`$. As this transition happens at a finite rate, the field cannot select the same minimum everywhere at the same time, and the Universe divides into domains in which $`\varphi `$ takes either positive or negative values. By continuity of the field, these domains must be separated by domain walls Kibble:1976sj  in which the field approximates the configuration given by Eq. 4 in the transverse direction. The initial size of the domains $`\widehat{\xi }`$ is controlled by the rate at which the transition occurs and how strongly damped the field is Zurek:1996sj .

The domain walls are mostly in the form of one infinite boundary separating percolating clusters of the two vacua Ryden:1989vj ; Coulson:1996nv ; Larsson:1997sp . The subsequent evolution of the field is controlled by the dynamics of this infinite domain wall. The wall has a tension of order $`M^3`$, and tries to straighten out and lose energy, and in finite volume eventually one or other vacuum will take over the whole space and the field thereby reaches equilibrium. One way of quantifying the approach to equilibrium is to measure the area density of the domain wall $`A`$, or equivalently the curvature scale of the wall $`\xi =1/A`$, where Hindmarsh:1996xv

$$A=\delta (\varphi )\left|\varphi \right|.$$
(5)
At late times, this quantity can only depend on time $`t`$ and the mass scale $`M`$. The fact that the wall obeys a Nambu-Goto equation independently of $`M`$ VilShe94  indicates that $`A=a/t`$ purely on dimensional grounds, where $`a`$ is a constant amplitude. This can be called the naive or “classical” scaling hypothesis for domain walls, which has been put on a more rigorous footing in Hindmarsh:1996xv ; Hin02 . The scaling hypothesis, and the theory of dynamic scaling, can be also be applied to other extended topological defects such as cosmic strings. More generally, one can define a scaling exponent $`b`$, such that $`A=a/t^b`$, and one of the goals of this paper is to measure both the exponent $`b`$ and the amplitude $`a`$ as accurately as possible.

The first numerical simulations of this system were performed by Press, Ryden and Spergel Press:1989yh ; Ryden:1989vj  who noted that in comoving coordinates the width of the wall shrinks as $`R^1`$, and so any numerical simulation with a lattice spacing fixed in comoving coordinates runs the risk of failing to resolve the domain wall at late times. They showed that ignoring the $`R`$ dependence of the $`\mu ^2`$ parameter did not substantially affect the dynamics of the walls, and we adopt the same approach. We also include a damping term to simulate cooling in the early stages of the evolution. Written in first order form, the field equations become

$`\dot{\varphi }`$ $`=`$ $`\pi `$ (6)
$`\dot{\pi }`$ $`=`$ $`^2\varphi \lambda \varphi \left(\varphi ^2\mu ^2\right)\left({\displaystyle \frac{\dot{R}}{R}}+\eta \right)\pi ,`$ (7)
The equations may be rescaled to $`\mu =1`$, $`\lambda =1`$, so that the width of the wall is 1. These classical equations can be thought of as effective equations representing the long-wavelength dynamics of the quantum field when the occupation number is high.

The results of Press:1989yh ; Ryden:1989vj  seemed to show that a good fit for the area scaling law was given by

$$A\mathrm{ln}(t)/t,$$
(8)
in both 2 and 3 dimensions, which is not predicted by the theory of dynamic scaling for domain walls Hindmarsh:1996xv ; Hin02 . It is therefore important to check these numerical simulations, and after a decade of development in computer technology one can do much larger simulations in order to eliminate transient effects. Indeed, our largest simulation is performed on 3D grid of $`4096^3`$, which gives us a dynamic range of roughly three orders of magnitude between $`t`$ and $`M^1`$, the characteristic response time of the field 111Such large simulations are possible only by having part of the simulation volume in memory at any one time, and keeping the rest on disk..

## III Numerical Simulations

### III.1 Evolution Algorithm

We used 2 and 3D cubic lattices with a 2nd order discretization of the Laplacian operator. The evolution of the discretized system was effected with the Verlet or leapfrog algorithm, common in molecular dynamics and offering a simple but effective way of discretizing the equations of motion (7). The scheme is

$$\begin{array}{ccc}\hfill \pi _{n+\frac{1}{2}}& =& \pi _{n\frac{1}{2}}+f(t,\varphi _n,\pi _n)\mathrm{\Delta }t\hfill \\ \hfill \varphi _{n+1}& =& \varphi _n+\pi _{n+\frac{1}{2}}\mathrm{\Delta }t,\hfill \end{array}$$
(9)
where $`f(t,\varphi _n,\pi _n)`$ is the right hand side of Eq. (7), and $`\pi _n=(\pi _{n+\frac{1}{2}}\pi _{n\frac{1}{2}})/2`$. As $`f`$ is linear in $`\pi _n`$, the equation can easily be rearranged so that $`\pi _{n+\frac{1}{2}}`$ is on the left hand side. The simulations on which the data is based all have $`\mathrm{\Delta }x=0.3`$ and $`\mathrm{\Delta }t=0.1`$.

We chose periodic boundary conditions, and restricted the length of the simulation to a maximum time equal to $`T_{tot}=N\mathrm{\Delta }x/2\mathrm{\Delta }t,`$ the time required for two signals emitted from the same point and travelling in opposite directions to interfere with each other. This represents the time it takes for the field to “notice” the finite dimension of the lattice it resides in.

### III.2 Initial Conditions

The main objective of the simulations was not to see the formation of the network, but the evolution at later times, and so it is not necessary to start with a proper thermal distribution. Instead, we try to imitate the configuration of the field in the potential at some temperature above the critical $`T>T_C`$, with a Gaussian distribution around zero in $`k`$-space. We created a Gaussian distribution by adding 12 uniformly distributed random numbers of unit variance, and algorithm that is quick and easy to parallelize. By setting $`\varphi _{k=0}=0`$ we set the average of the distribution to zero and the proceed with transforming the field to the $`x`$-representation. This is an alternative to simulating the actual phase transition by deforming the symmetry-breaking potential. The Gaussian distribution was chosen to give a pointwise spatial variance of about $`10^1\mu `$, and the value of the volume-averaged field was always less than 10$`^8.`$ Bigger ranges for the Gaussian distribution were found to lead to instabilities unless high dissipation was used. This is due to the occasional large values of the field, which have high energy densities due to the quartic potential. After initialisation the field is left to roll slowly to the two minima and start oscillating in them. A typical evolution of $`\left|\varphi \right|^2`$ can be seen in Fig. 1.

If the average field were not close to zero in the initial conditions distribute we would have found that the biased initial conditions would give a different time dependence in the area density Coulson:1996nv ; Larsson:1997sp . Indeed, if one phase is selected in the initial distribution then the evolution of the domain wall area density can be shown to follow the relation

$$A\frac{1}{t}e^{ct},$$
(10)
where $`c`$ a positive real number. Percolation theory predicts that for a phase occupying a fraction less than a critical value $`p_c,`$ the infinite wall disappears. The dominant phase takes over at a characteristic time scale, and quickly fills the entire simulation volume. In our case $`p0.5\pm 10^7,`$ far above the percolation threshold which for a cubic lattice can be shown to be $`p_c=0.31.`$

### III.3 Beginning the Evolution

In the beginning of the simulation dissipation is imposed in order for the field to sink in a controlled manner into the minima and to remove spurious high frequency modes. This is controlled by the parameter $`\eta `$ in (7). The dissipation aids the formation of the domain wall network and as soon as this is formed, dissipation is turned off and the system is left to evolve freely. A small amount of time, roughly the system’s characteristic time, is required for the field to adjust itself into the new equation with no dissipation and continue to evolve. For that reason any regression on data starts a few timesteps after the end of the dissipated period. Turning the dissipation off slowly seems to help the system adjust more quickly to the new conditions and smooths out the effects of this ‘adjustment’ period. More specifically, the end of the dissipation period $`T_D`$ and the beginning of the regression $`T_R`$ are specified in the program and their difference $`T_RT_D`$ is calculated. As soon as the system reaches $`T_D,`$ dissipation at timestep $`n`$ $`\eta _n`$ is decreased till zero using a Gaussian like function

$$\eta _n=\eta _\mathrm{i}\times \mathrm{exp}\left(\frac{A\left(nT_D\right)^2}{\left(T_RT_D\right)^2}\right)\times \mathrm{\Theta }(T_RT_D),$$
(11)
where $`n`$ is the timestep, $`\eta _\mathrm{i}`$ the initial dissipation and $`A`$ a real positive number controlling in more accuracy the length of this period.

Dissipations ranged from 0.1 to 0.3 in the simulations and affected the system for roughly 10%-30% of the total simulation time $`T_{tot},`$ with regression on the results starting at 20%-40% of the total time and the coefficient $`A`$ was taken to be in the range $`1`$ to $`1.5`$ .

### III.4 Wall Area Density Calculation

The main purpose of the simulations was to check the power law for the evolution of the wall area density. A crude way of calculating the total wall area is counting the places where adjacent lattice points have opposite signs and multiplying by $`\mathrm{\Delta }x`$ or $`\mathrm{\Delta }x^2`$ which is a rough estimation of the wall area at the “link”. More precisely, one should find an accurate discretization of the continuum area density operator (5).

The calculation is made by finding two adjacent lattice points where the field has opposite signs and calculating the gradient of the field at the “link”. For two such adjacent points $`\{i,j,k\}`$ and $`\{i1,j,k\}`$ the gradient is calculated as follows. The numerical approximation to the gradient at the direction of the link is just

$$\mathrm{\Delta }_i\varphi =\frac{\varphi _{i,j,k}\varphi _{i1,j,k}}{\mathrm{\Delta }x}$$

whereas the gradients at the remaining two directions are taken by averaging over the gradients of nearby links; the $`j`$ -component of the gradient for example would be

$$\mathrm{\Delta }_j\varphi =\frac{1}{2}\left(\frac{\varphi _{i1,j+1,k}\varphi _{i1,j1,k}}{2\mathrm{\Delta }x}+\frac{\varphi _{i,j+1,k}\varphi _{i,j1,k}}{2\mathrm{\Delta }x}\right).$$

One needs to account for the orientation of the area element at each link, otherwise one will end up overestimating the area Ryden:1989vj ; Scherrer:1998sq . For domain walls in 3D the problem can be resolved by simply multiplying the lattice area estimate by a factor $`\frac{2}{3}`$ Scherrer:1998sq , while a similar argument to that given in Scherrer:1998sq  gives the factor $`\pi /4`$ in 2D:

$`A`$ $`=`$ $`A_{Lat}{\displaystyle \frac{2}{3}}\text{(3D)}`$ (12)
$`A`$ $`=`$ $`A_{Lat}{\displaystyle \frac{\pi }{4}}\text{(2D)}`$ (13)
The Figures show $`A_{Lat}`$, while the Tables show $`A`$.

## IV Results

### IV.1 Wall area density

In all simulations the wall area data started to be taken shortly after the dissipation had been stopped. Fig. 2 shows the results from a 2D run, fitted to a

power law

$$A_{Lat}=at^b$$
(14)
by regression. Five simulations for the same parameters were run and an average was taken, with the fit taken on the averaged area. The results for 2D and 3D simulations are presented in Tables 1, 2, for Minkowski,

radiation and matter-dominated Friedmann-Robertson-Walker backgrounds respectively. The biggest 3D simulation ($`N=4096`$, $`\mathrm{\Delta }x=0.3`$, $`\mathrm{\Delta }t=0.08`$) gave $`a=0.980(\pm 0.017)`$, $`b=0.985(\pm 0.003).`$

There is a suggestion from the simulations that the wall area decreases slightly slower than $`b=1.`$ However, this deviation from the $`b=1`$ scaling could not be attributed to a relation of the form given in Eq. 8, as it has been suggested Press:1989yh . Plotting $`\mathrm{exp}(A_{Lat})\times t`$ against time shows that there seems to be no logarithmic term in the wall area evolution, Fig. 3.

Both the power law and the coefficient $`a`$ of Eq. 14. present challenges to the analytic method for the calculation of the wall area density of Ref. Hindmarsh:1996xv , which are compared in detail in Hin02 . The power law is very close to the predicted value of $`b=1`$, with good precision, but the coefficient $`a`$ showed larger fluctuations between runs.

## V Conclusions

The numerical integration of the equations of motion for the the $`\varphi ^4`$ model has given an insight to the dynamics involved in domain wall networks and has provided an accurate way to support the scaling solution predicted by theoretical computations. A power law with exponent very close to 1 was found to be the best solution according to the simulation data with no evidence for a logarithmic term suggested in previous studies Press:1989yh .

The fact that a domain wall network shows this dynamic scaling over approximately three orders of magnitude in the parameter $`M\xi `$, the ratio between the local curvature radius and the wall width, is remarkable, and has important implications for cosmic string networks. Firstly, it is clear from the visualisations in Fig. 4 (see also Sims )

that the energy in the domain walls is very quickly transferred into propagating modes of the field: the formation and collapse of closed loops (2D) or surfaces (3D), expected in the standard picture of the evolution of wall networks VilShe94 , is rare. Indeed, in 2D, self-intersections to form closed loops must be very rare; if two segments of wall approach each other they must be generically curved away from the point of closest approach, and therefore the acceleration is in the direction which would tend to increase the separation. Nonetheless, 2D walls scale perfectly well, so it seems plausible in that case that energy is being transferred directly into radiation. If the amplitude of the oscillations is large enough they could appear to form tiny loops or “protoloops” Moore:1998gp ; Moore:2001px .

We believe that our results add weight to the contention, first put forward in Vincent:1998cx , that extended defects (including cosmic strings in 3D) have a non-perturbative channel into massive radiation. At first sight this is difficult to square with the standard picture, in which walls and strings obey the Nambu-Goto equations of motion for large curvature radii, for in that case the total energy locked up in the extended defects is conserved, in the absence of an general relativistic effects such as an expanding background or gravitational radiation. It is certainly true that it is possible to find string trajectories which are very close to being solutions of the Nambu-Goto equations Moore:1998gp ; Olum:2000sg : however, the initial conditions have to be carefully prepared, and the existence of these trajectories does not preclude the existence of a non-perturbative radiative process for defect networks. Indeed, we maintain that our results are good evidence that there must be such a process.

###### Acknowledgements.

This work was conducted on the SGI Origin platform using COSMOS Consortium facilities, funded by HEFCE, PPARC and SGI. We also acknowledge support from the Sussex High Performance Computing Initiative.

## Appendix A Lattice Correction factor for 2D wall length density

In Ref. Scherrer:1998sq  it was shown that the naive lattice estimate of the area density of random domain walls in 3D overcounts the continuum value by a factor of 3/2 in the limit that the radii of curvature are large compared with the lattice spacing. In this section we perform the analagous calculation for two dimensions, finding it to be $`4/\pi `$.

The naive estimate is obtained by summing the length of all links containing a wall and dividing by the total volume. A link crossing a wall is defined to be one for which the values of the field $`\varphi `$ on the sites at either end have opposite signs. One can immediately see this will overestimate the length, as one is approximating a smooth curve by a sequence of line segments parallel to the lattice vectors $`𝐢`$ and $`𝐣`$.

In the continuum the centre of a domain wall is described by a 2-dimensional curve $`𝐫(\sigma )`$. Consider a segment of length $`l`$, with $`l\xi `$, where $`\xi `$ is the correlation length of the curve. Writing $`\mathrm{\Delta }𝐫=𝐫(l)𝐫(0)=x𝐢+y𝐣`$, we see that the lattice approximation to the length is

$$l_{\mathrm{Lat}}=|x|+|y|.$$
(15)
In the limit that $`l/\xi 0`$, the continuum value of the length is $`l=\sqrt{(x^2+y^2)}`$. Hence

$$l_{\mathrm{Lat}}=\frac{l_{\mathrm{Lat}}}{l}l=(|\mathrm{cos}\alpha |+|\mathrm{sin}\beta |)l$$
(16)
where $`\alpha `$ and $`\beta =\pi /2\alpha `$ are the direction cosines of the vector $`\mathrm{\Delta }𝐫`$. The ratio between the lattice estimate of the length and the true length is obtained by averaging over all possible orientations of the lattice relative to $`\mathrm{\Delta }𝐫`$:

$$\frac{l_{\mathrm{Lat}}}{l}=\frac{1}{2\pi }_0^{2\pi }𝑑\alpha \left(|\mathrm{cos}\alpha |+|\mathrm{sin}\alpha |\right)=\frac{4}{\pi },$$
(17)
as advertised.

This calculation can be extended to arbitrary $`l`$ through a more involved argument. Let us first define the correlation function

$$C(r)=\varphi (𝐫)\varphi (0),$$
(18)
which is assumed to be smooth at $`𝐫=0`$, so that

$$C(r)=C(0)\frac{1}{2}C^{\prime \prime }(0)r^2.$$
(19)
In Refs. OhtJasKaw82 ; Hindmarsh:1996xv  it is shown that the length density $`A`$ of the locus of zeroes of a Gaussian random field in 2 dimensions is given by

$$A=\frac{1}{2}\sqrt{\frac{C^{\prime \prime }(0)}{C(0)}}.$$
(20)
This is the value of the length density in the continuum.

On the lattice, we must consider the probability that the values of the field at opposite ends of a link have opposite signs, in which case we can say that the link is occupied by a segment of domain wall. Let us call this probability $`p_{\mathrm{occ}}`$, which is

$$p_{\mathrm{occ}}=2P(\varphi (𝐱)>0\mathrm{AND}\varphi (𝐱+𝐢\mathrm{\Delta }x)<0),$$
(21)
where the factor of 2 accounts for the opposite case $`\varphi (𝐱)>0\mathrm{AND}\varphi (𝐱+𝐢\mathrm{\Delta }x)<0`$. The lattice estimate of the length density is then

$$A_{\mathrm{Lat}}=\frac{2p_{\mathrm{occ}}}{\mathrm{\Delta }x},$$
(22)
where the lattice spacing is $`\mathrm{\Delta }x`$, and the factor 2 comes from the fact that there are twice as many links as sites in 2 dimensions.

Suppose we now define

$$P_{12}(C(r_{12}),V)=P(\varphi (𝐱_1)>V\mathrm{AND}\varphi (𝐱_2)<V),$$
(23)
where $`V`$ is an arbitrary threshold and $`r_{12}=|𝐱_1𝐱_2|`$. One can then almost trivially write

$$P_{12}(C(r_{12}),V)=_0^{r_{12}}\frac{P_{12}}{C(r)}\frac{dC(r)}{dr}𝑑r.$$
(24)
It can be shown HamGotWei86  that

$$\frac{P_{12}}{C_{12}}=\frac{1}{2\pi [C(0)^2C_{12}^2]^{\frac{1}{2}}}\mathrm{exp}\left(\frac{V^2}{C(0)+C_{12}}\right),$$
(25)
where $`C_{12}=C(r_{12})`$. Hence the lattice estimate of the area density is

$$A_{\mathrm{Lat}}=\frac{2}{\pi \mathrm{\Delta }x}\mathrm{cos}^1\left(\frac{C(\mathrm{\Delta }x)}{C(0)}\right).$$
(26)
Providing $`\mathrm{\Delta }x`$ is much smaller than the correlation length $`\xi `$, defined by $`\xi ^2=|C(0)/C^{\prime \prime }(0)|`$, we see that

$$A_{\mathrm{Lat}}=\frac{2}{\pi }\sqrt{\frac{C^{\prime \prime }(0)}{C(0)}},$$
(27)
and hence, from Eq. (20), that $`A_{\mathrm{Lat}}=4A/\pi `$.