text:
# Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection.

## 1 Introduction

In general, the production of astrophysical jets has been primarily explained by two possible mechanisms: magnetic forces associated with an accretion disk (Königl & Pudritz, 1999) and the interaction between an accretion disk and a magnetized star, in the case of protostellar jets (Shu et al., 1999). The first mechanism is a disk-wind driven model, where the gas is centrifugally flung out from the surface of a Keplerian disk when field lines thread the disk at an angle of $`60^{}`$ or less with respect to the surface. When the gas reaches the local Alfvén velocity, the ($`𝑱_𝒑\times 𝑩_\mathit{\varphi }`$) term becomes dominant, providing the collimation for the outflow. The second mechanism is the X-wind model, where the wind emanates near the region where the disk corotates with the central magnetic star. It is assumed that some magnetospheric lines, connecting the central star and the disk, are partially open. Excess angular momentum brought in by the disk is extracted by the magnetocentrifugally driven wind. The wind emanates from the inner edge of the disk and flows out along the open magnetic field lines. Since the first mechanism seeks to explain astrophysical jets independently of the detailed nature of the central object, from active galactic nuclei (AGNs) to young stellar objects (YSOs), it has the possibility of being a universal model.

Recent theories and models of astrophysical jets have mainly tried to provide a natural explanation for common features. These include the episodic outbursts and knots that are observed.

Numerical simulations of axisymmetric magnetized flows from accrection disks have previously been made using different models. Shibata & Uchida (1986) and Stone & Norman (1994) followed the internal dynamics of the disk. Bell & Lucek (1995) and Ustyugova et al. (1995) treat the disk as a boundary condition at the base of the wind, without following its internal dynamics, and assume it to be in Keplerian equilibrium. Following the latter model, Ouyed & Pudritz (1997a, b; 1999) made 2.5-dimensional, high-resolution, numerical magnetohydrodynamic (MHD) simulations of the onset and collimation of outflows from the surface of a Keplerian accretion disk around a central object. Using two different topologies for the magnetic field that threads the disk and the corona, they obtained steady or episodic jets, depending on the values of the constant (in time) ejection velocity from the disk to the corona. For the smallest values of the ejection velocity, their results showed episodic production of knots near the surface of the disk which propagate along the jet axis. They argued that these knots arose due to MHD shocks.

Krasnopolsky et al. (1999) made time-dependent, axisymmetric simulations of the magnetocentrifugal model jet formation. Their system consisted of a Keplerian disk and a central object. The main subtlety resides in the treatment of the disk-wind boundary, the disk surface. Some, but not all, of the flow variables can be fixed on this boundary. The boundary conditions imposed on the disk surface were limited to those necessary and sufficient to take into account information propagating upstream from the fast and Alfvén critical surfaces, avoiding overdetermination of the flow and the production of impulsive accelerations. They obtained a cold and steady jet, launched smoothly from the Keplerian disk and collimated by the magnetocentrifugal mechanism. The collimation is observed both in the shape of the field lines and in the density profiles, which become cylindrical, and thus jetlike, at large distances along the rotation axis. The structure of the flow is insensitive to the density or flow speed at the injection surface, as long as the mass flux is kept constant. As explained by the authors, steady-state jets obtained in earlier studies, using an essentially incomplete treatment of the disk boundary conditions, are qualitatively unchanged when this deficiency is rectified. It remains to be shown, however, whether the same is true for nonsteady jets.

Frank et al. (2000) made a series of simulations intended to address the issue of MHD jet propagation. Whereas some studies have used ad hoc initial conditions, they injected flows into a computational grid derived from models of collimated jets driven by magnetocentrifugal launching. The initial configurations in the jet were taken directly from the solution of the force balance perpendicular (the Grad-Shafranov equation) and parallel (the Bernoulli equation) to the magnetic surfaces generated by a magnetized rotator. Their simulations followed the evolution of jets composed of helical fields, embedded in hypersonic plasmas whose density and velocity varied with radius. They studied the propagation of jets driven by: (1) a purely Keplerian rotator; and (2) a Keplerian rotator with a sub-Keplerian boundary layer. Both rotators are exterior to a solid body rotator. Studying the adiabatic case for the last model (their multicomponent model) they obtained a jet in which strong instabilities developed in the core of the jet. Once the core is exposed, periodic pinches can be seen in the beam. Their stability analysis indicates that the origin of these periodic pinches could be current driven instabilities. The pinching instabilities on the axis have a wavelength of approximately half of the jet radius. However, the beam also has a second characteristic wavelength which runs along the surface of the core. This feature, which appears as an envelope encompassing the shorter wavelength modes, has a wavelength of approximately three jet radii. As explained by the authors, these results suggest that these instabilities could be at the origin of the knotty structure of a large number of jets as seen, for example, in HL Tau, HH1, HH30, and HH34.

Time-dependent MHD calculations of a modified version of the X-wind model, considering the star-disk field interaction, shows twisting of closed stellar field lines, current sheet formation, reconnection, and ejection of magnetic island with the plasma heated to $`10^8`$ K, consistent with X-ray flare observations (Hayashi et al., 1996; Goodson et al., 1997).

All the above previous studies assumed constant velocity injection from the disk. Here, we extend the previous investigations and study the effect of non-steady injection, in particular, periodic injection.

Adopting the disk as a boundary condition, we investigated how periodic perturbations at the disk surface can alter the formation and the propagation of the magneto-centrifugal jet. In the FU-Orionis phase of a protostar, for example, we know that there is a periodic ejection of material observed in the X-ray range. We do not include the physics of the disk that could generate a pulsed gas ejection. We assume that there is a periodic mass loading at the base of the corona of the disks, which introduces a periodic ejection of gas from the disk surface. Whereas Ouyed & Pudritz (1997a, b; 1999) used a constant ejection velocity from the disk, we study the temporal dependence of this velocity. We made 2.5-dimensional numerical simulations using the ZEUS-3D code (Stone & Norman, 1992a, b) in the axisymmetry option. A sinusoidal perturbation was introduced in the ejection velocity as a temporal function. Various values for the periods in the range from 10 to 80 rotation periods of the innermost radius of the accretion disk ($`r_i`$), were used in the simulations. The maximum value for the perturbation was chosen. This value is roughly equal to that which would produce a steady state jet if the gas were continuously ejected from the disk (Ouyed & Pudritz, 1997a, hereafter OP97a). The luminosity of the protostellar T Tauri stars (which have accrection disks) fluctuates rapidly with large amplitude. Similarly, quasars, which are also assumed to have accretion disks, also have luminosities which fluctuate rapidly with large amplitude. Therefore, we introduced in our simulations velocity perturbations that are a large percentage of the injection velocity. A region near the central object, $`80r_i`$ in the $`z`$-direction and $`20r_i`$ in the $`r`$-direction, was simulated.

Using different values for the period of the sinusoidal perturbation, we obtained the formation of regularly spaced structures along the jet axis. For small values of the period, the structures tend to dissipate along the jet axis; for medium values, they tend to persist and for large values, tend to fragment into smaller substructures.

In §2 the basic equations that are solved are presented. The model is described in §3 and the numerical results are presented and discussed in §4. Finally, in §5 we present a summary and the final conclusions.

## 2 Basic equations

Adopting cylindrical coordinates ($`r`$, $`\varphi `$, $`z`$), we placed the central object at the origin and took the $`z`$-axis to be perpendicular to the disk. The surface of the disk thus lies in the $`z=0`$ plane of our coordinate system.

The equations that describe the conservation of mass, momentum, energy and magnetic flux are:

$$\frac{\rho }{t}+\left(\rho 𝒗\right)=0,$$
(1)
$$\rho \left[\frac{𝒗}{t}+\left(𝒗\right)𝒗\right]+\left(p+p_A\right)+\rho \mathrm{\Phi }𝒋\times 𝑩=0,$$
(2)
$$\rho \left[\frac{e}{t}+\left(𝒗\right)e\right]+p\left(𝒗\right)=0,$$
(3)
$$\frac{𝑩}{t}\times \left(𝒗\times 𝑩\right)=0,$$
(4)
$$𝑩=0,$$
(5)
with $`\rho `$ the gas density, $`𝒗`$ the velocity, $`𝑩`$ the magnetic field, $`p`$ the gas pressure, $`e`$ the internal energy, $`𝒋=\times 𝑩/4\pi `$ the electric current density, and $`\mathrm{\Phi }=GM/\left(r^2+z^2\right)^{1/2}`$ the gravitational potential.

A polytropic gas $`p=K\rho ^\gamma `$ with the polytropic index $`\gamma =5/3`$ is assumed. Thus, ignoring radiative transfer effects, we do not solve the energy equation (3). We also assume that $`p_A`$ is an Alfvénic turbulent pressure. This provides, together with the thermal pressure, the main support against gravity for the cold corona. Here $`p_A=p/\beta _T`$, where $`\beta _T`$ is the ratio of the gas pressure to the Alfvénic turbulent pressure, assumed to be constant (OP97a).

All parameters are given in units of their values at the innermost radius of the accretion disk $`r_i`$. Thus, the normalized equation of motion (eq. ) is

$$\left[\frac{𝒗^{}}{t}+\left(𝒗^{}^{}\right)𝒗^{}\right]=\frac{1}{\delta _i}\frac{^{}\left(p^{}+p_A^{}\right)}{\rho ^{}}^{}\mathrm{\Phi }^{}+\frac{2}{\delta _i\beta _i}\frac{𝑱^{}\times 𝑩^{}}{\rho ^{}},$$
(6)
where the dimensionless variables are: $`r^{}=r/r_i`$, $`z^{}=z/r_i`$, $`v^{}=v/v_{K,i}`$, $`\rho ^{}=\rho /\rho _i`$, $`p^{}=p/p_i`$, $`B^{}=B/B_i`$, $`\mathrm{\Phi }^{}=1/\left(r^2+z^{}_{}{}^{}2\right)^{1/2}`$ and $`^{}=r_i`$. The parameter $`\delta _i=\rho _iv_{K,i}^2/p_i`$ is the ratio of the Keplerian to the thermal energy density and $`\beta _i=8\pi p_i/B_i^2`$ is the ratio of the gas to the magnetic pressure in the corona at $`r_i`$. At $`r_i`$, the Keplerian speed is $`v_{K,i}=\left(GM/r_i\right)^{1/2}`$ and $`B_i`$ is the poloidal field at $`r_i`$. For convenience, we drop the primes and refer only to the normalized variables. Time is given in units of $`t_i=r_i/v_{K,i}`$, so that the dimensionless time is

$$\tau =\frac{t}{t_i}.$$
(7)
## 3 The Model

In this model, the disk is a fixed boundary condition in pressure equilibrium with the corona. At the center there is a point mass that represents a star or a black hole. An initial potential (zero current) poloidal field threads the disk and the corona.

The computational domain defines the spatial region where the MHD equations are time evolved. ZEUS-3D, a grid code, divides the computational domain into cells, called “active” zones. Finite differencing the evolution near the boundaries of the grid requires that values for the dependent variables be specified beyond the computational domain. Therefore, at each boundary, “ghost” zones are added. Values for the dependent variables in the ghost zones are specified, using boundary conditions appropriate to the geometry and the physics of the problem being solved. Thus, the evolution equations are not solved for the ghost zones.

We use inflow boundary conditions at the disk surface and outflow boundary conditions on the remaining boundaries. Reflecting boundaries are used along the axis of symmetry, which coincides with the disk axis. The normal component (radial component) of the velocity and the magnetic field are set to zero on the boundary and reflected into the “ghost” zones. In addition, the 3-component (toroidal component) of the velocity and magnetic field is also inverted along the axis of symmetry into the “ghost” zones. The initial setup and boundaries are shown in Figure Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection. (Fig. 2 of OP97a).

In the corona, the initial toroidal field is zero. Magnetostatic equilibrium is numerically maintained to assure that an outflow can only arise due to the magnetocentrifugal mechanism. Resolving equation (6) for magnetostatic equilibrium, we obtain the expression for the coronal density distribution,

$$\rho =\left(1/\sqrt{r^2+z^2}\right)^{3/2},$$
(8)
which is used in our numerical simulations.

Since the surface of the accretion disk is a fixed boundary condition, there is no back reaction by the jet. We are, thus, time evaluating the coronal region ($`z0`$) only where the jet arises. However, we fix the density at the base of the corona (first column of cells), assuming that the disk is a suitable dynamic source of matter. The density distribution of the disk surface has the form

$$\rho _d=\eta _ir^{3/2},$$
(9)
where $`\eta _i=\rho _d/\rho _0`$ is the ratio between the disk density ($`\rho _d`$) and the density at the base of the corona ($`\rho _0\rho (z=0)`$).

The corona and the disk are threaded by the poloidal magnetic field in a potential configuration ($`𝑱=0`$), described by the stream function,

$$\mathrm{\Phi }_m=\sqrt{r^2+\left(z_d+z\right)^2}\left(z_d+z\right),$$
(10)
where $`z_d`$ is the dimensionless disk thickness. The field components are given as

$$B_z=\frac{1}{r}\frac{\mathrm{\Phi }_m}{r}\mathrm{and}B_r=\frac{1}{r}\frac{\mathrm{\Phi }_m}{z}.$$
(11)
Due to the torsion of the poloidal field at the disk surface, a toroidal magnetic field is produced at the base of the corona, which could propagate into the disk. We thus expect that the disk has a nonzero toroidal field,

$$B_\varphi =\frac{\mu _i}{r},$$
(12)
where the input parameter $`\mu _i=B_{\varphi i}/B_i`$ is the ratio of the toroidal to poloidal magnetic field at $`r_i`$ . There are three important initial timescales defined at the disk’s surface: the Kepler time $`\tau _K`$, the Alfvén time $`\tau _A`$ and the magnetic braking time $`\tau _B`$. Thus, we have (e.g., eq. \[4.64\] of OP97a)

$$(\tau _K,\tau _A,\tau _B)=(1,\sqrt{\delta _i\beta _i/2},\eta _i\sqrt{\delta _i\beta _i/2}).$$
(13)
To estimate the turbulent pressure needed to support the corona for a given gas pressure, we use the equation for the ratio of the turbulent to magnetic pressure $`\beta _T`$

$$\beta _T=\gamma /\left(\left(\gamma 1\right)\delta _i\gamma \right),$$
(14)
where, as in equation (6), $`\delta _i`$ is the ratio of the Kepler to thermal energy density. The gas is injected from the disk into the corona at a velocity $`𝒗_p=f_vv_K𝑩_p/B_p`$, where $`v_K=\left(1/r\right)^{1/2}`$ is the local Keplerian velocity at the disk surface (for $`r1`$)111For the radial region where the disk does not exist ($`r<1`$ (i.e., $`r<r_i`$)) we set $`𝒗_p`$ to zero., $`𝑩_p`$ is the local poloidal magnetic field that threads the disk, and $`f_v`$ is the temporal function that describes the temporal dependence of the injection velocity. In our simulation we introduced a pulsed injection velocity with a sinusoidal form,

$$f_v=v_{inj0}\left[1+\mathrm{sin}\left(2\pi \tau /T\right)\right],$$
(15)
where $`T`$ is the period of the perturbation, $`v_{inj0}`$ is an input parameter and $`\tau `$ is given by equation (7).

We, thus, have six free parameters in our simulation: $`\delta _i`$, $`\beta _i`$, $`\eta _i`$, $`\mu _i`$, $`v_{inj0}`$ and $`T`$.

## 4 Numerical results and discussion

We used the following parameters to investigate the time evolution of a jet from a Keplerian accretion disk with pulsed gas injection (similarly to OP97a):

$$(\delta _i,\beta _i,\eta _i,\mu _i,v_{inj0})=(100.0,1.0,100.0,1.0,10^3),$$
(16)
implying the initial timescales

$$(\tau _K,\tau _A,\tau _B)=(1.0,7.07,707.1).$$
(17)
Five simulations were made. In simulation A, there was no perturbation; in simulations B to E, we introduced perturbations with various values for the period $`T`$. The simulations were made in the domain $`(z,r)=(80,20)`$ with a resolution of $`(120,45)`$ for 400 time units $`\left(<\tau _B\right)`$. Along the z-axis, 30 zones span a distance of $`10r_i`$ outwards from the disk surface. An additional 90 ratioed zones span an additional distance of $`70r_i`$, with each subsequent zone increasing in size by a factor of 1.017. Twenty uniform zones span an initial five disk radius in the radial direction. Outside of these uniform grid zones, 25 zones are ratioed, with each subsequent zone increasing in size by a factor 1.065.

Figure Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection. shows the initial conditions in the atmosphere. The density distribution (eq. ) is shown at the top and the initial magnetic (potential) configuration, at the bottom of the figure. Initially, the velocity and the toroidal magnetic field are zero.

The results of the simulations are shown in a series of figures that are described below. Figures Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection. to Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection. show the distributions of density, poloidal velocity (vectors), poloidal and toroidal magnetic fields and toroidal velocity at 400 time units for the simulations A to E, respectively. The contour lines are isocontours in the $`rz`$ plane.

In simulation A (Fig. Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection.), no perturbation was introduced in the gas injection. Results in good agreement with those of OP97a were obtained. From the density distribution, shown at the top of this figure, we observe the formation of a highly collimated steady state jet along the z-axis. The collimation is due to the strong radial pinch force of a dominant toroidal field (fourth panel from top to bottom) that propagates into the corona. We obtain a maximum value for the poloidal velocity (second panel from top to bottom) on the order of the Keplerian velocity at $`r_i`$.

A sinusoidal perturbation in the gas injection velocity (eq. ) at the disk surface was introduced in the other four simulations. As we are interested in studying the behaviour of the jet subject to this perturbation at its base, we chose a progressive range for the period T. In simulations B to E, we use $`T=10,20,40,80`$, respectively.

Along the jet axis in simulation B for $`T=10`$ (Fig. Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection.), we observe the formation of regularly spaced structures up to about $`z=50r_i`$. The distance of separation between two consecutive structures is about $`5r_i`$ (i.e., $`\left(T/2\right)r_i`$). The axial extent of the regularly spaced structures advance with the flow speed until it reaches the limiting distance of about $`50r_i`$. This pattern persists even for times greater than 400 time units.

Our results were obtained using ratioed zones that allow for a smaller number of calculated zones and, therefore, a poorer resolution as compared with OP97a. To verify the sensitivity of our results to this poorer resolution, we made simulation C using only uniform zones at a resolution close to that used in the uniform zone region in the poorer resolution case. This simulation was run in the domain $`(z,r)=(80,20)`$ with the resolution of $`(250,80)`$ uniform zones. For the others cases (D – E), we used the poorer resolution.

Regularly spaced structures along the whole jet (Fig. Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection.) are seen in simulation C, where $`T=20`$. These structures persist during the entire simulation, even for times greater than 400 units. The distance of separation between two consecutive structures is about $`11r_i`$, or roughly $`\left(T/2\right)r_i`$. We obtained these same results when using ratioed zones in the poorer resolution case.

We observe the formation of regularly spaced structures along the entire jet as well as their fragmentation in simulation D for $`T=40`$. Few minor structures are observed between the major structures (Fig. Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection.: top). Fragmentation of the structures can be seen in the isocontours of the toroidal velocity in the bottom panel of Figure Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection.. These structures persist during the entire simulation, even for times greater than 400 units. The distance of separation between two major structures is about $`22r_i`$, i.e., roughly $`\left(T/2\right)r_i`$.

In simulation E for $`T=80`$, we observe the formation of regularly spaced structures along the entire jet (Fig. Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection.: top), with more pronounced fragmentation clearly seen in the plot of the toroidal velocity (Fig. Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection.: bottom). These structures persist during the entire simulation, even for times greater than 400 units. The distance of separation between two major structures is about $`40r_i`$, i.e., $`\left(T/2\right)r_i`$.

The axial distributions, for all the simulations, of density, axial velocity, toroidal velocity and toroidal magnetic field along a cut at $`r=5r_i`$, parallel to the disk’s axis at $`400`$ time units, are shown in Figure (Fig. Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection.). These quantities are normalized to their maximum values along this cut. We note from the axial velocity distribution (dotted line), that the flow is rapidly accelerated along the field line from the disk to roughly $`z=15r_i`$. After this point there is a “Hubble flow” constant acceleration, with $`v_z`$ proportional to $`z`$ for simulation A. However, for simulations B to E this constant acceleration character is modified, reflecting the pulsed nature of the jet. From the toroidal magnetic field distribution (dashed-dotted line) we observe the effect of torsion of the poloidal field lines due to the rotation of the disk. The regular structures, in simulations B – E, have their minimum mass concentration (solid line) at the maximum values of toroidal magnetic field (dashed-dotted line), reflecting jet constriction by this field. The shift between the toroidal magnetic field and the density in the pinched structures is a general result of long standing (e.g., Chan & Henriksen (1980)) and it has been associated with the periodic oscillations caused by the pinch instability due to the toroidal magnetic field.

For simulations A to E at $`\tau =400`$, the axial and radial mass fluxes are shown in Figure Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection.. The mass flux is defined as

$$\dot{m}=\rho 𝒗𝒅𝑨,$$
(18)
where $`A`$ is the area. The axial mass flux222The computation of the axial mass flux involves computing the product of $`v_z`$ and the area of each annulus for each radial zone $`[\pi (r_{j+1}^2r_j^2)]`$, where $`j`$ is an active zone and $`r_j`$ refers to the position of the face or edge of the zone. At each $`z`$, the axial mass flux is computed from the sum of the axial mass flux contributions of all 45 radial zones. refers to the quantity of mass per unit time which passes through a circular area $`A=\pi r^2`$ and is calculated for each $`z`$, using $`r=20r_i`$; the radial mass flux333The radial mass flux at each $`z`$ involves computing the product of $`v_r`$ and the area of the annulus for a given radial zone $`[2\pi r_j(z_{k+1}z_k)]`$, where $`j`$ and $`k`$ define each active zone, and $`r_j`$ and $`z_k`$ refers to the face/edge position of the zone, along the $`r`$ and $`z`$ direction, respectively. In particular, we computed the radial flux at each $`z`$ using $`r=20r_i`$. refers to the quantity of mass per unit time which passes through the area $`A=2\pi rz`$ and is shown for each z, using $`r=20r_i`$. These fluxes, given in units of $`\dot{m}_i=2\pi r_i^2\rho _iv_{Ki}`$, are important for the detection of the decollimation of the jet. Their values agree with the results of OP97a. At the time $`\tau =400`$, the axial fluxes at $`z=80r_i`$ are comparable to those at $`z0`$; the small difference owing to decollimation of the jet. As the axial mass flux decreases, the radial flux can be seen to increase by roughly the same amount. We can see the pulsed nature of the jet (simulations B to E) from the axial mass fluxes (solid lines) which oscillate around their average values. The radial mass fluxes in simulations B to E are similar to those of simulation A (without perturbation). In simulation B, the axial mass flux curve is not very modified by the density structures, as compared with the case for which there is no perturbation (simulation A); in simulations C to D, the axial mass flux curves are strongly modified by the density structures. Their average values, however, do not change.

Figure Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection. shows the sonic, Alfvénic and magnetosonic Mach numbers for all the simulations at $`\tau =400`$, in a cut at $`r=5r_i`$ along the $`z`$-axis. In the computation of the sonic Mach numbers, we use only the thermal pressure $`p`$ in the sound speed. To compute the magnetosonic Mach number, we use the effective sound speed which includes the thermal and the Alfvénic turbulent pressure, $`p_A`$444The magnetosonic speed is defined as $`\left(c_s^2+v_{Az}^2\right)^{1/2}`$ where $`c_s=\left(\gamma \left(p+p_A\right)/\rho \right)^{1/2}`$ is the effective sound speed and $`v_{Az}=B/\left(4\pi \rho \right)^{1/2}`$ is the axial Alfvén velocity.. We note in Figure Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection. that the jet is already supersonic and super-Alfvénic near the base of the corona and super-magnetosonic at $`z4r_i`$. Only where density structures exist are the Mach numbers (Fig. Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection.) modified. Their median values, however, remain roughly the same. The maximum values of Alfvénic Mach numbers (dashed-dotted line) are due to the decreased axial magnetic fields. Negative values correspond to a reversal of the axial magnetic field.

## 5 Summary and Conclusions

We made numerical 2.5-dimensional simulations of astrophysical jets produced from Keplerian disks subject to a sinusoidal perturbation of the gas injection velocity from the disk into the corona. Ratioed zones were used in order to obtain better resolution near the region of the generation of the jet with a poorer resolution outside.

Five simulations were made with different values of the period of the sinusoidal perturbation. No perturbation was introduced in simulation A. In simulations B to E, we used periods $`T=10,20,40`$ and $`80`$, respectively, in units of $`t_i=r_i/v_{Ki}`$, where $`r_i`$ is the innermost radius of the disk and $`v_{Ki}`$, its Keplerian velocity.

For simulation A, where there was no perturbation, we obtained a steady state highly collimated jet with a maximum axial velocity of the order of $`v_{Ki}`$ after a time $`\tau =400`$. In a cut at $`r=5r_i`$ parallel to the disk’s axis, the jet is supersonic and already super-Alfvénic near the base of the corona and super-magnetosonic at $`z4r_i`$. These results agree with those of OP97a, who used better resolution.

Introducing a sinusoidal perturbation in the gas injection velocity, the results common, to simulations B to E, are:

1. The formation of regularly spaced structures with a separation distance of roughly $`\left(T/2\right)r_i`$ along the jet axis.

2. The toroidal magnetic field is out of phase ($`90^{}`$) with the density curve along the jet axis.

3. The axial mass fluxes are modified by the periodic density structures, although the average remains roughly the same, as compared with the same curve for the case of no perturbation (simulation A). The axial fluxes at $`z=80r_i`$ are comparable to the axial fluxes at $`z0`$.

4. The radial mass fluxes are similar to those of simulation A (where there was no perturbation). However, differences arise when we examine a cut at $`r=5r_i`$. There is oscillatory behaviour in some of the simulations at this cut.

5. In a cut at $`r=5r_i`$ the Mach number distributions are modified only where density structures exist. The median values, however, remain roughly the same, and are similar to simulation A.

The results that are different in some of these simulations are:

a. In simulation B ($`T=10`$) the structures persist only up to an axial distance of $`z50r_i`$.

b. In simulations D ($`T=40`$) and E ($`T=80`$), fragmentation of the structures occur.

These simulations, introducing periodic perturbations in the accretion disk, produced collimated jets with a wavelength of $`\left(T/2\right)r_i`$, where $`T`$ is the temporal period of the perturbation. The creation of structures with a specific distance separation was determined by the period $`T`$ of the gas ejection from the disk surface. For example, we obtained a jet with a distance separation of about $`11r_i`$ associated with a period of 20 time units for the perturbation. However, for different values of $`T`$, an interesting effect occurs: the structures tend to dissipate at large axial distances for smaller values of $`T`$ and to fragment into smaller structures for greater values of $`T`$.

Using this same numerical approach, Ouyed & Pudritz(1999) obtained episodic jets, using steady gas injection from the disk with velocities which were smaller than $`10^3v_{K,i}`$ which we used. For injection velocities equal to this value, they obtained a highly collimated steady state jet555They obtained steady state jets with velocities equal to or greater than $`5\times 10^4v_{K,i}`$, as well as episodic jets with different distances of separation for smaller velocities, using the same free parameters that we used. with features similar to ours, with no perturbation in the velocity field of the disk. When we introduced perturbations with specific features in our steady state jet, our results showed the tendency to produce regular structures with separation distances $`11r_i`$, similar to the case where OP99 used a smaller velocity of $`10^4v_{K,i}`$.

Gardiner & Frank (2000), using a magnetocentrifugal approach for the jet launching, studied an analytical one-dimensional model, in which they investigated the behaviour of a jet with a variable output speed. They inferred that a pulsing jet with an embedded helical magnetic field will inevitably develop a periodic field structure, consisting of rarefied regions (where the helix is “combed out”, resulting in a more poloidal geometry) alternating with denser, toroidally dominated knots. To confirm these predictions, they performed axisymmetric (2.5-dimensional) numerical simulations of a radiative MHD jet. They used conditions appropriate for protostellar jets and introduced a pulsed velocity with a sinusoidal form at the jet nozzle. Their results showed that pulsing jets with initially helical fields will evolve into an alternating poloidal-toroidal-dominated field geometry. The strongest toroidal field are confined to the dense knots, while predominantly poloidal fields are in the rarefied regions.

Both our results and those of Gardiner & Frank (2000) show that episodic jets are produced by introducing a pulsed mass load at the base of the jet. In particular, in our simulations, as well as in the study made by Ouyed & Pudritz (1997b), knot formation is associated with an anticorrelation of density and toroidal field strength. However, it should be noted that our study is concerned with the near-field region (i.e., the field region near the origin of the jet), whereas other studies, such as Gardiner & Frank (2000), are concerned with the far-field region.

There is currently some debate over the origin of the knots in protostellar flows. It remains unclear as to whether the knots are due to the pulsing of the source of the jet or to instabilities in the beam. Our results show that for a given constant ejection velocity from the disk surface which produces a steady state jet, an episodic jet is produced when the given ejection velocity is pulsed. We showed that the temporal and spatial knotty structure of the jets depends on the period of the velocity perturbation at the base of the jet.

Ouyed & Pudritz (1997b, 1999) rule out the possibility of a pinch MHD instability causing the knotty structure. This instability depends upon the strength of $`\beta `$ and no such dependence was found by the authors. They found that only for $`\beta _i`$ around unity, were knots produced. Birkinshaw (1990) noted that the Kelvin-Helmholtz (KH) instability does not exist if the beam speed is less than the Alfvén speed. Hardee & Rosen (1999) showed that KH growth rates increase as the magnetosonic Mach number decreases, provided the jet is super-Alfvénic and the magnetized jet is nearly stabilized to the KH instability when the jet is sub-Alfvénic. Hardee & Rosen (2002) found that the jets can be stalilized to the KH helical and higher order asymmetric normal modes provided the velocity shear, $`\mathrm{\Delta }_uu_ju_e`$ between the jet and the external medium is less than a “surface” Alfven speed, $`v_{AS}\left[\left(\rho _J+\rho _E\right)\left(B_J^2+B_E^2\right)/\left(4\pi \rho _J\rho _E\right)\right]^{1/2}`$, where $`\rho _J(\rho _E)`$ is the jet (external medium) density and $`B_J(B_E)`$ is the jet (external medium) magnetic field. Ouyed et al. (2003) showed that the jet beyound the Alfvén surface becomes unstable to nonaxisymmetric KH instabilities, in three-dimensional simulations of jets from Keplerian accretion disks.

We also made a number of simulations using values of $`\beta _i`$ different from unity. Figure Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection. shows the results of one of these simulations, a periodic perturbation at the disk surface, using a period of $`T=20`$. We plotted the isocontours of the density. The results for $`\beta _i=0.1`$ are shown at the top of the figure. For this case, we see a smoothing of the structures along the jet. The isocontours of the density oscillate slightly near the jet axis. At the botton of this figure, the results for $`\beta _i=10`$ are shown. Here we observe the formation of structures from the jet axis up to a distance of about $`z=35r_i`$. The region near the jet axis shows irregular structures and is less well defined (compare with Figure Numerical Simulations of Astrophysical Jets from Keplerian Disks with Periodic Ejection., where $`\beta _i=1`$).

Our results are consistent with those of Hardee & Rosen (1999, 2002) and Ouyed et al. (2003). This is shown in the sequence of simulations in Figures 5 and 11. When $`\beta _i`$ is small (strong magnetic fields, and the flow is sub-Alfvénic), there is a not much growth of the instability. The super-Alfvénic case in which $`\beta _i`$ is close to unity (the flow speed is close to the Alfvén speed) is more unstable.

We simulated a region very near the central object ($`20r_i\times 80r_i`$). It is interesting to evaluate the spatial dimensions involved in our simulations. For a protostar of mass $`0.5M_{}`$, taking $`r_i=0.04`$ AU for the inner radius of the disk, we have $`80r_i10^5`$ pc. For a black hole of $`10^8M_{}`$ (the Schwartzschild radius is $`1`$ AU), assuming $`r_i=20.6`$ AU, we have $`80r_i10^2`$ pc. A distance of separation of $`11r_i`$ of the structures corresponds, for the protostar case, to $`10^6`$ pc and for the black hole case, to $`10^3`$ pc. These distances are much smaller than those between the Herbig-Haro objects observed in protostellar jets (about $`10^2`$ pc). They are also much smaller than the distances observed between the emission knots, frequently found in AGN jets. In contrast to other simulations (e.g., Frank et al.(2000)), we obtain structures very near the formation region of the jets, involving distances not yet resolved by current instruments.

We can not say, at the present moment, whether or not the origin of the knotty structures seen in many jets at large distances from the source are due to the structures observed in the present simulations for regions close to the source. Further studies are needed to investigate this question.

The instabilities seen by Ouyed & Pudritz (1997b, 1999), producing the knotty structures, occur at the fast magnetosonic surface and propagate downstream. These instabilities, as explained by the authors, are due to MHD shocks. Some conditions in the knot generating region (in particular, for high $`𝑩_\varphi /𝑩_z`$) seem to point towards a scenario in which the knot spacings could possibly be identified with the resonant wavelength pinch mode. As noted above, Ouyed & Pudritz (1997b, 1999) rule out the possibility of a pinch MHD instability causing the knotty structures since it has a strong $`\beta `$ dependence, which was not found by the authors. The authors only observed knotty structures when $`\beta `$ was around unity.

We believe that the knotty structures that we obtained in our simulations are due to the association of the pulsing mass load at the base of the jet and the instability due to MHD shocks, discussed in the studies of Ouyed & Pudritz (1997b, 1999). Our results indicate that a particular frequency of the pulsed velocity can excite the instability and produce persistent knotty structures along the jet.

The authors would like to thanks the editor, Steven Shore, and the anonymous referee for helpful comments. B. F. V. acknowledges the financial support of the Brazilian agency FAPESP under grant 95/99122. V. J. P. and R. O. would like to thank the Brazilian agency CNPq for partial support. We would like to thank the LCCA-USP for the possibility of using the CRAY J-90 computer in our simulations and especially the staff for their invaluable help. R. Ouyed and especially C. Fendt contributed to this paper with encouraging and valuable discussions, for which we are grateful. The authors also would like to thank the project PRONEX (No. 41.96.0908.00) for partial support. R. O. would like to thank the Brazilian agency FAPESP for partial support under grant 00/06770-2. V. J. P. would like to thank FAPESP for partial support under grant 2000/06769-4.


format: text

text:
# On the Asymptotic Number of Plane Curves and Alternating Knots

## 1. Introduction.

Our motivation for this work is the enumeration of topological equivalence classes of smooth open and closed curves in the plane (see Figure 1; precise definitions are given in Section 2.1).

The problem of characterizing closed curves was considered already by Gauss and has generated many works since then: see  and references therein. Our interest here is in the numbers $`a_p`$ and $`\alpha _p`$ of such open and closed curves with $`p`$ self-intersections, and more precisely we shall consider the asymptotic properties of $`a_p`$ and $`\alpha _p`$ when $`p`$ goes to infinity. The numbers $`a_p`$ were given up to $`p=10`$ in  and have been recently computed up to $`p=22`$ by transfer matrix methods . Asymptotically, as $`p`$ goes to infinity, one expects the relation $`a_p4\alpha _p`$ to hold (see below), so that we concentrate on the numbers $`a_p`$.

In the present paper we propose a physical reinterpretation of the numbers $`a_p`$ that leads to the following conjecture, and we present numerical results supporting it.

###### Conjecture 1.

There exist constants $`\tau `$ and $`c`$ such that

$$\alpha _p\underset{p\mathrm{}}{}\frac{1}{4}a_p\underset{p\mathrm{}}{}c\tau ^pp^{\gamma 2},$$

where

(1)
$$\gamma =\frac{1+\sqrt{13}}{6}0.76759\mathrm{}$$

From the data of  one has the numerical estimate: $`\tau 11.4`$. But the main point in Conjecture 1 lies not so much in the existence of $`\tau `$ as in the explicit value of $`\gamma `$. It should indeed be observed that $`\gamma `$, rather than $`\tau `$, is the interesting information in the asymptotic of $`a_p`$. Observe for instance that the value of $`\gamma `$ is left unchanged if one redefines the size of a closed curve as the number $`p^{}=2p`$ of arcs between crossings. More generally, as discussed in Section 8, the exponent $`\gamma `$ determines the branching behavior of generic large curves under the uniform distribution, and is *universal* in the sense that the same value is expected in the asymptotic of related families of objects like prime self-intersecting curves or alternating knots.

Conjecture 1 is similar in nature to the conjecture of Di Francesco, Golinelli and Guitter on the asymptotic behavior of the number of plane meanders . The two problems do not fall into the same universality class (in particular the predictions for the exponent $`\gamma `$ are different in the two problems). However our approach to design a numerical test is applicable also to the meander problem.

The rest of the paper is organized as follows. Precise definitions are given, and a more general family of drawings is introduced that play an important rôle in the identification of the associated physical model (Section 2). The physical background leading to the conjecture is then discussed (Section 3) and a numerically testable quantity is proposed (Section 4). The sampling method is briefly presented (Section 5) before the analysis of numerical data (Section 6). We conclude with some variants and corollaries of the conjecture (Section 8).

## 2. Plane curves and colored planar maps

### 2.1. Plane curves and doodles

For $`p`$ a positive integer, let $`𝒜_p`$ be the set of equivalence classes of self-intersecting loops $`\gamma `$ in the plane, that is: (i) $`\gamma `$ is a smooth mapping $`S^1^2`$; (ii) there are $`p`$ points of self-intersection, all of which are regular crossings; (iii) two loops $`\gamma `$ and $`\gamma ^{}`$ are equivalent if there exists an orientation preserving homeomorphism $`h`$ of the plane such that $`\gamma ^{}(S^1)=h(\gamma (S^1))`$.

Similarly let $`A_p`$ be the set of equivalence classes of self-intersecting open curves $`\gamma `$ in the plane: (i) $`\gamma `$ is a smooth mapping $`[0,1]^2`$ and $`\gamma (0)`$ and $`\gamma (1)`$ belong to the infinite component of $`^2\gamma ((0,1))`$; (ii) there are $`p`$ points of self-intersection, all of which are regular crossings; (iii) two open curves are equivalent if there exists an orientation preserving homeomorphism $`h`$ of the plane such that $`\gamma ^{}([0,1])=h(\gamma ([0,1]))`$ and $`\gamma ^{}(i)=h(\gamma (i))`$ for $`i=0,1`$.

Observe that, unlike closed curves, open curves are oriented from the initial point $`\gamma (0)`$ to the final point $`\gamma (1)`$. Moreover a unique closed curve is obtained from an open curve by connecting the final point to the initial one in counterclockwise direction around the curve. These definitions are illustrated by Figure 1.

In order to study the families $`𝒜_p`$ and $`A_p`$ and to obtain Conjecture 1 we introduce a more general class of drawings, that we call *doodles*. For given positive integers $`p`$ and $`k`$, let $`A_{k,p}`$ be the set of equivalence classes of $`(k+1)`$-uples $`\mathrm{\Gamma }=(\gamma _0,\gamma _1,\mathrm{},\gamma _k)`$ of curves drawn on the plane such that: (i) the curve $`\gamma _0`$ is an open curve of the plane: $`\gamma _0`$ is a smooth mapping $`[0,1]^2`$, and $`\gamma _0(0)`$ and $`\gamma _0(1)`$ belong to the infinite component of $`^2(\gamma _0((0,1))_i\gamma _i(S^1))`$; (ii) for $`i1`$, each $`\gamma _i`$ is a loop, that is a smooth map $`S^1^2`$; (iii) there are $`p`$ points of intersection (including possibly self-intersections) of these curves, all of which are regular crossings; (iv) the union of the curves is connected, (v) two doodles $`\mathrm{\Gamma }=(\gamma _0,\mathrm{},\gamma _k)`$ and $`\mathrm{\Gamma }^{}=(\gamma _0^{},\mathrm{},\gamma _k^{})`$ are equivalent if there exists an orientation preserving homeomorphism $`h`$ of the plane such that $`\gamma _0^{}([0,1])_i\gamma _i^{}(S^1)=h(\gamma _0([0,1])_i\gamma _i(S^1))`$ and $`\gamma _0^{}(x)=h(\gamma _0(x))`$, for $`x=0,1`$. In other terms, a doodle is made of an open curve intersecting a set of loops, that are considered up to continuous deformations of the plane. An example of doodle is given in Figure 2 (left-hand side).

### 2.2. Colored planar maps

An equivalent presentation of doodles is in terms of rooted planar maps . A *planar map* is a proper embedding of a connected graph in the plane considered up to homeomorphisms of the plane. It is *4-regular* if all vertices have degree four. It is *rooted* if one root edge is marked on the infinite face and oriented in counterclockwise direction. Equivalently the root edge can be cut into an in- and an out-going half-edge (also called *legs*) in the infinite face. There is an immediate one-to-one correspondence between doodles with $`p`$ crossings, and 4-regular planar maps with $`p`$ vertices and two legs. This correspondence is illustrated on Figure 2.

We shall consider the number $`a_{k,p}=\mathrm{card}A_{k,p}`$ of doodles with $`p`$ crossings and $`k`$ loops and more specifically we shall consider the asymptotic properties of $`a_{k,p}`$ as $`p`$ (and possibly $`k`$) goes to infinity. It turns out to be convenient to introduce the generating function $`a_p(n)`$ as $`k`$ varies:

(2)
$$a_p(n)=\underset{k=0}{\overset{\mathrm{}}{}}a_{k,p}n^k$$

The requirement that a doodle is connected implies that it cannot contain more loops than crossings so that $`a_p(n)`$ is a polynomial in the (formal) variable111The subsequent interpretation in terms of colored doodles and the strong tradition in the physics litterature are our (admittingly poor) excuses for the use of $`n`$ to denote a formal variable. $`n`$. For real valued $`n`$, $`a_p(n)`$ can be understood as a weighted summation over all doodles with $`p`$ crossings, and, more specifically for $`n`$ a positive integer, $`a_p(n)`$ can be interpreted as the number of *colored* doodles in which each loop has been drawn using a color taken from a set of $`n`$ distinct colors.

On the one hand, in the special case $`k=0`$, $`a_{0,p}=a_p`$ gives by definition the number of open self-intersecting plane curves. Observe that $`a_p`$ is also given by $`n=0`$ since $`a_p(0)=a_{0,p}`$. On the other hand, the generating functions of the $`a_p(n)`$ for other values of $`n`$, namely $`n=1,2,2`$, have been computed exactly (see respectively ). We elaborate now on the case $`n=1`$ since it will play a crucial role in what follows.

The number $`a_p(1)`$ counts the number of doodles with $`p`$ crossings irrespective of the number of loops $`k`$. In terms of maps, $`a_p(1)`$ is the number of rooted 4-regular planar maps with $`p`$ vertices. The number of such planar maps is known , from which one can compute the asymptotics:

(3)
$$a_p(1)=2\frac{3^p(2p)!}{p!(p+2)!}\underset{p\mathrm{}}{}\frac{2}{\sqrt{\pi }}\mathrm{\hspace{0.33em}12}^pp^{5/2}.$$

Observe in this case the power-law exponent $`5/2`$, which is *universal* for rooted planar maps in the sense that it is observed for a variety of families of rooted planar maps (see ). As opposed to this, the exponential growth factor $`12`$ is specific to the family of rooted 4-regular planar maps.

There is a physical interpretation of the power-law behavior $`p^{5/2}`$: it is given by two-dimensional gravity. This explanation begs to be generalized to any $`n`$, and we shall explore such a possibility now.

## 3. Two-dimensional quantum gravity and asymptotic combinatorics

The purpose of this section is to give the rationale behind our conjectures. We place the discussion at a rather informal level that we hope achieves a double purpose: on the one hand it should give an intuition of the path leading to our conjectures to the reader with zero-knowledge in quantum field theory (QFT), and on the other hand it should convince the expert of this (quantum) field. For our defense, let us observe that filling in more details would require a complete course on QFT, with the result of not getting much closer to a mathematical proof.

### 3.1. From planar maps to two-dimensional quantum gravity

The main idea of the physical interpretation of the numbers $`a_p(1)`$ is to consider planar maps as discretized random surfaces (with the topology of the sphere). As the number of vertices of the map grows large, the details of the discretization can be assimilated to the fluctuations of the metric on the sphere. To make this idea more precise, let us describe a way to associate a metric on the sphere to a given 4-regular map $`m`$: to each vertex of $`m`$ associate a unit square and identify the sides of these squares according to the edges of $`m`$ (arbitrary number of corners of squares get identified); the result is by construction a metric space with the topology of a sphere. Upon taking a 4-regular map uniformly at random in the set of map with $`p`$ edges, a random metric sphere with area $`p`$ is obtained.

Now, physics tells us that the metric is the dynamical field of general relativity *i.e.* gravity, and that this type of fluctuations in the metric are characteristic of a quantum theory. In our case it means that, as $`p`$ becomes large, the discrete nature of the maps can be ignored and there exists a scaling limit, the properties of which are described by two-dimensional euclidian quantum gravity. In particular any parameter of random planar maps that makes sense in the scaling should converge to its continuum analog. A fundamental parameter of this kind turns out to be precisely the number of (unrooted) planar maps: it is expected to scale to the partition function $`Z_\mathrm{g}(A)`$ of two-dimensional quantum gravity with spherical topology at fixed area $`A`$, through a relation of the form

(4)
$$\frac{1}{p}a_p(1)\underset{p\mathrm{}}{}Z_\mathrm{g}(A),\text{ with }A\text{ proportional to }p.$$

(Here the factor $`1/p`$ is due to the fact that the partition function does not takes the rooting into account.) The only thing we want to retain from $`Z_\mathrm{g}(A)`$ is that the power law dependance of its large area asymptotic takes the form $`A^{7/2}`$, in accordance with Formula (3). (Trying to give here a precise description of the partition function $`Z_\mathrm{g}`$ would carry us to far away, and anyway the arguments in this section are non rigourous.)

In the case $`n=1`$, this is the whole physical picture: a fluctuating but empty two-dimensional spacetime – there is no matter in it. What happens when $`n1`$? As already discussed, an appealing image is to consider that one must “decorate” the planar map by coloring each curve $`\gamma _i`$ with $`n`$ colors. Alternatively, the physicist’s view would be to consider that we have put a statistical lattice model (of crossing loops) on a random lattice (the planar map). This perfectly fits with the previous interpretation of the planar map as a fluctuating two-dimensional spacetime: as we learn from physics, in the limit of large size, adding a statistical lattice model amounts to coupling matter to quantum gravity. Matter is described by a quantum field theory (QFT) living on the two-dimensional spacetime. The parameters of the lattice model that survive in the scaling limit are recovered in the critical (long distance) behavior of this QFT, which in turn is described by a conformal field theory (CFT). Then, provided we can find a CFT describing the lattice model corresponding to a given $`n1`$, the analog of Relation (4) holds with the partition function $`Z_{\mathrm{g}+\mathrm{CFT}(n)}(A)`$ of this CFT coupled to gravity: in the large size limit,

(5)
$$\frac{1}{p}a_p(n)\underset{p\mathrm{}}{}Z_{\mathrm{g}+\mathrm{CFT}(n)}(A).$$

In this picture, the only thing we need to know about the CFT that describes the scaling limit of our model is its *central charge* $`c`$, which roughly counts its number of degrees of freedom. Indeed, the study of CFT coupled to gravity was performed in , resulting in the following fundamental prediction: the partition function $`Z_{\mathrm{g}+\mathrm{CFT}}(A)`$ of gravity dressed with matter has a power-law dependence on the area of the form $`A^{\gamma 3}`$ where the critical exponent $`\gamma `$ depends only on the central charge of the underlying CFT via (for $`c<1`$)

(6)
$$\gamma =\frac{c1\sqrt{(1c)(25c)}}{12}.$$

Returning to our asymptotic enumeration problem (not forgetting the extra factor $`p`$ which comes from the marked edge), we find:

(7)
$$a_p(n)\underset{p\mathrm{}}{}\mathrm{e}^{\sigma p+\left(\gamma 2\right)\mathrm{log}p+\kappa }$$

where $`\sigma `$, $`\kappa `$ are unspecified “non-universal” parameters, whereas the “universal” exponent $`\gamma `$ is given by Eq. (6) with the central charge $`c`$ of the a priori unknown underlying CFT$`(n)`$. The absence of matter, that is the case $`n=1`$, corresponds to a CFT with central charge $`c=0`$: one recovers $`\gamma 2=5/2`$ as expected from Eq. (3). In general, all parameters in Eq. (7) are functions of $`n`$; assuming furthermore that their dependence on $`n`$ is smooth in a neighborhood of $`n=1`$, one can recover by Legendre transform of $`\sigma (n)`$ the asymptotics of $`a_{k,p}`$ as $`k`$ and $`p`$ tend to infinity with the ratio $`k/p`$ fixed. Observe finally that the knowledge of the CFT could give more informations. For instance, the irrelevant operators of the CFT control subleading corrections to Eq. (7).

### 3.2. The identification of two candidate models

We now come to the issue of the determination of the CFT for a arbitrary $`n`$. An observation made in , based on a matrix integral formulation, is that this CFT must have an $`O(n)`$ symmetry (for $`n`$ positive integer – for generic $`n`$ this symmetry becomes rather formal and cannot be realized as a unitary action of a compact group on the Hilbert space). There exists a well-known statistical model with $`O(n)`$ symmetry, a model of (dense) non-crossing loops , whose continuum limit for $`|n|<2`$ is described by a CFT with central charge

(8)
$$c_\mathrm{I}=16(\sqrt{g}1/\sqrt{g})^2n=2\mathrm{cos}(\pi g),0<g<1$$

In  it was speculated that there is no phase transition between the model of non-crossing loops, which we call model I, and our model of crossing loops, and therefore the central charge is the same and given by Eq. (8). If this were the case, the study of irrelevant operators of this CFT would allow moreover to predict that subleading corrections to Eq. (7) have power-law behavior for all $`|n|<2`$ with exponents depending continuously on $`n`$.

However, another scenario is possible. In , it was suggested that $`O(n)`$ models, for $`n<2`$, possess in general a low temperature phase with spontaneous symmetry breaking of the $`O(n)`$ symmetry into a subgroup $`O(n1)`$. This is a well-known mechanism222The Mermin–Wagner theorem, which forbids spontaneous symmetry breaking of a continuous symmetry in two dimensions, only applies to $`n`$ integer greater or equal to $`2`$. in QFT (see e.g.  chapters 14, 30), which produces Goldstone bosons living on the homogeneous space $`O(n)/O(n1)=S^{n1}`$. In the low energy limit the bosons become free and the central charge is simply the dimension of the target space $`S^{n1}`$:

(9)
$$c_{\mathrm{II}}=n1n<2$$

For generic real $`n`$ ($`n<2`$) this is only meaningful in the sense of analytic continuation, but we assume it can be done and call it model II. This CFT possess a marginally irrelevant operator, leading to main corrections to leading behavior (7) of logarithmic type *i.e.* in $`\frac{1}{\mathrm{log}p}`$, $`\frac{\mathrm{log}\mathrm{log}p}{(\mathrm{log}p)^2}`$ etc.

It was furthermore argued in  that the critical phase of model II is generic in the sense that the low-energy CFT is not destroyed by small perturbations – the most relevant $`O(n)`$-invariant perturbation is the action itself, which corresponds to a marginally irrelevant operator for $`n<2`$. On the contrary, the model I of non-crossing loops is unstable to perturbation by crossings; some numerical work on regular lattices (at $`n=0`$)  tends to suggest that it flows towards model II.

Note that both Conjectures (8) and (9) supply the correct value $`c=0`$ for $`n=1`$ and $`c=1`$ for the limiting case $`n=2`$.333Actually, the two resulting $`c=1`$ theories are not identical: the one from model I seems to be the wrong one, although this is a subtle point on which we do not dwell here. Of course, in no way do we claim that these are the only possible scenarios which fit with known results – one might have a plateau of non-critical behavior ($`c=0`$) around $`n=1`$, for instance; or two-dimensional quantum gravity universality arguments might not apply at all in some regions of $`n`$ – but they seem the most likely candidates and therefore it is important to find a numerically accessible quantity which at least discriminates between the two conjectures.

## 4. The general conjectures and a testable parameter

The physical reinterpretation of doodles as a model on random planar lattices has led us to postulate that the weighted summation over doodles satisfies

$$a_p(n)c_0(n)\tau (n)^pp^{\gamma (n)2},$$

with the critical exponent $`\gamma (n)`$ given in terms of the central charge $`c(n)`$ by

$$\gamma (n)=\frac{c(n)1\sqrt{(1c(n))(25c(n))}}{12}.$$

Moreover we have presented two concurrent models which fix the value of $`c(n)`$. Since negative values of $`n`$ create additional technical difficulties (appearance of complex singularities in the generating function, cf ), we formulate the conjectures in the restricted range $`0n<2`$:

###### Conjecture 2 (Model I).

Colored doodles are in the universality class of dense non-crossing loops, so that for $`0n<2`$, $`n=2\mathrm{cos}(\pi g)`$, $`1/2g<1`$,

$$c(n)=16(\sqrt{g}1/\sqrt{g})^2.$$

###### Conjecture 3 (Model II).

Colored doodles are in the universality class of models with spontaneously broken $`O(n)`$ symmetry, so that for $`0n<2`$,

$$c(n)=n1.$$

Observe that Conjecture 3 implies Conjecture 1 for $`n=0`$, while Conjecture 2 would give $`c(0)=16(\sqrt{2}1/\sqrt{2})^2=2`$ and $`\gamma (0)=1`$. According to the discussion of the previous section, Conjecture 3 appears more convincing. In order to get a numerical confirmation, we look for a way to discriminate between the two.

Since the model at $`n=1`$ is much easier to manipulate, we look for such a quantity at $`n=1`$. Of course the known value of the exponent $`\gamma (1)`$ is a natural candidate but as already mentioned both conjectures agree on this: we propose instead the derivative of the exponent at $`n=1`$,

(10)
$$\gamma ^{}\frac{d}{dn}_{|n=1}\gamma (n).$$

The reason that it can easily be computed numerically is that it appears in the expansion of the average number of loops $`k_p`$ for a uniformly distributed random planar map with $`p`$ vertices. Indeed one easily finds

(11)
$$k_p=\frac{d}{dn}_{|n=1}\mathrm{log}a_p(n)\underset{p\mathrm{}}{=}\sigma ^{}p+\gamma ^{}\mathrm{log}p+\kappa ^{}+o(1)$$

Here we have assumed expansion (7) to be uniform with smoothly varying constants $`\sigma (n)`$, $`\gamma (n)`$, $`\kappa (n)`$ in some neighborhood of $`n=1`$, and written $`\sigma ^{}\frac{d}{dn}_{|n=1}\sigma (n)`$, $`\kappa ^{}\frac{d}{dn}_{|n=1}\kappa (n)`$.

The conjectures 2 and 3 provide the following predictions for $`\gamma ^{}`$:

(12)
$$\gamma ^{}=\{\begin{array}{cc}\frac{3\sqrt{3}}{4\pi }=0.413\mathrm{}\hfill & \text{in CFT I}\hfill \\ \frac{3}{10}=0.3\hfill & \text{in CFT II}\hfill \end{array}$$

The quantity $`k_p`$ is not known theoretically, so that we cannot immediately conclude in either direction. However it is possible to estimate it numerically using random sampling.

## 5. Sampling random planar maps

In this section we present the algorithm we use to sample a random map from the uniform distribution on rooted 4-regular planar maps with $`p`$ vertices. The problem of sampling random planar maps with various constraints under the uniform distribution was first approached in mathematical physics using markov chain methods . However these methods require a large and unknown number of iterations, and only approximate the uniform distribution. Another approach was proposed based on the original recursive decompositions of Tutte  but has quadratic complexity , and is limited as well to $`p`$ of order a few thousands.

We use here a more efficient method that was proposed in  along with a new derivation of Tutte’s formulas. The algorithm, which we outline here in the case of 4-regular maps, requires only $`O(p)`$ operations to generate a map with $`p`$ vertices and manipulates only integers bounded by $`O(p)`$. Moreover maps are sampled exactly from the uniform distribution. The only limitation thus lies in the space occupied by the generated map. In practice we were able to generate maps with up to 100 million vertices, with a generation speed of a million vertices per second.

The algorithm relies on a correspondence between rooted 4-regular planar maps and a family of trees that we now define. A *blossom tree* is a planted plane tree such that

* vertices of degree one are of two types: buds and leaves;
* each inner vertex has degree four and is incident to exactly one bud;
* the root is a leaf.

An example of blossom tree is shown on Figure 3. By definition a blossom tree with $`p`$ inner vertices has $`p+2`$ leaves (including the root) and $`p`$ buds. Observe that removing the buds of a blossom tree gives a planted complete binary tree with $`p`$ inner vertices, and that conversely $`3^p`$ blossom trees can be constructed out of given binary tree with $`p`$ inner vertices. Since the number of binary trees with $`p`$ inner vertices is well known to be the Catalan number $`\frac{1}{p+1}\left(\genfrac{}{}{0pt}{}{2p}{p}\right)`$, the number of blossom tree is seen to be

$$3^p\frac{1}{p+1}\left(\genfrac{}{}{0pt}{}{2p}{p}\right).$$

Let us define the *closure* of a blossom tree. An example is shown on Figure 3. Buds and leaves of a blossom tree with $`p`$ inner vertices form in the infinite face a cyclic sequence with $`p`$ buds and $`p+2`$ leaves. In this sequence each pair of consecutive bud and leaf (in counterclockwise order around the infinite face) are merged to form an edge enclosing a finite face containing no unmatched bud or leaf. Matched buds and leaves are eliminated from the sequence of buds and leaves in the infinite face and the matching process can be repeated until there is no more buds available. Two leaves then remain in the infinite face.

###### Proposition 1 ().

Closure defines a $`(p+2)`$-to-$`2`$ correspondence between blossom trees and rooted four-regular planar maps. In particular the number of rooted four-regular planar maps is

$$\frac{2}{p+2}\frac{3^p}{p+1}\left(\genfrac{}{}{0pt}{}{2p}{p}\right).$$

This proposition implies that to generate a random map according to the uniform distribution on rooted 4-regular planar maps with $`p`$ vertices one can generate a blossom tree according to the uniform distribution on blossom tree and apply closure. A synopsis of the sampling algorithms is given below. An implementation is available on the web page of G.S.

Random sampling of a rooted 4-regular maps with $`p`$ vertices.

*Step 1.* Generate a random complete binary tree $`T_1`$ according to the uniform distribution on complete binary trees with $`p`$ inner vertices. (This is done in linear time using *e.g.* prefix codes .)

*Step 2.* Convert $`T_1`$ into a random blossom tree $`T_2`$ from the uniform distribution on blossom trees with $`p`$ inner vertices: independantly add a bud on each vertex in a uniformly chosen position among the three possibilities.

*Step 3.* Use a stack (a.k.a. a last-in-first-out waiting line) to realise the closure of $`T_2`$ in linear time: Perform a counterclockwise traversal of the infinite face until $`p`$ buds and leaves have been matched; when a bud is met, put $`b`$ into the stack; when a leaf $`\mathrm{}`$ is met and the stack is non empty, remove the last bud entered in the stack and match it with $`\mathrm{}`$.

*Step 4.* Choose uniformly the root between the two remaining leaves.

## 6. Simulation results

The algorithm described in the previous section allows to generate random rooted 4-regular planar maps with $`p`$ vertices and two legs, with uniform probability. One can compute various quantities related to the map thus generated and then average over a sample of maps, as always in Montecarlo simulations. Here the main quantity of interest is the number of loops of the map. If we generate $`N`$ maps of size $`p`$ so that the $`i`$th map has $`k_{p,i}`$ loops, $`1iN`$, then $`\frac{1}{N}_{i=1}^Nk_{p,i}`$ has an expectation value of $`k_p`$ and a variance of $`\frac{1}{N}k^2_p`$, where $`k_p=\frac{d}{dn}_{|n=1}\mathrm{log}a_p(n)`$ and $`k^2_p=\frac{d^2}{dn^2}_{|n=1}\mathrm{log}a_p(n)`$ (the latter can of course itself be estimated as the expectation value of $`\frac{1}{N1}_{i=1}^Nk_{p,i}^2\frac{1}{N(N1)}(_{i=1}^nk_{p,i})^2`$). According to Eq. (7), both $`k_p`$ and $`k^2_p`$ are of order $`p`$ for $`p`$ large. However, we are interested in corrections to the leading behavior of $`k_p`$ which are of order $`\mathrm{log}p`$, cf Eq. (11), so that we need to keep the absolute error small. This implies that the size of the sample $`N`$ should scale like $`p`$, or that the computation time grows quadratically as a function of $`p`$.

In practice we have produced data for $`p=2^{\mathrm{}}`$ with $`\mathrm{}24`$, the sample size being of the order of up to $`10^7`$. To ensure a good sampling we used the “Mersenne twister” pseudo-random generator , which is both fast and unbiased. The last few values of $`\mathrm{}`$ are only given to show where the statistical error begins to grow large due to limited memory and computation time. Let us call $`k_{\mathrm{}}`$ the numerical value found for $`k_{p=2^{\mathrm{}}}`$. The results obtained are shown on table 1.

First, as a rough check of the asymptotic behavior, let us define $`u_{\mathrm{}}=2k_{\mathrm{}}k_{\mathrm{}+1}`$. If expansion (11) is correct, then $`u_{\mathrm{}}`$ must display an affine behavior as a function of $`\mathrm{}`$: $`u_{\mathrm{}}=(\mathrm{}1)\gamma ^{}\mathrm{log}2+\kappa ^{}+O(1/\mathrm{})`$. Indeed, as one can see on Fig. 4, this is the case.

By comparison with the proposed asymptote it seems clear that $`\gamma ^{}`$ is close to $`0.3`$. To make this statement more precise, one can try to fit the set of the $`k_{\mathrm{}}`$ to $`\sigma ^{}p+\gamma ^{}\mathrm{log}p+\kappa ^{}`$, where $`\mathrm{}`$ ranges from $`\mathrm{}=\mathrm{}_{min}`$ to $`\mathrm{}=\mathrm{}_{max}=24`$ and $`\mathrm{}_{min}`$ is varied. The results are reported on Tab. 2. Unfortunately the confidence level remains fairly low until $`\mathrm{}_{min}`$ becomes so high that statistical error is huge, which tends to indicate strong corrections to the proposed fit.

It is important to understand that if Conjecture 2 were true, then the corrections to asymptotic behavior would be power-law – starting with $`p^{1/2}`$. This means that the procedure used in Tab. 2 should converge quickly to the correct values of $`\sigma ^{}`$, $`\gamma ^{}`$, $`\kappa ^{}`$ (to check this we have performed a similar analysis with a model of non-crossing loops on random planar maps and obtained fast convergence with high accuracy – 2 digits on $`\gamma ^{}`$). Here the range of values of $`\gamma ^{}`$ seems to be $`0.29`$$`0.34`$, far from the value predicted by Conjecture 2. It is therefore our view that the numerical data render Conjecture 2 extremely unlikely.

On the other hand, the value $`0.3`$ predicted by Conjecture 3 remains possible. The fluctuations observed even for very high $`p`$ would be caused by the logarithmic corrections present in model II due to the marginally irrelevant operator, as mentioned in Section 3.2. This operator is expected to induce a correction in $`1/\mathrm{log}p`$ (which is in principle computable exactly using quantum field theory techniques, since it is universal; progress on this will be reported elsewhere), plus higher corrections, all of which remain significant in our range of data. This would also explain why it is so hard to extract useful information from the first few (exact) values of $`a_p(n)`$ given in .

In conclusion, and in view of the theoretical as well as numerical evidence, our belief is that Conjecture 3 is indeed correct.

## 7. Variants and corollaries

First observe that planar maps have in general no symmetries. More precisely the fraction of planar maps with $`p`$ edges that have a non trivial automorphism group goes to zero exponentially fast under very mild assumption on the family considered . If this (very plausible) property holds then a typical closed curve will be obtained by closing $`d`$ different open curves, where $`d`$ is the degree of the outer face. But the average degree of faces in any fixed 4-regular planar map is four. Thus the relation $`a_p4\alpha _p`$.

Second let us give a property illustrating the importance of the critical exponent $`\gamma `$ as opposed to the actual value of $`\tau `$. A closed plane curve $`C`$ is said to be *$`\alpha `$-separable*, for $`0<\alpha 1`$ a constant, if there exist two simple points $`x`$ and $`y`$ of $`C`$ such that $`\mathrm{\Gamma }\{x,y\}`$ is not connected and both connected components contain at least $`p^\alpha `$ crossings. The pair $`(x,y)`$ is called a cut of $`C`$. In other terms, $`C`$ is $`\alpha `$-separable if it is obtained by gluing the endpoints of two big enough open plane curves (up to homeomorphisms of the sphere).

###### Corollary 1.

Assume Conjecture 1 is valid, and consider a uniform random closed plane curve $`\mathrm{\Gamma }_p`$ with $`p`$ crossings. The probability that $`\mathrm{\Gamma }_p`$ is 1-separable decays at least like $`p^\gamma p^{0.77}`$. More generally, if $`\alpha >1/(1\gamma )=(7\sqrt{13})/60.56`$, the probability that $`\mathrm{\Gamma }_p`$ is $`\alpha `$-separable goes to zero as $`p`$ goes to infinity.

For comparison, $`\gamma =1/2`$ and $`1/(1\gamma )=2/3`$ for doodles, which are thus easier to separate.

Indeed let us compute the expected number of inequivalent cuts of a closed plane curve with $`p`$ crossings. When considered up to homeomorphisms of the sphere, close plane curves with a marked cut are in one-to-one correspondence with pairs of open plane curves. Hence, with a factor $`p`$ for the choice of infinite face,

(13)
$$p\underset{p^{}=q}{\overset{pq}{}}\frac{a_p^{}a_{pp^{}}}{\alpha _p}<cstp\underset{p^{}=q}{\overset{pq}{}}\frac{(p^{})^{\gamma 2}(pp^{})^{\gamma 2}}{p^{\gamma 2}}=O(pq^{\gamma 1}).$$

In particular if $`qp^{1/(1\gamma )}`$ this expectation goes to zero as $`p`$ goes to infinity.

It is typical that in the computation of probabilistic quantities, like in Equation (13), the exponential growth factors cancel, leading to behaviors that are driven by polynomial exponents. This explains the interest of in these *critical exponents* and gives probabilistic meaning to their apparent *universality*. As a final illustration of this point let us present two variants of Conjecture 1: (Definitions of prime self-intersecting curves and alternating knots can be found in .)

###### Conjecture 4.

The number $`\alpha _p^{}`$ of closed prime self-intersecting curves with $`p`$ crossings and the number $`\alpha _p^{\prime \prime }`$ of prime alternating knots with $`p`$ crossings lie in the same universality class as closed self-intersecting curves: there are constants $`\tau ^{}`$, $`\tau ^{\prime \prime }`$, $`c^{}`$, $`c^{\prime \prime }`$ such that

$$\alpha _p^{}c^{}\tau ^{}{}_{}{}^{p}p^{\gamma 2},\alpha _p^{\prime \prime }c^{\prime \prime }\tau ^{\prime \prime }{}_{}{}^{p}p^{\gamma 3},$$

where $`\gamma `$ is given in Conjecture 1.

Observe that knot diagrams are naturally considered up to homeomorphisms of the sphere , while we have considered plane curves up to homeomorphisms of the plane. This explains the discrepancy of a factor $`p`$ in Conjecture 4 for $`\alpha _p^{\prime \prime }`$, since one of the $`p+2`$ faces of a spherical diagram must be selected to puncture the sphere and put the diagram in the plane.

## 8. Conclusion

We have given arguments supporting Conjecture 1 for the asymptotic number of plane curves with a large number of self-intersections, as well as the more general Conjecture 3. The numerical results provided in Section 6 support Conjecture 1 only indirectly since they are related to another specialization of Conjecture 3 (derivative at $`n=1`$ versus $`n=0`$). However the alternative proposal is not compatible with either of these new numerical results (as is the case of Conjecture 2) or earlier ones.

Our method to test the conjecture could be applied to other models like open curves with endpoints that are not constrained to stay in the infinite face, or the meanders studied by Di Francesco et al.

### Acknowledgements

P.Z.-J. would like to thank J. Jacobsen for pointing out Refs.  and  to him.


text:
# The Branch Set of a QuasiregularMappingSupported by NSF grant DMS 9970427. I thank Mario Bonk and Alex Eremenko for their criticism on earlier versions of this article. My warmest thanks go to Mario Bonk, Seppo Rickman, and Dennis Sullivan for collaboration, mentoring, and friendship.

## 1. Branched coverings

A continuous mapping $`f:XY`$ between topological spaces is said to be a branched covering if $`f`$ is an open mapping and if for each $`yY`$ the preimage $`f^1(y)`$ is a discrete subset of $`X`$. The branch set $`B_f`$ of $`f`$ is the closed set of points in $`X`$ where $`f`$ does not define a local homeomorphism.

Nonconstant holomorphic functions between connected Riemann surfaces are examples of branched coverings. From the Weierstrassian (power series) point of view this property of holomorphic functions is almost immediate. It is a deeper fact, due to Riemann, that the same conclusion can be drawn from the mere definition of complex differentiability, or, equivalently, from the Cauchy-Riemann equations. Most of this article discusses the repercussions of this fact.

## 2. Quasiregular mappings

In a 1966 paper , Reshetnyak penned a definition for mappings of bounded distortion or, as they are more commonly called today, quasiregular mappings. These are nonconstant mappings $`f:\mathrm{\Omega }𝐑^n`$ in the Sobolev space $`W_{loc}^{1,n}(\mathrm{\Omega };𝐑^n)`$, where $`\mathrm{\Omega }𝐑^n`$ is a domain and $`n2`$, satisfying the following requirement: there exists a constant $`K1`$ such that

$$|f^{}(x)|^nKJ_f(x)$$
(2.1)
for almost every $`x\mathrm{\Omega }`$, where $`|f^{}(x)|`$ denotes the operator norm of the (formal) differential matrix $`f^{}(x)`$ with $`J_f(x)`$=det$`f^{}(x)`$ its Jacobian determinant. One also speaks about $`K`$-quasiregular mappings if the constant in (2.1) is to be emphasized.111The definition readily extends for mappings between connected oriented Riemannian $`n`$-manifolds.

Requirement (2.1) had been used as the analytic definition for quasiconformal mappings since the 1930s, with varying degrees of smoothness conditions on $`f`$. Quasiconformal mappings are by definition quasiregular homeomorphisms, and Reshetnyak was the first to ask what information inequality (2.1) harbours per se. In a series of papers in 1966–69, Reshetnyak laid the analytic foundations for the theory of quasiregular mappings. The single deepest fact he discovered was that quasiregular mappings are branched coverings (as defined above). It is instructive to outline the main steps in the proof for this remarkable assertion, which akin to Riemann’s result exerts significant topological information from purely analytic data. For the details, see, e.g., , , .

To wit, let $`f:\mathrm{\Omega }𝐑^n`$ be $`K`$quasiregular. Fix $`y𝐑^n`$ and consider the preimage $`Z=f^1(y)`$. One first shows that the function $`u(x)=\mathrm{log}|f(x)y|`$ solves a quasilinear elliptic partial differential equation

$$\text{div}𝒜(x,u(x))=0,𝒜(x,\xi )\xi |\xi |^n,$$
(2.2)
in the open set $`\mathrm{\Omega }Z`$ in the weak (distributional) sense. In general, $`𝒜`$ in (2.2) depends on $`f`$, but its ellipticity only on $`K`$ and $`n`$. For holomorphic functions, i.e., for $`n=2`$ and $`K=1`$, equation (2.2) reduces to the Laplace equation $``$div$`u=0`$.

Now $`u(x)`$ tends to $`\mathrm{}`$ continuously as $`x`$ tends to $`Z`$. Reshetnyak develops sufficient nonlinear potential theory to conclude that such polar sets, associated with equation (2.2), have Hausdorff dimension zero. It follows that $`Z`$ is totally disconnected, i.e., the mapping $`f`$ is light. This is the purely analytic part of the proof. The next step is to show that nonconstant quasiregular mappings are sense-preserving. This part of the proof mixes analysis and topology. What remains is a purely topological fact that sense-preserving and light mappings between connected oriented manifolds are branched coverings.

Initially, Reshetnyak’s theorem served as the basis for a higher dimensional function theory. In the 1980’s, it was discovered by researchers in nonlinear elasticity. In the following, we shall discuss more recent, different types of applications.

## 3. The branch set

Branched coverings between surfaces behave locally like analytic functions according to a classical theorem of Stoïlow. By a theorem of Chernavskiĭ, for every $`n2`$, the branch set of a discrete and open mapping between $`n`$-manifolds has topological dimension at most $`n2`$. For branched coverings between 3-manifolds, the branch set is either empty or has topological dimension 1 , but in dimensions $`n5`$ there are branched coverings between $`n`$-manifolds with branch set of dimension $`n4`$, cf. Section 7.222See  for a recent survey on dimension theory and branched coverings.

The branch set of a quasiregular mapping is a somewhat enigmatic object in dimensions $`n3`$. It can be very complicated, containing for example many wild Cantor sets of classical geometric topology , . There is currently no theory available that would explain or describe the geometry of allowable branch sets, cf. Problems 2 and 4 in Section 7.

In the next three sections, we shall discuss the problem of finding bi-Lipschitz parametrizations for metric spaces. It will become clear only later how this problem is related to the branch set.

## 4. Bi-Lipschitz parametrization of spaces

A homeomorphism $`f:XY`$ between metric spaces is bi-Lipschitz if there exists a constant $`L1`$ such that

$$L^1d_X(a,b)d_Y(f(a),f(b))Ld_X(a,b)$$

for each pair of points $`a,bX`$. It appears to be a difficult problem to decide when a given a metric space $`X`$ can be covered by open sets each of which is bi-Lipschitz homeomorphic to an open set in $`𝐑^n`$, $`n2`$. If this is the case, let us say, for brevity and with a slight abuse of language, that $`X`$ is locally bi-Lipschitz equivalent to $`𝐑^n`$.

Now a separable metrizable space is a Lipschitz manifold (in the sense of charts) if and only if it admits a metric, compatible with the given topology, that makes the space locally bi-Lipschitz equivalent to $`𝐑^n`$ . The problem here is different from characterizing Lipschitz manifolds among topological spaces, for the metric is given first, cf. , , , .

To get a grasp of the difficulty of the problem, consider the following example: There exist finite $`5`$-dimensional polyhedra that are homeomorphic to the standard $`5`$-sphere $`𝐒^5`$, but not locally bi-Lipschitz equivalent to $`𝐑^5`$. This observation of Siebenmann and Sullivan  is based on a deep result of Edwards , which asserts that the double suspension $`\mathrm{\Sigma }^2H^3`$ of a $`3`$-dimensional homology sphere $`H^3`$, with nontrivial fundamental group, is homeomorphic to the standard sphere $`𝐒^5`$. (See also .) One can think of $`X=\mathrm{\Sigma }^2H^3`$ as a join $`X=𝐒^1H^3`$, and it is easy to check that the complement of the suspension circle $`𝐒^1`$ in $`X`$ is not simply connected. Consequently, every homeomorphism $`f:X𝐒^5`$ must transfer $`𝐒^1`$ to a closed curve $`\mathrm{\Gamma }=f(𝐒^1)`$ whose complement in $`𝐒^5`$ is not simply connected. A general position argument and Fubini’s theorem imply that, in this case, the Hausdorff dimension of $`\mathrm{\Gamma }`$ must be at least 3. Hence $`f`$ cannot be Lipschitz. In fact, $`f`$ cannot be Hölder continuous with any exponent greater than $`1/3`$. It is not known what other obstructions there are for a homeomorphism $`X𝐒^5`$, cf. \[16, Questions 12–14\].

See  and  for surveys on parametrization and related topics.

## 5. Necessary conditions

What are the obvious necessary conditions that a given metric space $`X`$ must satisfy, if it were to be locally bi-Lipschitz equivalent to $`𝐑^n`$, $`n2`$? Clearly, $`X`$ must be an $`n`$-manifold. Next, bi-Lipschitz mappings preserve Hausdorff measure in a quantitative manner, so in particular $`X`$ must be $`n`$-rectifiable in the sense of geometric measure theory; moreover, locally the Hausdorff $`n`$-measure should assign to each ball of radius $`r>0`$ in $`X`$ a mass comparable to $`r^n`$. Let us say that $`X`$ is *metrically $`n`$-dimensional* if it satisfies these geometric measure theoretic requirements.

It is not difficult to find examples of metrically $`n`$-dimensional manifolds that are not locally bi-Lipschitz equivalent to $`𝐑^n`$. The measure theory allows for cusps and folds that are not tolerated by bi-Lipschitz parametrizations. Further geometric constraints are necessary; but, unlike in the case of the measure theoretic conditions, it is not obvious what these constraints should be. A convenient choice is that of local linear contractibility: locally each metric ball in $`X`$ can be contracted to a point inside a ball with the same center but radius multiplied by a fixed factor. 333See  for analytic implications of this condition.

Still, a metrically $`n`$-dimensional and locally linearly contractible metric $`n`$-manifold need not be locally bi-Lipschitz equivalent to $`𝐑^n`$. The double suspension of a homology 3-sphere with nontrivial fundamental group as described in the previous section serves as a counterexample. In 1996, Semmes ,  exhibited examples to the same effect in all dimensions $`n3`$, and recently Laakso  crushed the last hope that the above conditions might characterize at least $`2`$-dimensional metric manifolds that are locally bi-Lipschitz equivalent to $`𝐑^2`$. However, unlike the examples of Edwards and Semmes, Laakso’s metric space cannot be embedded bi-Lipschitzly in any finite dimensional Euclidean space. Thus the following problem remains open:

###### Problem 1

Let $`X`$ be a topological surface inside some $`𝐑^N`$ with the inherited metric. Assume that $`X`$ is metrically 2-dimensional and locally linearly contractible. Is $`X`$ then locally by-Lipschitz equivalent to $`𝐑^2`$?

In conclusion, perhaps excepting the dimension $`n=2`$, more necessary conditions are needed in order to characterize the spaces that are locally bi-Lipschitz equivalent to $`𝐑^n`$.444There are interesting and nontrivial sufficient conditions known, but these are far from being necessary , , , , . The idea to use Reshetnyak’s theorem in this connection originates in two papers by Sullivan , , and is later developed in . Recall that in this theorem topological conclusions are drawn from purely analytic data. Now imagine that such data would make sense in a space that is not a priori Euclidean. Then, if one could obtain a branched covering mapping into $`𝐑^n`$, manifold points would appear, at least outside the branch set. We discuss the possibility to develop this idea in the next section.

## 6. Cartan-Whitney presentations

Let $`X`$ be a metrically $`n`$-dimensional, linearly locally contractible $`n`$-manifold that is also a metric subspace of some $`𝐑^N`$. Suppose that there exists a bi-Lipschitz homeomorphism $`f:Xf(X)𝐑^n`$. Then $`f`$ pulls back to $`X`$ the standard coframe of $`𝐑^n`$, providing almost everywhere defined (essentially) bounded differential 1-forms $`\rho _i=f^{}dx_i`$, $`i=1,\mathrm{},n`$. To be more precise here, by Kirzsbraun’s theorem, $`f`$ can be extended to a Lipschitz mapping $`\overline{f}:𝐑^N𝐑^n`$, and the 1-forms

$$\rho _i=\overline{f}^{}dx_i=d\overline{f}_i,i=1,\mathrm{},n,$$
(6.1)
are well defined in $`𝐑^N`$ as flat 1-forms of Whitney. Flat forms are forms with $`L^{\mathrm{}}`$-coefficients such that the distributional exterior differential of the form also has $`L^{\mathrm{}}`$-coefficients. The forms in (6.1) are closed, because the fundamental relation $`d\overline{f}^{}=\overline{f}^{}d`$ holds true for Lipschitz maps.

According to a theorem of Whitney \[45, Chapter IX\], flat forms $`(\rho _i)`$ have a well defined trace on $`X`$, and on the measurable tangent bundle of $`X`$, essentially because of the rectifiability.555There is a technical point about orientation which we ignore here \[17, 3.26\]. Because $`f=\overline{f}|X`$ has a Lipschitz inverse, there exists a constant $`c>0`$ such that

$$(\rho _1\mathrm{}\rho _n)c>0$$
(6.2)
almost everywhere on $`X`$, where the Hodge star operator $``$ is determined by the chosen orientation on $`X`$.

Condition (6.2) was turned into a definition in . We say that $`X`$ admits local Cartan-Whitney presentations if for each point $`pX`$ one can find an $`n`$-tuple of flat 1-forms $`\rho =(\rho _1,\mathrm{},\rho _n)`$ defined in an $`𝐑^N`$-neighborhood of $`p`$ such that condition (6.2) is satisfied on $`X`$ near the point $`p`$.

###### Theorem 1

 Let $`X𝐑^N`$ be a metrically $`n`$-dimensional, linearly locally contractible $`n`$-manifold admitting local Cartan-Whitney presentations. Then $`X`$ is locally bi-Lipschitz equivalent to $`𝐑^n`$ outside a closed set of measure zero and of topological dimension at most $`n2`$.

To prove Theorem 1, fix a point $`pX`$, and let $`\rho =(\rho _1,\mathrm{},\rho _n)`$ be a Cartan-Whitney presentation near $`p`$. The requirement that $`\rho `$ be flat together with inequality (6.2) can be seen as a quasiregularity condition for forms.666In fact, (6.2) resembles a stronger, Lipschitz version of (2.1) studied in , , . We define a mapping

$$f(x)=_{[p,x]}\rho $$
(6.3)
for $`x`$ sufficiently near $`p`$, where $`[p,x]`$ is the line segment in $`𝐑^N`$ from $`p`$ to $`x`$, and claim that Reshetnyak’s program can be run under the stipulated conditions on $`X`$. In particular, we show that for a sufficiently small neighborhood $`U`$ of $`p`$ in $`X`$, the map $`f:U𝐑^n`$ given in (6.3) is a branched covering which is locally bi-Lipschitz outside its branch set $`B_f`$, which furthermore is of measure zero and of topological dimension at most $`n2`$. It is important to note that $`\rho `$ is not assumed to be closed, so that $`df\rho `$ in general.

In executing Reshetnyak’s proof, we use recent advances of differential analysis on nonsmooth spaces , , , as well as the theory developed simultaneously in . Incidentally, we avoid the use of the Harnack inequality for solutions, and therefore a deeper use of equation (2.2); this small improvement to Reshetnyak’s argument was found earlier in a different context in .

Theorem 1 provides bi-Lipschitz coordinates for $`X`$ only on a dense open set. In general, one cannot have more than that. The double suspension of a homology 3-sphere, as discussed in Section 4, can be mapped to the standard 5-sphere by a finite-to-one, piecewise linear sense-preserving map. By pulling back the standard coframe by such map, we obtain a global Cartan-Whitney presentation on a space that is not locally bi-Lipschitz equivalent to $`𝐑^5`$. Similar examples in dimension $`n=3`$ were constructed in , , by using Semmes’s spaces , . On the other hand, we have the following result:

###### Theorem 2

Let $`X𝐑^N`$ be a metrically $`2`$-dimensional, linearly locally contractible $`2`$-manifold admitting local Cartan-Whitney presentations. Then $`X`$ is locally bi-Lipschitz equivalent to $`𝐑^2`$.

Theorem 2 is an observation of M. Bonk and myself. We use Theorem 1 together with the observation that, in dimension $`n=2`$, the branch set consists of isolated points, which can be resolved. The resolution follows from the measurable Riemann mapping theorem together with the recent work by Bonk and Kleiner . While Theorem 2 presents a characterization of surfaces in Euclidean space that admit local bi-Lipschitz coordinates, we do not know whether the stipulation about the existence of local Cartan-Whitney presentations is really necessary (compare Problem 1 and the discussion preceding it).

For dimensions $`n3`$, it would be interesting to know when there is no branching in the map (6.3). In , we ask if this be the case when the flat forms $`(\rho _i)`$ of the Cartan-Whitney presentation belong to a Sobolev space $`H_{loc}^{1,2}`$ on $`X`$. The relevant example here is the map $`(r,\theta ,z)(r,2\theta ,z)`$, in the cylindrical coordinates of $`𝐑^n`$, which pulls back the standard coframe to a frame that lies in the Sobolev space $`H_{loc}^{1,2ϵ}`$ for each $`ϵ>0`$. Indeed, it was shown in  that in $`𝐑^n`$ every (Cartan-Whitney) pullback frame in $`H_{loc}^{1,2}`$ must come from a locally injective mapping.

## 7. Other recent progress and open problems

In his 1978 ICM address, Väisälä  asked whether the branch set of a $`C^1`$-smooth quasiregular mapping is empty if $`n3`$. It was known that $`C^{n/(n2)}`$-smooth quasiregular mappings have no branching when $`n3`$. The proof in \[Ri, p. 12\] of this fact uses quasiregularity in a rather minimal way. In this light, the following recent result may appear surprising:

###### Theorem 3

 For every $`ϵ>0`$ there exists a degree two $`C^{3ϵ}`$-smooth quasiregular mapping $`f:𝐒^3𝐒^3`$ with branch set homeomorphic to $`𝐒^1`$.

We are also able to improve the previous results as follows:

###### Theorem 4

 Given $`n3`$ and $`K1`$, there exist $`ϵ=ϵ(n,K)>0`$ and $`ϵ^{}=ϵ^{}(n,K)>0`$ such that the branch set of every $`K`$-quasiregular mapping in a domain in $`𝐑^n`$ has Hausdorff dimension at most $`nϵ`$, and that every $`C^{n/(n2)ϵ^{}}`$-smooth $`K`$-quasiregular mapping in a domain in $`𝐑^n`$ is a local homeomorphism.

The second assertion in Theorem 4 follows from the first, by way of Sard-type techniques. The first assertion was known earlier in a local form where $`ϵ>0`$ was dependent on the local degree . Our improvement uses  together with the work  by Rickman and Srebro.

The methods in  fall short in showing the sharpness of Theorem 4 in dimensions $`n4`$ in two technical aspects. First, we would need to construct a quasiconformal homeomorphism of $`𝐑^n`$ to itself that is uniformly expanding on a codimension two affine subspace; moreover, such a map needs to be smooth outside this subspace. In $`𝐑^3`$, it is easier to construct a mapping with expanding behavior on a line; moreover, every quasiconformal homeomorphism in dimension three can be smoothened (with bounds) outside a given closed set .

We finish with some open problems related to branching and quasiregular mappings. The problems are neither new nor due to the author.

###### Problem 2

What are the possible values for the topological dimension of the branch set of a quasiregular mapping?

By suspending a covering map $`H^3𝐒^3`$, where $`H^3`$ is as in Section 4, and using Edwards’s theorem, one finds that there exists a branched covering $`𝐒^5𝐒^5`$ that branches exactly on $`𝐒^1𝐒^5`$. It is not known whether there exists a quasiregular mapping $`𝐒^5𝐒^5`$ with similar branch set. If no such map existed, we would have an interesting implication to a seemingly unrelated parametrization problem; it would follow that no double suspension of a homology 3-sphere with nontrivial fundamental group admits a quasisymmetric homeomorphism onto the standard 5-sphere, cf. , \[16, Question 12\].

By work of Bonk and Kleiner , the bi-Lipschitz parametrization problem in dimension $`n=2`$ is equivalent to an analytic problem of characterizing, up to a bounded factor, the Jacobian determinants of quasiconformal mappings in $`𝐑^2`$. An affirmative answer to Problem 1 in Section 5 would give an affirmative answer to the following problem.

###### Problem 3

(Compare \[16, Question 2\]) Is every $`A_1`$-weight in $`𝐑^2`$ locally comparable to the Jacobian determinant of a quasiconformal mapping?

An $`A_1`$-weight is a nonnegative locally integrable function whose mean-value over each ball is comparable to its essential infimum over the ball. See , , ,  for further discussion of this and related problems.

###### Problem 4

\[16, Question 28\] Is there a branched covering $`f:𝐒^n𝐒^n`$, for some $`n3`$, such that for every pair of homeomorphisms $`\varphi ,\psi :𝐒^n𝐒^n`$, the mapping $`\varphi f\psi `$ fails to be quasiregular?

Branched coverings constructed by using the double suspension are obvious candidates for such mappings. In \[15, 9.1\], we give an example of a branched covering $`f:𝐒^3𝐒^3`$ such that for every homeomorphism $`\psi :𝐒^3𝐒^3`$, $`f\psi `$ fails to be quasiregular. The example is based on a geometric decomposition space arising from Bing’s double .

We close this article by commenting on the lack of direct proofs for some fundamental properties of quasiregular mappings related to branching. For example, it is known that for each $`n3`$ there exists $`K(n)>1`$ such that every $`K(n)`$-quasiregular mapping is a local homeomorphism , \[28, p. 232\]. All known proofs for this fact are indirect, exploiting the Liouville theorem, and in particular there is no numerical estimate for $`K(n)`$. It has been conjectured that the winding mapping $`(r,\theta ,z)(r,2\theta ,z)`$ is the extremal here (cf. Section 6). Thus, if one uses the inner dilatation $`K_I(f)`$ of a quasiregular mapping, then conjecturally $`K_I(f)<2`$ implies that $`B_f=\mathrm{}`$ for a quasiregular mapping $`f`$ in $`𝐑^n`$ for $`n3`$ \[29, p. 76\].

Ostensibly different, but obviously a related issue, arises in search of Bloch’s constant for quasiregular mappings. Namely, by exploiting normal families, Eremenko  recently proved that for given $`n3`$ and $`K1`$, there exists $`b_0=b_0(n,K)>0`$ such that every $`K`$-quasiregular mapping $`f:𝐑^n𝐒^n`$ has an inverse branch in some ball in $`𝐒^n`$ of radius $`b_0`$. No numerical estimate for $`b_0`$ is known. More generally, despite the deep results on value distribution of quasiregular mappings, uncovered by Rickman over the past quarter century, the affect of branching on value distribution is unknown, cf. \[29, p. 96\].