text:
# The neutron charge form factor and target analyzing powers from ³(He)⃗(e⃗,e'n)–scattering

## Abstract

The charge form factor of the neutron has been determined from asymmetries measured in quasi–elastic $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}(\stackrel{}{\mathrm{e}},\mathrm{e}^{}\mathrm{n})`$ at a momentum transfer of 0.67 $`(\mathrm{GeV}/\mathrm{c})^2`$. In addition, the target analyzing power, $`A_y^o`$, has been measured to study effects of final state interactions and meson exchange currents.

Introduction: The form factors of the nucleon are fundamental observables. Precise data allow for sensitive tests of the theory of the strong interaction — quantum chromodynamics (QCD) — in the non-perturbative regime. A complete study of the theoretical concepts requires measurements not only for the proton but also for the neutron. Accurate data at low momentum transfer are also required for the calculation of nuclear form factors.

Due to the lack of a free neutron target only neutrons bound in light nuclei can be studied. In this case, determinations of the charge, $`G_{en}`$, and magnetic, $`G_{mn}`$, form factor from elastic or quasi–elastic cross section data via the Rosenbluth technique do not lead to data with the desired precision. The subtraction of the proton contribution, theoretical corrections due to the unfolding of the nuclear structure and corrections to final state interaction (FSI) and meson–exchange currents (MEC) limit the accuracy to $``$30%.

Measurements of precise data of the neutron form factors became possible by means of alternative techniques exploiting polarized electron beams and polarized targets or recoil polarimeters. The technique to determine $`G_{en}`$ with a precision of $`<`$10% relies on asymmetry measurements in quasi–free $`(e,e^{}n)`$ coincidence experiments in which the asymmetry is given by the interference term and is proportional to $`G_{en}G_{mn}`$ in the Plane Wave Impulse Approximation (PWIA). The small contribution of $`G_{en}`$ is amplified by the large value of $`G_{mn}`$ and a measurement of the asymmetry allows for significant improvements in the precision . The continuous wave (cw) electron beams available today allowed for the determination of $`G_{mn}`$ with accuracies of $``$2% .

Because of its low binding energy, the deuteron is usually employed for studies of neutron properties. However, for polarization studies $`{}_{}{}^{3}\mathrm{He}`$ is particularly suitable due to the fact that for the major part of the ground state wave function the spins of the two protons are coupled antiparallel, so that spin dependent observables are dominated by the neutron . In addition, at least at low $`\mathrm{Q}^2`$, corrections due to nuclear structure effects, FSI, and MEC can be calculated using modern three-body calculations. These calculations allow for a quantitative description of the three-nucleon system with similar precision as for the deuteron .

The asymmetry in double polarization experiments is determined with

$$A(\theta ^{},\varphi ^{})=\frac{1}{P_eP_t}\frac{N^+N^{}}{N^++N^{}}$$
(1)
where $`\theta ^{},\varphi ^{}`$ are the polar and the azimuthal angle of the target spin direction with respect to the three momentum transfer $`\stackrel{}{\mathrm{q}}`$. The polarizations of beam and target are given by $`P_e`$ and $`P_t`$ and the normalized $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}(\stackrel{}{\mathrm{e}},\mathrm{e}^{}\mathrm{n})`$ events for positive (negative) electron helicity are $`N^+(N^{})`$. With the target spin orientation parallel and perpendicular to $`\stackrel{}{\mathrm{q}}`$ two independent asymmetries $`A_{}=A(0^o,0^o)`$ and $`A_{}=A(90^o,0^o)`$ can be measured. In PWIA $`G_{en}`$ can then be determined via

$$G_{en}^{PWIA}=\frac{b}{a}G_{mn}\frac{(P_eP_tV)_{}}{(P_eP_tV)_{}}\frac{A_{}}{A_{}}$$
(2)
with the kinematical factors $`a`$ and $`b`$ . The factor $`V`$ accounts for a possible dilution due to contributions with vanishing asymmetry. As $`P_e`$, $`P_t`$, and $`V`$ do not depend on the target spin orientation they cancel in principle in the determination of $`G_{en}^{PWIA}`$. In practice, $`A_{}`$ and $`A_{}`$ are measured in sequence, as such $`P_e`$ and/or $`P_t`$ may change during the two asymmetry measurements. It will be discussed below that such changes are measured and accounted for.

In order to study the FSI–corrections necessary for the determination of $`G_{en}`$ the target analyzing power $`A_y^o`$ provides an experimental quantity that is sensitive to these effects. For an unpolarized beam and the target spin aligned perpendicular to the scattering plane the target analyzing power can be measured with

$$A_y^o=\frac{1}{P_t}\frac{N^{}N^{}}{N^{}+N^{}}$$
(3)
where $`N^{},(N^{})`$ are the normalized $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}(\stackrel{}{\mathrm{e}},\mathrm{e}^{}\mathrm{N})`$ events for target spin aligned parallel (antiparallel) to the normal of the scattering plane. For coplanar scattering $`A_y^o`$ is identical to zero in PWIA due to the combination of time reversal invariance and hermiticity of the transition matrix . Thus, a non–zero value of $`A_y^o`$ signals FSI and MEC effects and its measurement provides a check of the calculation of these effects.

The present paper reports about a new determination of $`G_{en}`$ from measurements of $`A_{}`$ and $`A_{}`$ of $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}(\stackrel{}{\mathrm{e}},\mathrm{e}^{}\mathrm{n})`$–scattering at a four momentum transfer of $`\mathrm{Q}^2`$=0.67$`(\mathrm{GeV}/\mathrm{c})^2`$. The same kinematics is chosen as for the measurements by Rohe et al. . In addition, the same technique and almost the same apparatus is employed which allows to combine the data reducing the statistical error bar of $`G_{en}`$ by almost a factor of two.

Consequently, the improved precision requires a careful determination of FSI and MEC effects. Target analyzing powers $`A_y^o`$ have been measured at the same $`\mathrm{Q}^2`$ and at $`\mathrm{Q}^2`$=0.37$`(\mathrm{GeV}/\mathrm{c})^2`$ in order to properly determine FSI and MEC corrections of the combined result.

Experimental Setup: At the Mainz Microtron (MAMI)  longitudinally polarized electrons with a polarization of $``$0.8 were produced with a strained layer GaAsP crystal at a typical current of $`10\mu A`$ . The polarized cw electron beam was accelerated to an energy of 854.5 MeV and guided to the three–spectrometer hall . The $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}`$–target consisted of a spherical glass container with two cylindrical extentions sealed with oxygen–free 25$`\mu `$m Cu–windows. Coating the glass container with Cs led to relaxation times of about 80h. The Cu–windows were positioned outside of the acceptance of the spectrometer ($``$5 cm) and shielded with Pb–blocks to minimize background from beam–window interactions. The $`{}_{}{}^{3}\mathrm{He}`$–target was polarized via metastable optical pumping to a typical polarization of 0.5 and compressed to an operating pressure of 4 bar .

Spectrometer A with a solid angle of 28 msr and a momentum acceptance of 20% was used to detect the quasi–elastically scattered electrons at a scattering angle of 78.6o. The recoiling nucleons were detected in coincidence with an array of plastic scintillator bars placed at 32.2o, the direction of $`\stackrel{}{\mathrm{q}}`$ for the maximum of the quasi–elastic peak.

The hadron detector consisted of an array of four layers of five plastic scintillator bars with dimensions $`50\times 10\times 10`$ cm3 preceeded by two 1 cm thick $`\mathrm{\Delta }`$E detectors for particle identification. The detector was placed at a distance of 160 cm from the target, resulting in a solid angle of 100 msr. The neutron efficiency during the present experiment was determined to 18.3%. The entire detector was shielded with 10 cm Pb except for an opening towards the target were the Pb–shield was reduced to 2 cm.

The entire $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}`$–target was enclosed in a rectangular box of 2 mm thick $`\mu `$-metal and iron. The box served as an effective shield for the stray field of the magnetic spectrometers and provided a homogeneous magnetic guiding field of $`410^4`$ T produced by three independent pairs of coils. With additional correction coils a relative field gradient of less than $`510^4`$ cm-1 was achieved. The setup also allowed for an independent rotation of the target spin in any desired direction with an accuracy of 0.1o by remote control.

The product of target and beam polarization was monitored during the data taking via determination of the asymmetry for elastic $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}`$$`(\stackrel{}{e},e)`$–scattering for which the form factors, hence the asymmetries, are accurately known . The analysis of these data resulted in a polarization product of $`0.279\pm 0.010`$ for runs with $`A=|A_{}|`$ and $`0.282\pm 0.003`$ for $`A=|A_{}|`$. The different error bar results from the sensitivity of elastic scattering to the target spin direction.

In addition, the time dependence of the polarization of the target cell was continuously measured during the experiment by Nuclear Magnetic Resonance, while the absolute polarization was measured by the method of Adiabatic Fast Passage . The mean target polarization from these measurements was $`0.356\pm 0.015`$. ¿From the elastic scattering data and the target polarization measurements a beam polarization of $`P_e=0.788\pm 0.036`$ was extracted which agreed well with the determination with a Møller polarimeter ($`0.827\pm 0.017`$).

Determination of $`G_{en}`$: To determine $`G_{en}`$ the asymmetries $`A_{}`$ and $`A_{}`$ of $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}(\stackrel{}{\mathrm{e}},\mathrm{e}^{}\mathrm{n})`$ have been measured. The same kinematics was chosen as in  with the motivation to combine the two measurements hereby decreasing the statistical error bar of $`G_{en}`$.

In the analysis of the data the neutron is identified with a cut on the coincidence time and the absence of a hit in the $`\mathrm{\Delta }`$E–amplitude spectrum for two consecutive $`\mathrm{\Delta }`$E–detectors. Neutrons from (p,n) charge exchange in the Pb-shielding contribute in first order to the dilution factor $`V`$, but the effect cancels in the determination of $`G_{en}`$ through equation 2.

In order to minimize the dependence on the target polarization, data were accumulated alternatively for $`A_{}`$ and $`A_{}`$ at regular intervals by corresponding rotations of the target spin. The polarization ratio that enters in the determination of $`G_{en}`$ (see equation 2) was unity within 2.6%.

Experimental corrections have been determined via Monte Carlo simulation of the experiment based on PWIA with the momentum distribution compiled by Jans et al. . Accounting for energy loss via bremsstrahlung and for asymmetric angle and momentum acceptances of the spectrometer and the hadron detector can be reliably done . The dominant correction is due to the asymmetric acceptance of the electron spectrometer which leads to an angle shift of $`\stackrel{}{\mathrm{q}}`$. The resultant effect is an enhancement of the measured $`A_{}`$–value due to the contribution proportional to $`G_{mn}^2`$. Bremsstrahlung and missing energy lead to a similar effect. The total correction from these effects amounts to $`7.4\pm `$3.0%.

Finally, the value for the magnetic form factor required for the determination of $`G_{en}`$ is taken from the parameterization by Kubon et al.  with $`G_{mn}`$= (1.037$`\pm `$0.012)$`\mu _nG_D`$ where $`\mu _n`$ is the magnetic moment of the neutron in units of nuclear magnetons and $`G_D`$ the dipole form factor. The resulting experimental value is $`G_{en}^{PWIA}=0.0416\pm 0.0102_{stat}\pm 0.0024_{syst}`$.

This value is in good agreement with the value by Rohe et al. . A weighted average of the two values leads to $`G_{en}^{PWIA}=`$0.0468 $`\pm 0.0064_{stat}`$ $`\pm 0.0027_{syst}`$ which corresponds to a reduction of the statistical error bar by almost a factor of two.

Target analyzing power: The target analyzing power $`A_y^o`$ has been measured for $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}(\mathrm{e},\mathrm{e}^{}\mathrm{n})`$ and $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}(\mathrm{e},\mathrm{e}^{}\mathrm{p})`$ at $`\mathrm{Q}^2`$=0.67$`(\mathrm{GeV}/\mathrm{c})^2`$ (the kinematics of the $`G_{en}`$–measurement) and at 0.37$`(\mathrm{GeV}/\mathrm{c})^2`$. The measurement at 0.37$`(\mathrm{GeV}/\mathrm{c})^2`$ was performed by lowering the beam energy to 600 MeV as the geometrical constraints of the target shielding box and the hadron detector did not permit a change of the scattering or recoil angle. An unpolarized beam was used and the target spin was aligned perpendicular to the scattering plane and reversed every 2 minutes.

The analysis of the $`A_y^o`$–data is very similar to the one described above. Electrons are accepted for energy transfers $`\omega =225290`$ MeV (314–408 MeV) for the low (high) $`\mathrm{Q}^2`$–point. The hadron is identified with a cut on the coincidence time and the $`\mathrm{\Delta }`$E–amplitude spectrum.

Contrary to the determination of $`G_{en}^{PWIA}`$, dilution effects do not cancel for $`A_y^o`$ and have to be determined. The 2cm Pb–absorber of the hadron detector leads to misidentified proton/neutron events due to charge exchange scattering in the Pb–absorber. The 3 times larger $`ep`$ cross section and the 5 times larger efficiency of the hadron detector for protons leads to a dilution effect that is negligible for $`A_{y(e,e^{}p)}^o`$ but must be taken into account for $`A_{y(e,e^{}n)}^o`$.

The correction was measured by replacing the $`{}_{}{}^{3}\mathrm{He}`$–gas in the target with hydrogen and tagging the recoil protons with the elastically scattered electrons. The fraction of protons, misidentified as neutrons amounts to $`0.18\pm 0.01`$ ($`0.13\pm 0.01`$) for the low (high) $`\mathrm{Q}^2`$–point. An additional contribution results from the uncorrelated background in the coincidence time spectrum determined to 0.056(0.025) for the $`(e,e^{}n)`$–events.

The $`A_{y(e,e^{}n)}^o`$–values have been corrected according to

$$A_{y(e,e^{}n)}^o=\frac{(A_{ytotal}^oxA_{yback}^o)}{(1x)}$$
(4)
with $`x`$ the total fraction of background events, $`A_{yback}^o`$ its analyzing power and $`A_{ytotal}^o`$ the analyzing power of the total $`(e,e^{}n)`$–yields.

The corrected experimental results for $`A_y^o`$ are shown in table 1. Total error bars are given. The errors are dominated by statistics with a small contribution of systematic errors due to false asymmetries and polarization measurements.

For both $`(e,e^{}n)`$ and $`(e,e^{}p)`$ the agreement of the experimental values at $`\mathrm{Q}^2`$=0.37$`(\mathrm{GeV}/\mathrm{c})^2`$ with the result of a complete calculation by Golak et al.  is quite satisfactory. Neglecting the contribution of MEC in the calculation has little effect on $`A_y^o`$ for $`(e,e^{}n)`$. On the other hand a calculation for $`(e,e^{}n)`$ was also performed setting the proton form factors $`G_{ep}`$ and $`G_{mp}`$ to zero. As can be seen from table 1 the effect is quite dramatic suggesting that 98% of the FSI–effect measured with $`A_{y(e,e^{}n)}^o`$ results from a coupling of the virtual photon to the proton followed by a $`(p,n)`$ charge exchange reaction in the three–body–system.

A similar theoretical study is not possible at $`\mathrm{Q}^2`$=0.67$`(\mathrm{GeV}/\mathrm{c})^2`$ due to the non–relativistic nature of present day calculations and the fact that the transferred energy is well above the pion production threshold.

FSI corrections of $`G_{en}`$: For the same reason the FSI–effects in $`A_{}`$/$`A_{}`$ which are needed as corrections to determine $`G_{en}`$ can not be calculated at this $`\mathrm{Q}^2`$ using today’s non–relativistic Faddeev codes. However, we will discuss in the following that with the measurements of $`A_y^o`$ and the measurements by Carasco et al.  a reliable estimate of the effects can be made. In this approach, we first determine $`G_{en}^{PWIA}`$ which accounts for relativistic kinematics — the only significant relativistic effect to consider at this $`\mathrm{Q}^2`$  — and then apply the FSI corrections based on the acquired experimental information.

Two effects contribute to the FSI correction and have to be considered at first order. First, the photon couples to one of the protons followed by a charge exchange process in the three–body–system simulating an $`(e,e^{}n)`$–event. At $`\mathrm{Q}^2`$=0.37$`(\mathrm{GeV}/\mathrm{c})^2`$ the complete Faddeev calculation by Golak et al.  which successfully predicted $`A_y^o`$ predicts a total FSI effect for the ratio $`A_{}`$/$`A_{}`$ of 25%. The calculations also show that the charge exchange process which is responsible for 98% of $`A_y^o`$ amounts to 60% of the total FSI–effect in $`A_{}`$/$`A_{}`$.

The ratio of the elementary cross sections $`\sigma _{ep}/\sigma _{en}`$, which is a measure for the probability of the photon coupling to a proton or a neutron is similar at $`\mathrm{Q}^2`$ of 0.37$`(\mathrm{GeV}/\mathrm{c})^2`$and 0.67$`(\mathrm{GeV}/\mathrm{c})^2`$. We therefore assume that the charge exchange process also contributes with 60% to the FSI effect in $`A_{}`$/$`A_{}`$ at 0.67$`(\mathrm{GeV}/\mathrm{c})^2`$.

With the experimental knowledge of $`A_y^o`$ at both $`\mathrm{Q}^2`$–values the contribution of the charge exchange processes in $`A_{}`$/$`A_{}`$ can be determined with the ratio of the experimental $`A_y^o`$–values scaling the effect to 0.67$`(\mathrm{GeV}/\mathrm{c})^2`$. This results in a FSI correction of 3% in $`G_{en}`$.

Second, the photon couples to the neutron followed by a rescattering process in the three–body–continuum which may also lead to FSI–effects. The effect for this type of FSI in $`A_{}`$ and $`A_{}`$ of $`(e,e^{}p)`$ has been discussed in detail by Carasco et al. . The results of a calculation which treats only the FSI between the two (slow) spectator nucleons agree well with the experimental data. The same calculation has been used to compute $`A_{}`$ and $`A_{}`$ of $`(e,e^{}n)`$. Contrary to the significant FSI effect for $`A_{}`$ and $`A_{}`$ of $`(e,e^{}p)`$ observed in  the corresponding contribution is small for the asymmetries of $`(e,e^{}n)`$. The resulting FSI–effect of this process averaged over the accepted phase space is 0.4% in $`G_{en}`$.

Thus, we conclude that the total FSI–correction to $`G_{en}`$ at 0.67$`(\mathrm{GeV}/\mathrm{c})^2`$ is small (of the order of 3.4%) and dominated by charge–exchange processes. The correction is accounted for in the final result with a relative uncertainty of 50% added in quadrature to the quadratic sum of the experimental uncertainties of the combined $`G_{en}`$–result. The final value of $`G_{en}`$=$`0.0484\pm 0.0071`$ is shown in figure 1. This result is in excellent agreement with the $`G_{en}`$–values deduced from the quadrupole form factor of elastic $`ed`$–scattering .

Conclusions: In the present experiment, $`G_{en}`$ has been measured via the double polarization reaction $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}(\stackrel{}{\mathrm{e}},\mathrm{e}^{}\mathrm{n})`$. It has greatly improved the accuracy of our knowledge of $`G_{en}`$ from such measurements at $`\mathrm{Q}^2`$=0.67$`(\mathrm{GeV}/\mathrm{c})^2`$. The applied contribution from FSI is estimated as (3.4$`\pm `$1.7)% at this high $`\mathrm{Q}^2`$ which is considerably smaller than the statistical uncertainty. The good agreement of the final value of $`G_{en}`$=$`0.0484\pm 0.0071`$ with data from other double polarization experiments corrected for FSI is very satisfactory. The value for $`G_{en}`$ also agrees well with $`G_{en}`$–values extracted from the quadrupole form factor determined in elastic $`ed`$–scattering.

Acknowledgments: This work was supported by the Schweizerische Nationalfonds, Deutsche Forschungsgemeinschaft (SFB 443), the Polish Committee for Scientific Research, the Foundation for Polish Science and the NSF. The numerical calculations have been performed on the Cray SV1 of the NIC, Jülich, Germany.


format: text

text:
# When Crossings Count — Approximating the Minimum Spanning TreeA preliminary version of the paper appeared in the 16th ACM Symposium of Computational Geometry, 166–175, 2000.

## 1 Introduction

Given a set of lines in the plane a natural measure of distances between any two points is the number of lines one has to cross to reach from one point to the other. This is a discrete distance measure that can be used to approximate the Euclidean distance and other distance measures. However, since this measure is defined by an arrangement of lines it is not locally defined and is thus computationally cumbersome. Finding the minimum spanning tree (MST) of a set of points, so that the number of intersections between the tree and the given set of lines is minimized, quantify how the set of points interact with the set of lines; see Figure 1. In fact, when the set of lines is the set of all possible lines, then this MST is the standard Euclidean MST \[AF97\] (here one minimizes the average number of edges of the MST crossed when picking a random line). Such an MST is related to a spanning tree of low stabbing number (STLSN) \[Wel92, Aga91\]. While the spanning tree of low stabbing number guarantee that any line intersects at most $`O(\sqrt{n})`$ edges of the spanning tree, the MST guarantees that the overall number of intersections between the tree and a given set of lines is minimized. Thus, if we have the set of lines in advance then the MST will have overall less intersections than the STLSN. The spanning tree of low-stabbing number was used in several applications, see for example \[Aga91, MWW91\]. In particular, having such an MST enables one: (i) to answer half-plane range queries in an efficient manner using a near linear space \[GHS91\], (ii) bound the complexity of the faces of the arrangement of lines that contain the points \[HS01\], and (iii) traverse between the points in an efficient way, so that the number of updates needed is minimized. (Imagine traversing among the points and maintaining the set of half-planes that contain the current point. Each time one crosses a line an update operation is performed.)

Computing the MST for the general case of arcs can be done in $`O(n^2\mathrm{log}n)`$ time by performing wavefront propagation from each of the points (see Section 2). As for approximation algorithms, Har-Peled and Sharir \[HS01\] gave recently an approximation algorithm for the case of arcs, computing a Steiner tree in expected running time $`w=O(\lambda _{t+2}(n+𝒲_{opt})\mathrm{log}(n))`$, where $`t`$ is the maximum number of intersections between a pair of arcs, $`\lambda _{t+2}()`$ is the maximum length of a Davenport-Schinzel sequence of order $`t`$, and $`𝒲_{opt}`$ is the weight of the optimal Steiner tree.111It is easy to verify that if we have triangle inequality then the Steiner tree weight is at least half the weight of the MST. The algorithm outputs a tree of weight $`w`$ (and thus gives roughly $`O(\mathrm{log}n)`$-approximation).

In this paper, we present two results:

* A near linear time $`(1+\epsilon )`$-approximation algorithm to the minimum-spanning tree under the crossing metric in the planar case.
* We show to embed the crossing metric among hyperplanes into a Hamming distance in high dimensions. As a result, we show how one can apply known subquadratic approximation algorithms for problems involving point-sets and hyperplanes in high dimensions (MST, clustering, matching, etc).

The connection between the crossing metric, and points in high dimension follows by interpreting the input points as points in abstract VC-space \[PA95\] induced by the lines. Namely, we associate with each point in the plane, an $`n`$-dimensional binary vector, where $`i`$-th coordinate indicate on which side of the $`i`$-th line the point lies. In this way, we mapped our input points into points lying on the $`n`$-dimensional hypercube. The crossing metric is no more than the Hamming distance between the mapped points. We can now deploy the techniques of \[IM98\] to those mapped points, yielding an approximation algorithm for the MST problem. Bringing down the running time to be subquadratic requires some additional work.

Specifically, we show how to compute a mapping of the points into space of dimension $`O(\mathrm{log}^7n)`$; this embedding can be computed in $`\stackrel{~}{O}(n^{4/3})`$ time222Here and in the rest of this paper $`f(n)=\stackrel{~}{O}(g(n))`$ iff $`f(n)=g(n)(1/\epsilon )^{O(1)}\mathrm{log}^{O(1)}n`$, and $`f(n)=O_\epsilon (g(n))`$ iff $`f(n)=O(g(n)/\epsilon ^{O(1)})`$, for $`n`$ points, so that we get a $`(1+\epsilon )`$ gap property for a specified range of distances is preserved.

As a result, we can solve several approximation problems for this metric, among them is the MST problem. In fact, our near-linear approximate MST algorithm in the plane can be roughly viewed as an unraveling of the corresponding MST approximation algorithm in high dimensions. Similar bounds can be derived for $`d>2`$ dimensions. See Section 4 for details.

The paper is organized as follows: In Section 2, we describe how one can compute the exact MST using wavefront propagation. In Section 3, we present the planar $`(1+\epsilon )`$-approximation algorithm for the MST. Next, in Section 4, we describe the embedding into points in high dimension and demonstrate its usage for computing an approximate MST. Concluding remarks are given in Section 5.

## 2 Minimum Spanning Tree by Continuous Dijkstra

In this section, we present a simple algorithm for computing the crossing MST. It relies on a simple direct solution interpreted as a geometric algorithm. We also present a “weight sensitive” algorithm (Lemma 2.7) that computes portions of the MST in time proportional to its overall weight.

In the following, we assume that we are given a set $`L`$ of lines and a set $`P`$ of points in the plane. For simplicity, we assume $`|P|=|L|=n`$.

###### Definition 2.1

For a set $`L`$ of lines, the crossing metric is defined to be the minimum number of lines of $`L`$ that one has to cross as one moves between two prespecified points. Thus, for a pair of points $`p,q\mathrm{I}\mathrm{R}^2`$ the crossing distance between $`p`$ and $`q`$, denoted by $`𝒟_L(p,q)`$, is the number of lines of $`L`$ that intersects the segment $`pq`$. If $`L`$ is a set of arcs, a similar crossing metric is defined, although the “shortest path” in this case is no longer necessarily a straight segment.

###### Definition 2.2

For a set $`L`$ of lines, and a set $`P`$ of points in the plane, let $`𝒯_{opt}(P,L)`$ denote a minimum spanning tree of $`P`$ under the crossing metric induced by $`L`$, and let $`𝒲_{opt}(P,L)`$ denote the weight of $`𝒯_{opt}(P,L)`$.

Let $`𝒜=𝒜(L)`$ denote the planar arrangement induced by the lines of $`L`$. Let $`G_{adj}=G_{adj}(𝒜)`$ be the adjacency graph of $`𝒜`$; namely, each face of $`𝒜`$ is a vertex, and two vertices are connected if the two corresponding faces share an edge. Let $`V`$ be the set of vertices of $`G_{adj}`$ that corresponds to the faces of $`𝒜`$ that contains points of $`P`$. Clearly, the crossing MST of $`P`$ in $`𝒜`$, corresponds to the MST of $`V`$ in the graph $`G_{adj}`$ (here, each edge has associated weight $`1`$).

Computing the MST of $`V`$ in $`G_{adj}`$ can be done by performing a simultaneous flooding of $`G_{adj}`$ from the vertices of $`V`$. Indeed, we compute in the $`i`$-th iteration all the vertices of $`G_{adj}`$ that are in distance $`i`$ from any vertex of $`V`$. This can be easily done using a modified BFS. In the beginning, the flood front is made out of $`n`$ connected components. Every time two connected components of the flood front collide, we discovered a new edge of the MST. This edge connects the two vertices that induced the two parts of the wavefront that collided. This is a somewhat non-standard algorithm for computing the MST, but one can easily verify that it indeed computes the MST of $`V`$ in $`G_{adj}`$.

This flooding algorithm has a natural geometric interpretation: Let $`_{2i}`$ denote the set of all faces of $`𝒜`$ that are in (crossing) distance at most $`i`$ from any point of $`P`$. Clearly, $`_0`$ is the set of faces of $`𝒜`$ that contain points of $`P`$. The algorithm works in $`n/2`$ phases. We do a wavefront propagation in $`G_{adj}`$, starting from all the vertices that correspond to the marked faces (i.e., faces of $`𝒜`$ that contain points of $`P`$). In each iteration, we propagate the wavefront from the faces of $`_{2i2}`$ into the faces of $`_{2i}`$. It is easy to verify that a connected component of the flood corresponds to a connected component of the wavefront of $`_{2i}`$. (Note, that two faces of $`_{2i}`$ might be adjacent but belong to different wavefronts as the wavefronts did not cross the separating edge yet and thus were not merged into a single wavefront.) The connected components are maintained implicitly by a union-find data-structure. In particular, during the $`i`$-th iteration of the wavefront propagation in $`G_{adj}`$, when two different connected components of the wavefront collide, it corresponds to two points of $`P`$ with crossing distance equal to $`2i1`$ or $`2i`$ from each other.

In particular, if there is an edge of the MST of weight $`2i1`$ or $`2i`$ it would be discovered when the corresponding wavefronts collide. The $`i`$-th iteration of the wavefront propagation, corresponds to the detection of edges of weight $`2i1`$ and $`2i`$ in the MST. For the MST applications, we first handle all relevant edges of weight $`2i1`$, and later all such edges of weight $`2i`$. This requires a somewhat careful implementation, and we omit the the technical but straightforward details. See Figure 2.

Note, that the wavefront propagation can be done without constructing $`G_{adj}`$ in advance, and one can compute parts of $`G_{adj}`$ on the fly as needed (i.e., we need to compute only the parts of $`G_{adj}`$ that are covered by the wavefront, or are about to be covered). Of course, in the worst case, the whole graph $`G_{adj}`$ would be computed, which takes $`O(n^2\mathrm{log}n)`$ time (this corresponds to computing the whole arrangement $`𝒜(L)`$).

###### Lemma 2.3

Given a set $`L`$ of $`n`$ lines, and a set $`P`$ of $`n`$ points, a minimum spanning tree $`𝒯_{opt}(P,L)`$ of $`P`$ under the crossing metric $`𝒟_L`$ can be computed in $`O(n^2\mathrm{log}n)`$ time.

###### Remark 2.4

In the algorithm of Lemma 2.3 we did not use the fact that $`L`$ is a set of lines. The same algorithm will work for the case where $`L`$ is a set of arcs. Since we do not have the triangle inequality in this case, the edges of the MST are no longer line segments, but rather a Jordan arcs. (For example, imagine that the set $`L`$ is a single segment and we would like to connect two points that are separated by this segment. This can be done with no crossing by going “around” this segment.)

To be able to generate parts of $`G_{adj}`$ incrementally, as we perform the wavefront propagation, we need a way to compute the relevant portions of $`𝒜(L)`$ on the fly.

###### Theorem 2.5 (\[HS01\])

Let $`L`$ be a set of $`n`$ lines, as above, and $`P`$ a set of $`m`$ points in the plane. Then one can compute, in expected $`O(\left(n+w+m\right)\alpha (n)\mathrm{log}n)`$ time, a Steiner tree $`\widehat{}`$ of $`P`$, so that the expected weight of $`\widehat{}`$ is $`O((n+w)\alpha (n)\mathrm{log}n)`$, where $`w=𝒲_{opt}(P,L)`$ and $`\alpha (n)`$ is the inverse of the Ackermann function. Alternatively, one can compute the $`m`$ faces that contain the points of $`P`$ in the same time bound.

###### Lemma 2.6 (\[Aga91\])

There exists a Steiner tree $`^{}`$ of $`P`$, so that $`𝒲_{opt}(P,L)=O(n\sqrt{n})`$, and this is tight in the worst case (even for the case the arcs are lines).

In the worst case, Theorem 2.5 is inferior to implicit point-location data-structures \[AMS98\] (which can perform the implicit point-location needed in roughly $`O(n^{4/3})`$ time for $`m=n`$), as implied by Lemma 2.6 (as the weight of the MST is $`\mathrm{\Omega }(n^{3/2})`$ in the worst case, and this is the time to compute the relevant portions of the arrangement using the algorithm of Theorem 2.5). However, the running time of the algorithm of Theorem 2.5 is sensitive to the overall weight of the MST. This would be crucial for our algorithm.

###### Lemma 2.7

Given a set $`L`$ of $`n`$ lines, a set $`P`$ of $`n`$ points, and a parameter $`i`$, one can compute, in expected $`O(i(n+𝒲_{opt})\alpha ^2(n)\mathrm{log}n)`$ time, a minimum spanning forest of $`P`$ under the crossing metric $`𝒟_L`$, that connects all the points of $`P`$ in distance at most $`2i`$ from each other, where $`𝒲_{opt}=𝒲_{opt}(P,L)`$.

###### Proof.

The wavefront propagation on $`G_{adj}`$ can be done using an implicit representation of the arrangement of $`𝒜(L)`$. Namely, we compute the set $`_i`$ of faces of $`𝒜(L)`$ in distance $`i`$ from the points of $`P`$. Observe that the complexity of $`_i`$ is $`O((n+𝒲_{opt})i\alpha (n/i))`$. Indeed, the points of $`P`$ can be connected by an arc $`\gamma =𝒯_{opt}(P,L)`$ having $`O(𝒲_{opt})`$ intersections with the lines of $`L`$, and let $`𝒜^{}`$ be the arrangement resulting from $`𝒜`$ by creating a tiny gate for each intersection of $`\gamma `$ with the lines of $`L`$. The zone of $`\gamma `$ in $`𝒜(L)`$ corresponds to a single face $`F`$ of $`𝒜^{}`$, and the faces of $`_i`$ are contained in the set of faces in distance $`i`$ from $`F`$. By \[dBDS95\], the complexity of this region is $`O((n+𝒲_{opt})i\alpha (n/i))`$ (this is a bound on the complexity of all the vertices in distance $`i`$ from the face $`F`$.).

Clearly, the faces of $`_i`$ have a spanning tree of weight $`O((n+𝒲_{opt})i\alpha (n/i))`$, and so it can be computed in an online fashion in $`O((n+𝒲_{opt})i\alpha ^2(n)\mathrm{log}n)`$ expected time, by Theorem 2.5. ∎

## 3 Approximation Algorithm for the Planar Case

The algorithm is depicted in Figure 3, Figure 4 and Figure 5. We next describe the algorithm and its analysis in more detail.

Lemma 2.7 provides us with an algorithm for approximating the MST in roughly quadratic time in the worst case. To get a near linear running time, we simulate the Dijkstra algorithm by performing the wavefront propagation in an approximate fashion.

###### Definition 3.1

A metric $`𝒟^{}`$ $`\epsilon `$-approximates a metric $`𝒟`$, if for any $`p,q,r,sP`$ such that $`𝒟^{}(p,q)𝒟^{}(r,s)`$ then $`𝒟(p,q)(1+\epsilon )𝒟(r,s)`$.

###### Definition 3.2

For a set $`F`$ of segments in the plane, and a metric $`𝒟`$, let $`weight_𝒟(F)=_{eF}𝒟(e)`$ denote the total weight of $`F`$ under the metric $`D`$.

The proof of the following lemma is straightforward, and is included only for the sake of completeness.

###### Lemma 3.3

Let the metric $`𝒟^{}`$ be an $`\epsilon `$-approximation to the metric $`𝒟`$ over a point-set $`P`$. Let $`T^{}`$ be an MST of $`P`$ under $`𝒟^{}`$. Then, $`weight_𝒟(T^{})(1+\epsilon )weight_𝒟(T)`$, where $`T`$ is the MST of $`P`$ under $`𝒟`$, and $`weight(T)`$ is the total weight of the edges of $`T`$.

###### Proof.

Let $`e_1^{},\mathrm{},e_{n1}^{}`$ be the the edges of $`T^{}`$ sorted by their weight $`𝒟^{}(e_1^{})\mathrm{}𝒟^{}(e_{n1}^{})`$. Let $`T_0=T`$, and let $`T_i`$ be the tree resulting from removing the heaviest edge (according to $`𝒟^{}`$) from the cycle present in $`T_{i1}\left\{e_i^{}\right\}`$ (if $`e_i^{}`$ is already in $`T_{i1}`$ we do nothing). Let $`e_i`$ denote this removed edge. Clearly, $`𝒟^{}(e_i^{})𝒟^{}(e_i)`$ and, by definition, $`𝒟(e_i^{})(1+\epsilon )𝒟(e_i)`$. Namely, we replaced an edge $`e_i`$ by an edge $`e_i^{}`$ which is heavier by a factor of $`(1+\epsilon )`$. In the end of the process $`T_{n1}`$ is just $`T^{}`$, and $`weight_𝒟(T^{})_{i=1}^{n1}(1+\epsilon )weight_𝒟(e_i)(1+\epsilon )weight_𝒟(T)`$. ∎

Lemma 3.3 suggest that if we can find a computationally cheaper approximate metric than $`𝒟_L(,)`$, then we can use it to compute the MST. A natural way to do that, is to randomly sample a subset $`RL`$, and use $`𝒟_R(,)`$ as the approximate metric. However, it is easy to verify that $`𝒟_R`$ is an $`\epsilon `$-approximate metric to $`𝒟_L`$ only if $`L=R`$.

###### Definition 3.4

Let $`𝒟^{},𝒟`$ be two metrics, $`\epsilon >0`$, and $`l`$ be prescribed parameters. The metric $`𝒟^{}`$ is an $`(\epsilon ,l)`$-approximation to $`𝒟`$, if for any $`p,q,r,sP`$, such that (i) $`𝒟(p,q),𝒟(r,s)l`$, and (ii) $`𝒟^{}(p,q)𝒟^{}(r,s)`$, we have $`𝒟(p,q)(1+\epsilon )𝒟(r,s)`$.

Namely, $`𝒟^{}`$ $`\epsilon `$-approximates $`𝒟`$ for distances not smaller than $`l`$.

###### Definition 3.5

For $`l,\epsilon `$, let $`\nu (l,\epsilon )=\mathrm{max}(128c_{samp}\frac{\mathrm{log}n}{l\epsilon ^2},1)`$, where $`c_{samp}`$ is an appropriate constant. Let $`𝒮(L,l,\epsilon )`$ be a random subset of $`L`$ generated by picking independently each line of $`L`$ with probability $`\nu (l,\epsilon )`$.

Let $`\rho (l,\epsilon )=\nu (l,\epsilon )l=128c_{samp}\frac{\mathrm{log}n}{\epsilon ^2}`$. The value $`\rho (l,\epsilon )`$ is the expected crossing distance in $`𝒜(𝒮(L,l,\epsilon ))`$ between two points $`p,qP`$ such that $`𝒟_L(p,q)=l`$.

###### Lemma 3.6

Let $`L`$ be a set of $`n`$ lines in the plane, $`l`$ a positive integer number, $`\epsilon >0`$, and let $`R=𝒮(L,l,\epsilon )`$ be a random subset of $`L`$.

For any two points $`p,q`$ of distance $`𝒟_L(p,q)l`$ from each other we have

$$𝒟_L(p,q)\frac{n}{r(1\epsilon /4)}𝒟_R(p,q)(1+\epsilon )𝒟_L(p,q),$$

with probability $`1n^{c_0}`$.

Furthermore, $`𝒟_R(,)`$ is an $`(\epsilon ,l)`$-approximation to $`𝒟_L(,)`$ with high probability.

###### Proof.

Indeed, let $`X_{pq}=D_R(p,q)`$. We have,

$`\mu =E[X_{pq}]=𝒟_L(p,q)\nu (l,\epsilon )128𝒟_L(p,q)c_{samp}{\displaystyle \frac{\mathrm{log}n}{l\epsilon ^2}}{\displaystyle \frac{128c_{samp}\mathrm{log}n}{\epsilon ^2}}.`$
By Chernoff inequality \[MR95, MPS98\], we have that

$`P\left[\left|X_{pq}\mu \right|>{\displaystyle \frac{\epsilon }{4}}\mu \text{}\right]`$ $``$ $`2\left({\displaystyle \frac{e^{\epsilon /4}}{\left(1+\frac{\epsilon }{4}\right)^{1+\epsilon /4}}}\right)^\mu =2\mathrm{exp}\left(\mu \left({\displaystyle \frac{\epsilon }{4}}\left(1+{\displaystyle \frac{\epsilon }{4}}\right)\mathrm{log}\left(1+{\displaystyle \frac{\epsilon }{4}}\right)\right)\right)`$
$``$ $`2\mathrm{exp}\left(\mu \left({\displaystyle \frac{\epsilon }{4}}(1+{\displaystyle \frac{\epsilon }{4}})\left({\displaystyle \frac{\epsilon }{4}}{\displaystyle \frac{\epsilon ^2}{32}}\right)\right)\right)`$
$``$ $`2\mathrm{exp}\left(\mu {\displaystyle \frac{\epsilon ^2}{64}}\right)\mathrm{exp}\left({\displaystyle \frac{128c_{samp}\mathrm{log}n}{\epsilon ^2}}{\displaystyle \frac{\epsilon ^2}{64}}\right)n^{c_{samp}},`$
since $`\mathrm{log}(1+x)xx^2/2`$, for $`0x1`$. In particular, this implies that with high probability $`\mu (1\epsilon /4)X_{pq}\mu (1+\epsilon /4)`$. Namely, with high probability we have

$`𝒟_L(p,q)`$ $``$ $`{\displaystyle \frac{X_{pq}}{\nu (l,\epsilon )(1\epsilon /4)}}{\displaystyle \frac{\nu (l,\epsilon )(1+\epsilon /4)}{\nu (l,\epsilon )(1\epsilon /4)}}𝒟_L(p,q)={\displaystyle \frac{1+\epsilon /4}{1\epsilon /4}}𝒟_L(p,q)`$
$``$ $`(1+\epsilon )𝒟_L(p,q).`$
Consider now four points $`p,q,r,s`$, such that $`𝒟_L(p,q),𝒟_L(s,t)l`$ and $`𝒟_R(p,q)𝒟_R(r,s)`$. By the above discussion, we have with high probability

$$𝒟_L(p,q)\nu (l,\epsilon )(1\epsilon /4)𝒟_R(p,q)𝒟_R(r,s)(1+\epsilon )𝒟_L(r,s)\nu (l,\epsilon )(1\epsilon /4).$$

Namely, $`𝒟_L(p,q)(1+\epsilon )𝒟_L(r,s)`$. Namely, $`𝒟_R(,)`$ is an $`(\epsilon ,l)`$-approximation to $`𝒟_L(,)`$ with probability $`1\left(\genfrac{}{}{0pt}{}{n}{2}\right)n^{c_{samp}}`$. ∎

Lemma 3.6 and Lemma 3.3 suggest that we compute the MST by computing an appropriate random sample $`R`$ (by using a threshold $`l`$), and deploy the algorithms of Section 2 to compute the MST of $`P`$ in $`𝒜(R)`$. Such an MST would be an approximate MST. There are two main problems with this approach: (i) For short distances (i.e., $`l=1`$), just starting the wavefront propagation (i.e., Lemma 2.7) is prohibitively expensive (it roughly takes $`O(𝒲_{opt}(P,L))`$ time which might be $`\mathrm{\Omega }(n^{3/2})`$), (ii) For long distances (i.e., $`il`$), the wavefront propagation becomes, again, prohibitly expensive (i.e. $`\stackrel{~}{O}(ni)`$) by Lemma 2.7.

###### Corollary 3.7

Let $`U`$ be the total weight of all the edges of $`𝒯`$ having weigh less than $`\epsilon 𝒲_{opt}(P,L)/(10n)`$. Then $`U\epsilon 𝒲_{opt}(P,L)/10`$.

Lemma A.6 describes how we can approximate $`𝒲_{opt}(P,L)`$ to within a polylogarithmic factor using random sampling in near linear time. Since the algorithm of this lemma is very similar to the techniques used below, we defer its description to the appendix. Equipped with such approximation $`M`$, we know by Corollary 3.7 that we do not “care” about edges of the MST of length smaller than $`l_0=O(\epsilon M/(npolylog(n)))`$. In particular, we can generate a random sample $`R_0`$ which provides an $`(\epsilon ,l_0)`$-approximation to $`𝒟_L(,)`$. Thus, we can approximate the MST by computing the MST of $`𝒯_{opt}(P,R_0)`$.

This, however, does not address the second problem. Indeed, computing the MST of $`𝒯_{opt}(P,R_0)`$ might still be too expensive, as the following lemma testifies.

###### Lemma 3.8

Given a set $`L`$ of $`n`$ lines, a set $`P`$ of $`n`$ points, and parameters $`l,i,\epsilon ,U`$, such that $`l=\mathrm{\Omega }\left(𝒲_{opt}(P,L)/(nU)\right)`$ and let $`R=𝒮(L,l,\epsilon )`$ be a random sample of $`L`$. Then, one can compute, in expected $`\stackrel{~}{O}(iUn)`$ time, a minimum spanning forest of $`P`$ under the crossing metric $`𝒟_R`$, that connects all the points of $`P`$ in distance at most $`2i`$ from each other.

###### Proof.

Let $`X`$ denote the size of $`R`$. Clearly, The expected value of $`X`$ is

$$E[X]=n\nu (l,\epsilon )=128nc_{samp}\frac{\mathrm{log}n}{l\epsilon ^2}=O\left(\frac{Un^2\mathrm{log}n}{\epsilon ^2𝒲_{opt}(P,L)}\right),$$

by Definition 3.5. Let $`\gamma =𝒯_{opt}(\gamma ,L)`$. Let $`Y=weight(\gamma ,R)`$. Clearly,

$$E[Y]=weight(P,R)=𝒲_{opt}(P,L)\nu (l,\epsilon )=O\left(\frac{Un\mathrm{log}n}{\epsilon ^2}\right).$$

Namely, $`E[𝒲_{opt}(P,R)]E[Y]=O\left(\frac{Un\mathrm{log}n}{\epsilon ^2}\right)`$. The running time bound now follows immediately by applying the algorithm of Lemma 2.7 to $`P`$ and $`R`$. ∎

The algorithm of Lemma 3.8 first performs wavefront propagation for distances in $`𝒜(R)`$ which are smaller than $`\rho (l,\epsilon )`$. For such distances $`𝒜(R)`$ does not provide reliable estimate (i.e., ordering) of the crossing distances between points. However, once the distances propagated exceed $`\rho (l,\epsilon )`$, we know by Lemma 3.6 that the distances are now $`(\epsilon ,l)`$-approximated correctly. The main importance of the algorithm of Lemma 3.8 is that the algorithm has near linear running time for small values of $`U`$ and $`i`$.

Using Lemma 3.8 together with Corollary 3.7 implies that we can compute a spanning forest for the “short” edges of $`𝒯_{opt}(P,L)`$ in near linear time.

###### Lemma 3.9

Given a set $`P`$ of $`n`$ points in the plane, and a set $`L`$ of $`n`$ lines in the plane. One can compute a spanning forest $`F`$ of $`P`$, such that the weight of $`F`$ is $`\epsilon 𝒲_{opt}(P,L)/10`$. Furthermore, every pair of points of $`P`$ in distance $`\mathrm{\Omega }(𝒲_{opt}(P,L)\epsilon /(n\mathrm{log}^3n))`$ belong to the same connected components of $`F`$. The running time of this algorithm is $`\stackrel{~}{O}\left(n\right)`$.

###### Proof.

Using the algorithm of Lemma A.6, compute in $`\stackrel{~}{O}(n)`$ time, a number $`M`$ such that $`𝒲_{opt}(P,L)M=O(n\alpha (n)\mathrm{log}^2n+𝒲_{opt}(P,L)\alpha (n)\mathrm{log}n)`$. In particular, let

$$l_{short}=\frac{\epsilon M}{c_1n\mathrm{log}^3n}\frac{\epsilon }{40n}𝒲_{opt}(P,L),$$
(1)
for $`c_1`$ large enough. On the other hand, $`l_{short}=\mathrm{\Omega }(𝒲_{opt}(P,L)/(Un))`$, where $`U=O((\mathrm{log}^3n)/\epsilon )`$.

We now compute a spanning forest for $`P`$, using Lemma 3.8 with $`l_{short}`$ and $`U`$ as specified and $`i=2\rho (l,\epsilon )`$. The running time of this algorithm is

$$\stackrel{~}{O}\left(iUn\right)=\stackrel{~}{O}\left(\rho (l_{short},\epsilon )n\right)=\stackrel{~}{O}\left(\frac{\mathrm{log}n}{\epsilon ^2}n\right)=\stackrel{~}{O}\left(n\right).$$

Clearly, $`F`$ has at most $`n`$ edges, and all the points of $`P`$ in distance $`l_{short}`$ are in the same connected component of $`F`$ by Lemma 3.6.

Furthermore, for any edge $`pq`$ of $`F`$, we have that with high probability $`𝒟_L(p,q)2(1+\epsilon )l_{short}4l_{short}`$ by Lemma 3.6. In particular, $`weight(F,L)4nl_{short}(\epsilon /10)𝒲_{opt}(P,L)`$. ∎

Lemma 3.9 implies that we can compute a cheap spanning forest of $`P`$ in near linear time that “captures” all the light edges of the MST. Next, we can compute the rest of the edges of the MST using Lemma 3.8 repeatedly.

###### Lemma 3.10

Given a set $`P`$ of $`n`$ points in the plane, and a set $`L`$ of $`n`$ lines in the plane, a parameter $`\epsilon >0`$, and a spanning forest $`F`$ of $`P`$, such that every pair of points of $`P`$ in distance $`l`$ belong to the same connected components of $`F`$, where $`l=\mathrm{\Omega }(𝒲_{opt}(P,L)\epsilon /(n\mathrm{log}^3n))`$. Then, one can compute a spanning forest $`F^{}`$ of $`P`$ such that all the points of $`F`$ in distance $`2l`$ belong to the same connected component of $`F^{}`$. The forest $`F^{}`$ can be computed in $`\stackrel{~}{O}\left(n\right)`$ expected time.

###### Proof.

We use the same algorithm of Lemma 3.9, with the modification that when calling to the algorithm of Lemma 3.8, we pass on $`F`$, such that the algorithm ignore generated edges that belong to the same connected component of $`F`$. It is again clear, that only edges of length between $`l`$ and $`2(1+\epsilon )l`$ would be added to the spanning forest. The exact details of how to specify $`U`$ and $`i`$ are similar to Lemma 3.9, and are omitted. ∎

Our algorithm for computing the MST works by using Lemma 3.9. This results in a spanning forest $`F_0`$ of the points of $`P`$, and a value $`l_{short}`$ as specified by Equation (1). We now use Lemma 3.10 repeatedly $`O(\mathrm{log}n)`$ times, in the $`i`$-th iteration handling distances between $`2^{i1}l_{short}`$ to $`22^il_{short}(1+\epsilon )`$, for $`i=1,\mathrm{},O(\mathrm{log}n)`$), till we handle all distances $`n`$. Namely, in the $`i`$-th iteration, we compute a spanning forest $`F_i`$ of all points in distance $`2^il_{short}`$ from each other using Lemma 3.10 using $`F_{i1}`$ as our “starting” spanning forest.

Clearly, the expected running time of the resulting algorithm is $`\stackrel{~}{O}\left(n\right)`$. What is not clear, is that the resulting MST is indeed an $`\epsilon `$-approximate MST.

###### Lemma 3.11

With high probability, the tree $`T`$ computed by the above algorithm is an $`\epsilon `$-MST of $`P`$ in $`𝒜(L)`$.

###### Proof.

All the edges generated by the algorithm of Lemma 3.9, in the first stage of the algorithm, have total weight $`(\epsilon /20)𝒲_{opt}(P,L)`$ with high probability.

Let $`𝒯_{opt}(P,L)`$ be the optimal spanning tree. If $`T`$ is not an $`\epsilon `$-approximate MST, then $`weight_{𝒟_L}(T)>(1+\epsilon )weight_{𝒟_L}(𝒯_{opt})`$. In particular, there must be an edge of $`𝒯_{opt}`$ which its insertion into $`T`$ would results in substantially lightly spanning tree. Formally, for an edge $`e`$, let $`T(e)`$ be the tree resulting from $`T`$ by inserting $`e`$ into $`T`$, and removing from $`T`$ the heaviest (according to $`𝒟_L`$) edge on the new cycle that was created, and let $`\mathrm{out}(T,e)`$ denote this “ejected” edge.

Arguing as in the proof of Lemma 3.3, it must be that there exists an edge $`\varphi =pq`$ of $`𝒯_{opt}`$ such that

$$(1+\epsilon )𝒟_L(\varphi )<𝒟_L(\mathrm{out}(T,\varphi )),$$

and $`𝒟_L(\varphi )>𝒲_{opt}(P,L)/(20n)`$.

Let $`i`$ be the index such that $`2^{i1}l_{short}𝒟_L(\varphi )2^il_{short}`$. With high probability, we know that after the $`i`$-th iteration $`p`$ and $`q`$ are in the same connected component of $`F_i`$. Assume that $`p`$ and $`q`$ were not in the same connected component of $`F_{i1}`$ (the other case is easier and as such is omitted).

Let $`T^{\prime \prime }`$ be the spanning forest maintained by the algorithm just after $`p`$ and $`q`$ were present in the same connected component. With high probability, for any edge $`e^{\prime \prime }`$ of $`T^{\prime \prime }`$, we have $`𝒟_L(e^{\prime \prime })(1+\epsilon )𝒟_L(\varphi )`$, since the random sample $`R_i`$ we used in the $`i`$-iteration is $`(2^{i1}l_{short},\epsilon )`$-approximation to $`𝒟_L`$.

But then, it is not possible that the algorithm added $`\mathrm{out}(T,\varphi )`$ to the spanning tree $`T^{\prime \prime }`$, as all the edges on the cycle in $`T^{\prime \prime }\left\{\varphi \right\}`$ are lighter than $`(1+\epsilon )𝒟_L(\varphi )`$. A contradiction. ∎

We summarize our result:

###### Theorem 3.12

Given a set $`P`$ of $`n`$ points in the plane, $`L`$ a set of $`n`$ lines, and $`\epsilon >0`$ a parameter. Then one can compute a spanning tree $`T`$ of $`P`$, in $`\stackrel{~}{O}\left(n\right)`$ expected time, such that $`weight(T,L)(1+\epsilon )𝒲_{opt}(P,L)`$. The result is correct with high probability.

## 4 Approximation Algorithms for the Intersection Metric via Embeddings

Let $`P=\{p_1,\mathrm{},p_n\}`$ be a given set of $`n`$ points, and $`L=\{l_1,\mathrm{},l_m\}`$ be a set of $`m`$ lines, where $`m=n^{O(1)}`$. As mentioned earlier, the metric $`𝒟_L`$ is computationally cumbersome. One possible way to overcome this problem, is to embed this metric into a more convenient metric (while introducing a small distortion error).

In this section, we show a somewhat weaker result. We show how to embed the points of $`P`$ into $`O(\mathrm{log}^7n)`$-dimensional space in $`\stackrel{~}{O}(n+m+n^{2/3}m^{2/3})`$ time, so that a specific distance gap in the crossing metric, is mapped to a corresponding gap in the target space.

We first observe that the crossing distance between two points $`p`$ and $`q`$, can be computed by interpreting this distance as a Hamming distance on the hypercube in $`m`$ dimensions induced by the lines. Namely, each line $`l`$ contribute a coordinate — a point gets a ’1’ in this coordinate if it is on one side of $`l`$, and a ’0’ if it is on the other side of $`l`$. Formally, let $`l^+`$ denote the open half-plane defined by a line $`l`$ that contains the origin, and $`l^{}`$ denote the other open plane. For a point $`p\mathrm{I}\mathrm{R}^2`$, let $`\stackrel{}{v}_L(p)=(b_1,\mathrm{},b_m)`$ be a $`m`$-bit vector so that $`b_i=1`$ iff $`pl_i^+`$. It is easy to verify that $`𝒟_L(p,q)=d_H(\stackrel{}{v}_L(p),\stackrel{}{v}_L(q))`$, where $`d_H`$ is the Hamming distance.

###### Definition 4.1

Let $`RL`$, let $`f_R:\mathrm{I}\mathrm{R}^2𝖹𝖹`$ be the mapping that maps a point $`p`$ in the plane to its face ID in the arrangement $`𝒜(R)`$. Formally, we assign for each face in the arrangement $`𝒜(R)`$ a unique integer (say, and integer between $`1`$ and $`O(|R|^2)`$). The mapping $`f_R`$ maps a point $`p`$ in the plane to the integer identifying the face that contains $`p`$. (Note, that is does not uniquely define $`f_R()`$ as we did not specify how we assign the IDs to the faces.)

For a set $`=(R_1,\mathrm{},R_\mu )`$ of subsets of $`L`$, let $`f_{}:\mathrm{I}\mathrm{R}^2𝖹𝖹^\mu `$ be the mapping $`f_{}(p)=(f_{R_1}(p),f_{R_2}(p),\mathrm{},f_{R_\mu }(p))`$. For two points $`p,q\mathrm{I}\mathrm{R}^2`$, let $`d_H(f_{}(p),f_{}(q))`$ be the Hamming distance between $`f(p)`$ and $`f_{}(q)`$. Namely, this is the number of coordinates, where the two vectors $`f_{}(p)`$ and $`f_{}(q)`$ disagree.

One can view $`f_{}`$ as an embedding of the crossing metric $`𝒟_L`$ to the Hamming space $`𝖹𝖹^\mu `$.

###### Lemma 4.2

Given a set $`P`$ of $`n`$ points in the plane, a set $`L`$ of lines in the plane, a parameter $`\epsilon >0`$ and a parameter $`r`$. One can compute a set $``$ of $`\mu `$ subsets of $`L`$, such that for the embedding $`f_{}:\mathrm{I}\mathrm{R}^2𝖹𝖹^\mu `$, we have that, with high probability, for any $`p,qP`$ it holds:

* If $`𝒟_L(p,q)r`$, then $`d_H(f(p),f(q))M`$,
* If $`𝒟_L(p,q)(1+\epsilon )r`$ then $`d_H(f(p),f(q))(1+\epsilon )(1a/\mathrm{log}n)M`$,

where $`M`$ and $`a`$ are appropriate constants and $`\mu =O(\mathrm{log}^4n)`$.

###### Proof.

For sake of simplicity of exposition, we assume that $`m/r\mathrm{log}n`$, where $`m=|L|`$. If this is not correct, we can add “fictitious” lines to $`L`$ that have all the points of $`P`$ on one side of them. If we pick such a line to a set of $``$, we can ignore it when we compute the face IDs.

For a parameter $`\alpha `$ to be specified shortly, let $`k=\alpha m/r`$, $`R`$ be a sample of $`k`$ lines out of $`L`$ (performed with replacement), and let $`p,q`$ be two points of $`P`$. Let $`\rho =𝒟_L(p,q)/n`$. The probability that $`p,q`$ will be in two different faces of $`𝒜(R)`$ is

$$U(\rho )=1(1\rho )^k,$$

as this is the probability that not all the lines will miss the segment connecting $`p`$ and $`q`$.

Our target is to approximate the value of $`U(\rho )`$ so we could decide whether $`p,q`$ are close or far. Indeed, if $`U(\rho )U((1+\epsilon )r/m)`$ then $`𝒟_L(p,q)(1+\epsilon )r`$, and if $`U(\rho )U(r/m)`$ then $`𝒟_L(p,q)r`$.

To do so, we generate a set of subsets $`=(R_1,\mathrm{},R_\mu )`$, by random sampling as described above, where $`\mu `$ would be specified shortly. Now we consider the quality of the distance approximation provided by the embedding333A similar analysis (in the context of Hamming spaces) appeared already in \[Ind00\]; in our case, however, we have to put more care into the analysis, since we want $`\epsilon `$ and $`\epsilon ^{}`$ to be very close.. Let $`X(p,q)`$ denote the random variable which is the number of arrangements of $`𝒜(R_1),\mathrm{},𝒜(R_\mu )`$ that have $`p,q`$ in different faces. Note, that $`X(p,q)`$ is equal to the Hamming distance between $`f_{}(p)`$ and $`f_{}(q)`$, and it thus the distance between the images of $`p`$ and $`q`$ in the new space. Clearly, as $`\mu `$ tends to infinity, $`X(p,q)/\mu `$ tends to $`U(\rho )`$. Using Chernoff inequality, we can quantify the quality of approximation provided by $`\mu `$. Specifically, let $`z=U(r/m)`$ and $`Z=U((1+\epsilon )r/m)`$; in the following we will make sure that $`Z<1/2`$. Then, from the Chernoff bound \[MR95, MPS98\] it follows that for any $`\alpha >0`$ if $`\mu =C\frac{\mathrm{log}n}{z\alpha ^2}`$ for some constant $`C`$, then with high probability:

* if $`𝒟_L(p,q)r`$ then $`X(p,q)/\mu z(1+\alpha )`$
* if $`𝒟_L(p,q)r(1+\epsilon )`$ then $`X(p,q)/\mu Z(1\alpha )`$

Therefore, the mapping $`f_{}`$ converts the distance gap $`r:(1+\epsilon )r`$ into the gap $`z(1+\alpha )\mu :Z(1\alpha )\mu `$. We next fine tune $`k`$ (the size of each sample) so that the resulting gap will be as large as possible. (Intuitively, the larger the target gap is, the easier it is to detect it in later stages.) Therefore, in the following we focus on finding $`k`$ such that the ratio

$$\mathrm{\Delta }=\frac{Z(1\alpha )\mu }{z(1+\alpha )\mu }$$

is as large as possible. To this end, we observe that

$`z`$ $`=`$ $`U\left({\displaystyle \frac{r}{m}}\right)=1\left(1{\displaystyle \frac{r}{m}}\right)^k1e^{rk/m}\left(1{\displaystyle \frac{(rk/m)^2}{k}}\right)=1e^\alpha \left(1{\displaystyle \frac{\alpha ^2}{k}}\right)`$
$``$ $`1e^\alpha \left(1\alpha ^2\right)\alpha ^2+(1e^\alpha )\left(1\alpha ^2\right)\alpha ^2+\alpha \left(1\alpha ^2\right)\alpha (1+\alpha )`$
since $`\left(1{\displaystyle \frac{t}{n}}\right)^ne^t\left(1{\displaystyle \frac{t^2}{n}}\right)`$ \[MR95\], $`k=\frac{\alpha m}{r}`$, and $`x1e^x`$. Furthermore,

$`Z`$ $`=`$ $`U((1+\epsilon )r/m)=1(1(1+\epsilon )r/m)^k1e^{(1+\epsilon )rk/m}=1e^{(1+\epsilon )\alpha }`$
$``$ $`(1+\epsilon )\alpha ((1+\epsilon )\alpha )^2(1+\epsilon )\alpha (1(1+\epsilon )\alpha )`$
since $`(1t/n)^ne^t`$ \[MR95\] and $`1e^xxx^2/2xx^2`$.

Therefore

$`{\displaystyle \frac{Z}{z}}`$ $``$ $`{\displaystyle \frac{(1+\epsilon )\alpha (1(1+\epsilon )\alpha )}{\alpha (1+\alpha )}}(1+\epsilon )(1(1+\epsilon )\alpha )(1\alpha )(1+\epsilon )(1(2+\epsilon )\alpha ),`$
since $`1/(1+x)(1x)`$. Thus, if we set $`\alpha `$ to be $`1/\mathrm{log}n`$, then the distance gap becomes (at least)

$$\mathrm{\Delta }=\frac{Z(1\alpha )\mu }{z(1+\alpha )\mu }(1+\epsilon )(1(2+\epsilon )\alpha )(1\alpha )^2(1+\epsilon )\left(1\frac{a}{\mathrm{log}n}\right),$$

where $`a`$ is an appropriate constant. Also, note that the resulting value of $`z`$ is

$$z=1(1r/m)^k1e^{p_0k}=1e^\alpha \alpha \alpha ^2/2=\mathrm{\Omega }(1/\mathrm{log}n)$$

and $`\mu =(C\mathrm{log}n)/(z\alpha ^2)=C\mathrm{log}^2n/\alpha ^2=O(\mathrm{log}^4n)`$. Finally, since $`m/r\mathrm{log}n`$, we have that $`k=\alpha (m/r)=(1/\mathrm{log}n)(m/r)1`$ (i.e., the sample size $`k`$ is at least $`1`$). ∎

###### Lemma 4.3

Given a set $`P`$ of $`n`$ points, and a set $`L`$ of $`m`$ lines, one can compute the function $`f_{}()`$, of Lemma 4.2, for all the points of $`P`$ in $`\stackrel{~}{O}((m^{2/3}n^{2/3}+m+n))`$ expected time.

###### Proof.

We have to compute for each point of $`P`$ the face that contains it in each of the arrangements $`𝒜(R_1),\mathrm{},𝒜(R_\mu )`$, where $`\mu =O(\mathrm{log}^4n)`$. Or alternatively, compute all the faces of $`𝒜(R_1),\mathrm{},𝒜(R_\mu )`$ that contains points of $`P`$. For a single arrangement $`A_i`$ this can be done in $`O(m^{2/3}n^{2/3}\mathrm{log}^{2/3}(m/\sqrt{n})+(m+n)\mathrm{log}m)`$ expected time \[AMS98\]. Since there are $`\mu `$ coordinates (i.e., arrangements), the result follows. ∎

Thus, we showed how to embed $`𝒟_L`$ into $`\mu `$-dimensional Hamming space $`\mathrm{\Sigma }^\mu `$ in $`\stackrel{~}{O}(n+m+n^{2/3}m^{2/3})`$ time, mapping a $`(1+\epsilon )`$ gap between close and far points into a gap of size $`(1+\epsilon )(1O(1)/\mathrm{log}n)`$, where $`\mu =O(\mathrm{log}^4n)`$ and $`\mathrm{\Sigma }𝖹𝖹`$ is the set of face labels we use (i.e., $`|\mathrm{\Sigma }|=O(m^2)`$. By using standard embedding techniques (e.g. see \[KOR00\]) we can embed the Hamming space $`\mathrm{\Sigma }^\mu `$ into $`\{0,1\}^D`$ with $`D=O(\mu \mathrm{log}|\mathrm{\Sigma }|\mathrm{log}^2n)=O(\mathrm{log}^6n\mathrm{log}m)`$, preserving the gap up to another factor $`(1O(1)/\mathrm{log}n)`$. This gives an embedding of $`𝒟_L`$ into $`D=O(\mu \mathrm{log}m\mathrm{log}^2n)`$-dimensional binary Hamming cube, with error $`(1O(1)/\mathrm{log}n)`$. Thus it is sufficient for us to maintain $`c`$-nearest neighbor in $`\{0,1\}^D`$ where $`c=(1+\epsilon )(1O(1)/\mathrm{log}n)`$, which takes $`\stackrel{~}{O}(n^{1/c})=\stackrel{~}{O}(n^{1/(1+\epsilon /2)})`$ time per operation \[IM98\].

We conclude:

###### Theorem 4.4

By performing a $`\stackrel{~}{O}(n+m+n^{2/3}m^{2/3})`$-time preprocessing, one can reduce the problem of maintaining dynamic $`(1+\epsilon )`$-approximate nearest neighbor for any $`n`$-point crossing metric over $`m`$ lines, to the problem of maintaining dynamic $`(1+\epsilon )(1O(1)/\mathrm{log}n)`$-approximate nearest neighbor in Hamming space with $`O(\mathrm{log}^7n)`$ dimensions (assuming $`m=n^{O(1)}`$). The latter can be solved in $`\stackrel{~}{O}(n^{1/(1+\epsilon /2)})`$ time per operation.

### 4.1 Embedding of the Crossing Metric over $`\mathrm{I}\mathrm{R}^d`$

In this Section, we extend the methods from the previous section to the crossing metric defined by $`(d1)`$-dimensional hyperplanes in $`\mathrm{I}\mathrm{R}^d`$, for any fixed $`d2`$. To this end, it is sufficient to design an efficient procedure, which given a set of $`n`$ points $`P=p_1,\mathrm{},p_n`$ and a set of $`m`$ hyperplanes $`=\{H_1,\mathrm{},H_m\}`$, assigns a symbol $`\sigma _i\mathrm{\Sigma }𝖹𝖹`$ to each $`p_i`$ in such a way that $`\sigma _i\sigma _j`$ iff there exists $`H_k`$ which separates $`p_i`$ from $`p_j`$. Unfortunately, the idea from the previous section does not give subquadratic time algorithm for $`d>2`$, since even in $`d=3`$ the complexity of $`n`$ cells in an arrangement formed by $`n`$ planes could be $`\mathrm{\Omega }(n^2)`$. Fortunately, for our purpose, we do not need to compute the actual cells containing $`p_i`$s. Rather, it is just sufficient to find the labels for those cells, or more specifically, a function $`h:P\mathrm{\Sigma }`$ such that $`h(p)=h(q)`$ iff $`p`$ and $`q`$ belong to the same arrangement cell.

Abusing notations, we denote by $`H_k(p)`$ the function returning $`1`$ if $`p`$ lies on one side of $`H_k`$ and zero otherwise. We use the following hashing function

$$h(x)=\left(\underset{i}{}a_iH_i(x)\right),$$

where $`a_1\mathrm{}a_m`$ are independent and identically distributed random variables with uniform distribution over $`\{0,\mathrm{},n^c\}`$, where $`c`$ is a constant to be specified shortly. Note, that if $`p,q\mathrm{I}\mathrm{R}^d`$ lie in two different full-dimensional faces of $`𝒜()`$, then, as noted above, there must be a hyperplane $`H_k`$, so that $`H_k(p)H_k(q)`$, and say that $`H_k(p)=1`$. That is, $`h(p)=h^{}(p)+a_k`$ and $`h(q)=h^{}(q)`$, where $`h^{}(x)=_{ik}a_iH_i(x)`$. Since the $`a_i`$ were picked independently, it follows that $`h(p)=h(q)`$ only if $`h^{}(p)h^{}(q)=a_k`$. But the probability of that to happen is $`1/n^c`$. We conclude, that the probability of two points belonging to two different faces to be mapped to the same value by $`h()`$ is $`1/n^c`$. Thus, since we have $`O(n^2)`$ pairs of points to consider in our algorithm, it follows that the probability of the hashing to fail is $`n^{2c}`$ which can be made to be arbitrarily small by picking $`c`$ to be large enough.

Namely, we associate a weight $`a_i`$ with each half-space induced by a hyperplane $`H_i`$. For each point $`p_j`$, we compute the total weight of all the half-spaces that contain it, and all the points having the same total weight are associated with the same label. Computing the weight of a point $`p_j`$ falls into the class of problems known as intersection-searching \[Aga97\]. In particular, one can construct a data-structure in $`O(m^{1+\delta })`$ time, so that one can answer intersection-searching queries in $`O((n/m^{1/d})\mathrm{log}^{d+1}n)`$ time, where $`\delta >0`$ is arbitrarily small constant. As the algorithm needs to perform a linear number of such queries, we set $`m=n^{2d/(d+1)}`$. Thus, the algorithm computes the required labels in $`O(n^{2d/(d+1)+\delta })`$ time. We conclude:

###### Theorem 4.5

By performing a $`O(n^{2d/(d+1)+\delta })`$-time preprocessing, where $`\delta >0`$ is arbitrary constant, one can reduce the problem of maintaining dynamic $`(1+\epsilon )`$-approximate nearest neighbor for any $`n`$-point crossing metric over $`n`$ hyperplanes in $`\mathrm{I}\mathrm{R}^d`$, to the problem of maintaining dynamic $`(1+\epsilon )(1O(1)/\mathrm{log}n)`$-approximate nearest neighbor in Hamming space with $`O(\mathrm{log}^7n)`$ dimensions.

###### Remark 4.6

Note, that the constants in the bounds of Theorem 4.5 depend exponentially (or worse) on the dimension $`d`$.

###### Remark 4.7

As indicated in the introduction, having such a embedding, enable one to use a large collection of subquadratic approximation algorithms for the intersection metric, including dynamic amortized $`\stackrel{~}{O}(n^{4/3}+n^{1+1/c})`$-time (for $`d=2`$) $`c`$-approximation algorithms for bichromatic closest pair \[Epp95\] and $`\stackrel{~}{O}(n^{4/3}+n^{1+1/c})`$-time algorithms for: $`c`$-approximate diameter and discrete minimum enclosing ball \[GIV01\], $`O(c)`$-approximate facility location and bottleneck matching \[GIV01\]. Similar (i.e., subquadratic time) results hold for any $`d>2`$.

### 4.2 Computing an MST Using the Embedding

We next describe how to use the embedding described in the previous two sections, for getting an $`(1+\epsilon )`$-approximation algorithm for the MST under crossing metric. Note that everything described in this section is well known \[IM98\], and we provide it only for the sake of completeness. Also, the resulting algorithm is slower in the planar case than the algorithm of Section 3.

Computing the minimum spanning tree under the intersection metric, using the Kruskal’s algorithm, boils down to maintaining the bichromatic nearest-neighbor pair (under the intersection metric) between two sets $`P_1,P_2P`$, under insertions and deletions. A consequence of Eppstein result \[Epp95\] is the following:

###### Theorem 4.8 (\[Epp95\])

Given a dynamic data-structure for nearest-neighbor queries, where each insertion / deletion / query operation takes $`T(n)`$ time, then one can compute the MST in $`O(nT(n)\mathrm{log}^2n)`$ time.

It is easy to verify that if we get a $`(1+\epsilon )`$-approximation to the MST if we use an $`(1+\epsilon )`$-approximate dynamic nearest-neighbor data-structure (Eppstein, personal communication, 1999).

Namely, we need a data-structure that support dynamic approximation nearest-neighbor queries. After applying the embedding described above, we use the $`\epsilon ^{}`$-PLEB data-structure of \[IM98\] to maintain a $`(1+\epsilon ^{})`$-approximate nearest neighbor in the embedded space. Specifically, we construct an $`\epsilon `$-PLEB in the embedded points. In this way, we obtain an $`\epsilon `$-PLEB for our original points (i.e., we embedded a gap to a gap, so that a close point in the embedded space, corresponds to a close point in the crossing metric) data-structure that for a query $`p`$ return us a point of $`qP`$ so that $`𝒟_L(p,q)(1+\epsilon )r`$, if there exits a point $`q^{}P`$ so that $`𝒟_L(p,q^{})r`$.

Thus, by constructing $`\mathrm{log}_{1+\epsilon }n`$ such data-structures, we can use binary search on those data-structures to find and $`(1+\epsilon )`$-approximate nearest neighbor to a query point. Namely, this data-structure can be used to answer approximate nearest neighbor queries for the intersection metric. For the whole scheme to work, we need those data-structures to be dynamic; i.e., support insertions and deletions of points. Fortunately, the only part of the algorithm that needs to be dynamic is the second stage that uses the data-structure of \[IM98\] which is dynamic.

We conclude:

###### Theorem 4.9

Given a set $`P`$ of $`n`$ points in the plane, and a set $`L`$ of $`n`$ lines, one can compute in $`\stackrel{~}{O}\left(n^{4/3}+n^{1+1/(1+\epsilon )}\right)`$ time, a spanning tree of $`P`$ of weight $`(1+\epsilon )𝒲_{opt}(P,L)`$. The result returned by the algorithm is correct with high probability. For $`d>2`$ dimensions, such an MST can be approximated in $`\stackrel{~}{O}\left(n^{2d/(d+1)+\delta }+n^{1+1/(1+\epsilon )}\right)`$ time, where $`\delta >0`$ is an arbitrary constant.

## 5 Conclusions

We presented the first $`(1+\epsilon )`$-algorithm for approximating the minimum spanning tree under the crossing metric in the plane. We also presented a subquadratic time approximation algorithms for a variety of other problems, obtained by embedding the crossing metric into higher dimensional space. The techniques used in our paper seems to be new to low-dimension computational geometry, and we believe that they might be useful for other problems in computational geometry.

There are several interesting open problems for further research:

* Can the result be extended to other cases: segments or arcs instead of lines?
* Can a similar approximation algorithm be found for the case of minimum weight triangulation under the crossing metric?

### Acknowledgments

The authors wish to thank Pankaj Agarwal, Boris Aronov and Micha Sharir for helpful discussions concerning the problems studied in this paper and related problems.

## Appendix A A Rough Approximation to the Weight of the MST in Near Linear Time

In this appendix, we show how to approximate the weight of the minimum spanning tree up to roughly a factor of $`O(\alpha (n)\mathrm{log}n)`$ if its weight is at least linear. In Section 3, we presented a near linear time algorithm for $`(1+\epsilon )`$-approximation for the minimum spanning tree, that relies on this approximation algorithm.

Underlining the approximation algorithm, is the observation that an MST for a random sample of the lines of $`L`$ provides a rough approximation to the weight of the MST of $`L`$. If the weight of the MST of the sample is near linear, we can approximate it up to a $`O(\alpha (n)\mathrm{log}n)`$, using the following algorithm.

###### Lemma A.1

Given a set $`R`$ of $`r`$ lines, $`P`$ a set of $`n`$ points, and $`W`$ a prescribed parameter, one can decide whether $`𝒲_{opt}(P,R)`$ is large; namely, $`𝒲_{opt}(P,R)=\mathrm{\Omega }((r+n+W)\alpha (n)\mathrm{log}n)`$. The algorithm takes $`O((r+n+W)\alpha (n)\mathrm{log}^2n)`$ expected time. Furthermore, if $`𝒲_{opt}(P,R)W`$, the algorithm will report that its weight is large with probability at most $`n^c`$, where $`c`$ is an appropriate constant.

###### Proof.

Use the algorithm of Theorem 2.5 and execute it $`O(\mathrm{log}n)`$ times on $`P`$ and $`R`$. If the running time of the $`i`$-th execution of the algorithm exceeds $`\mathrm{\Omega }((r+n+W)\alpha (n)\mathrm{log}n)`$ abort it, and move on to the next execution. If $`𝒲_{opt}(P,R)W`$, then the algorithm of \[HS01\] provides a spanning tree of expected weight $`O((r+n+W)\alpha (n)\mathrm{log}n)`$ with the same bound on the expected running time. Thus, if in $`O(\mathrm{log}n)`$ executions the algorithm returns always that $`𝒲_{opt}`$ is large, we can conclude that with probability $`1n^c`$ the weight of $`𝒲_{opt}(P,R)`$ is not $`W`$. ∎

Lemma A.1 shows that we can approximate the weight of the MST in near linear time if its weight is near linear. However, if it is heavier, we will use random sampling to keep the running time under control.

Let $`RL`$ be a random sample of lines out of $`L`$, where each line is picked independently with probability $`r/n`$. Clearly, the probability of an intersection point $`u`$ (between a connected set $`\gamma `$ and a line of $`L`$), to be present in $`𝒜(R)`$ is $`r/n`$ (this is the probability that the line of $`L`$ passing through $`u`$ will be chosen to be in the random sample).

###### Definition A.2

For a curve $`\gamma `$, and a set of lines $`L`$, let $`weight(\gamma ,L)`$ denote the weight of $`\gamma `$ in the arrangement $`𝒜(L)`$. This is the number of intersections of $`\gamma `$ with the lines of $`L`$.

###### Lemma A.3

Let $`R`$ be a sample of lines of $`L`$ (chosen as described above), then with high probability:

$$𝒲_{opt}(P,L)\frac{n}{r}\left(c_0n\mathrm{log}n+2𝒲_{opt}(P,R)\right),$$

and with probability $`0.9`$ we have $`\frac{n}{r}\frac{𝒲_{opt}(P,R)}{10}𝒲_{opt}(P,L)`$, where $`c_0`$ is an appropriately large constant.

###### Proof.

Let $`𝒯_{opt}^L=𝒯_{opt}(P,L)`$, and let $`W_R=weight(𝒯_{opt}^L,R)`$ be the weight of $`𝒯_{opt}^L`$ under the crossing metric of $`R`$. Clearly, $`E[W_R]=𝒲_{opt}(P,L)\frac{r}{n}`$. Thus, we know that with probability $`0.9`$ we have $`W_R10𝒲_{opt}(P,L)\frac{r}{n}`$ (by Markov inequality), and with probability $`0.9`$, we have that $`𝒲_{opt}^R=𝒲_{opt}(P,R)W_R10𝒲_{opt}(P,L){\displaystyle \frac{r}{n}}`$.

Let $`p,qP`$ be two points, and let $`X_{pq}`$ be the distance between $`p,q`$ in the arrangement $`𝒜(R)`$. If the distance between $`p,q`$ is large, that is $`U=𝒟_L(p,q)c_0(n/r)\mathrm{log}n`$ (where $`c_0`$ is a large enough constant), then one can show using Chernoff inequality, that with high probability, we have:

$$\frac{U}{2}X_{pq}\frac{n}{r}2U.$$

On the other hand, by the above argument, each edge $`e=pq`$ of $`𝒯_{opt}^R=𝒯_{opt}(P,R)`$ either intersects at most $`c_0(n/r)\mathrm{log}n`$ lines of $`L`$, or alternatively, the number of lines of $`L`$ intersected by $`e`$ is smaller than $`2(n/r)X_e`$, where $`X_e`$ is the number of lines of $`R`$ that $`e`$ intersects. Thus, with high probability, we have

$`𝒲_{opt}(P,L)`$ $``$ $`weight(𝒯_{opt}^R,L)={\displaystyle \underset{e=pq𝒯_{opt}^R}{}}𝒟_L(p,q){\displaystyle \underset{e𝒯_{opt}^R}{}}\left(c_0{\displaystyle \frac{n}{r}}\mathrm{log}n+2X_e{\displaystyle \frac{n}{r}}\right)`$
$`=`$ $`c_0{\displaystyle \frac{n^2\mathrm{log}n}{r}}+𝒲_{opt}(P,R){\displaystyle \frac{2n}{r}}.`$
###### Remark A.4

We can make both probabilities in Lemma A.3 large by repeating the experiment $`O(\mathrm{log}n)`$ times, and picking the smallest $`W(P,R)`$ computed. With high probability, we have

$$\frac{n}{r}\frac{𝒲_{opt}(P,R)}{10}𝒲_{opt}(P,L)\frac{n}{r}\left(c_0n\mathrm{log}n+2𝒲_{opt}(P,R)\right).$$

In particular, if $`𝒲_{opt}(P,R)>c_0n\mathrm{log}n`$, we get that $`3𝒲_{opt}(P,R)\frac{n}{r}`$ is a constant factor approximation to $`𝒲_{opt}(P,L)`$.

###### Lemma A.5

Let $`r`$ be a prescribed parameter, and $`𝒲_{opt}=𝒲_{opt}(P,L)`$. Then, an algorithm can decide whether

* $`𝒲_{opt}`$ is small - namely $`𝒲_{opt}\frac{10c_0n^2\mathrm{log}n}{r}`$.
* $`𝒲_{opt}`$ is large - $`𝒲_{opt}=\mathrm{\Omega }(\frac{n^2}{r}\alpha (n)\mathrm{log}^2n)`$.
* $`𝒲_{opt}`$ is in between. Any of the two above answers are valid.

The algorithm takes $`O(n\alpha (n)\mathrm{log}^4n)`$ time, and returns a correct result with high probability.

###### Proof.

We pick $`m=O(\mathrm{log}n)`$ samples $`R_1,\mathrm{},R_m`$ by picking each line with probability $`r/n`$ into the sample. For each sample, we check whether $`𝒲_{opt}(P,R_i)10c_0n\mathrm{log}n`$, using the algorithm of Lemma A.1. This will require $`O(n\alpha (n)\mathrm{log}^3(n))`$ time for each sample, and $`O(n\alpha (n)\mathrm{log}^4(n))`$ overall.

If the algorithm of Lemma A.1 returned not large for any sample $`R`$, we know that $`𝒲_{opt}(P,R)=O(n\alpha (n)\mathrm{log}^2n)`$. And by Lemma A.3, we know that $`𝒲_{opt}(P,L)=O\left(\frac{n^2\alpha (n)\mathrm{log}^2n}{r}\right)`$ with high probability. ∎

Now, we can perform a binary search to approximate the weight of $`𝒲_{opt}(P,L)`$.

###### Lemma A.6

One can compute in $`O(n\alpha (n)\mathrm{log}^5n)`$ time a value $`M`$, so that

$$𝒲_{opt}(P,L)M=O(n\alpha (n)\mathrm{log}^2n+𝒲_{opt}(P,L)\alpha (n)\mathrm{log}n).$$

###### Proof.

Use Lemma A.5, set $`r_0=n`$. In the $`i`$-th iteration check whether $`𝒲_{opt}=\mathrm{\Omega }\left(\frac{n^2}{r_i}\alpha (n)\mathrm{log}^2n\right)`$, by using the algorithm of Lemma A.5. If it is, we set $`r_{i+1}=r_i/2`$, and repeat the process. We stop as soon as this check fails. Then, we know that with high probability

$$\frac{10c_0n^2\mathrm{log}n}{r_{i1}}𝒲_{opt}(P,L)=O\left(\frac{n^2\alpha (n)\mathrm{log}^2n}{r_i}\right)=M,$$

implying that $`M`$ is the required approximation. ∎

###### Remark A.7

Note, that if algorithm of Lemma A.6 stops after the first iteration, then $`𝒲_{opt}=O(n\alpha (n)\mathrm{log}^2n)`$. In such a case the approximation we get is much worse then logarithmic. However, this is to some extent the easiest case: Without any sampling we get a spanning tree of near linear (or sub linear) weight.


text:
# Local induction approximation in the theory of superfluid turbulence. Numerical consideration

## 1 Introduction.

It is well known that chaotic vortex structures, appearing in volume of superfluid helium under the particular conditions influence on hydrodynamic and thermodynamic properties of HeII. At this time, any explicit macroscopic equations to consider this influence are absent. In hydrodinamic approaches usually the Vinen equation is used. The system of hydrodynamic equations for superfluid turbulence, where the Vinen equation is incorporated, was obtained by some authors. (The details can be found for example in .) In turn, the Vinen equation is a phenomenological relation and it’s validity is, in principle, unknown. From this point of view, for an adequate description of various phenomena in HeII, it is desirable to consider exact microscopic equations of vortex dynamics, based on the Biot-Savart law. But any analytical progress in this approach is rather difficult (practically zero at this time). Numerical investigations also have some problems. The Biot-Savart law implies the nonlocal interaction and therefore large computer resourses are demanded. Since the works  to the present moment, the local induction approximation (LIA) of the Biot-Savart law is often used. But the validity of the LIA to describe vortex dynamics is the object of discussions until now.

This report devoted to the numerical investigation of the LIA and present the continuation of the early analytical work . In  was shown that the LIA can not describe a flux of amplitudes of harmonics across the spectrum in the one-dimesional Fourier-representation. This result seems rather unusual, taking into account that the LIA is a nonlinear differential equation. The aim of the paper is to check this with help of numerical modeling.

## 2 The problem statement.

The dynamical equation of the evolution of a vortex filament has the following LIA-form:

$$\frac{d\stackrel{}{s}(\xi ,t)}{dt}=\beta \stackrel{}{s^{}}\times \stackrel{}{s^{\prime \prime }}+\nu \stackrel{}{s^{\prime \prime }}$$
(1)
Here $`\stackrel{}{s}(\xi ,t)`$ is a radius vector of a vortex line point labeled by the variable $`\xi `$, $`\stackrel{}{s^{}}`$ is the derivative on the parameter $`\xi `$; $`t`$ is time; the quantity $`\beta `$ is the coefficient of the nonlinearity $`\beta =\frac{\kappa }{4\pi }\mathrm{log}\frac{R}{r_0}`$, with the circulation $`\kappa `$ and the cutting parameters $`R`$ (the external size, i.e. the averaged radius of curvature) and $`r_0`$ (the vortex core size). The coefficient of dissipation $`\nu `$ appears in the eq.(1) when external counterflow is absent . Below one supposed that $`\xi `$ is the arclength.

Amplitudes of one-dimensional Fourier harmonics are defined as:

$$\stackrel{}{s}_k(t)𝑑\xi s(t,\xi )e^{ik\xi }$$
(2)
The numerical calculations was executed for two cases. Firstly, the spectral evolution of some initial configuration was investigated. At that it was supposed $`\nu =0`$ in the eq.(1) and thus only nonlinear term was presented in the dynamical equation. The evolution of harmonic amplitudes was defined by this term only. Secondly, a random force term of the Gaussian type was added to the right hand side of eq.(1). Thus eq.(1) was modified as:

$$\frac{d\stackrel{}{s}(\xi ,t)}{dt}=\beta \stackrel{}{s^{}}\times \stackrel{}{s^{\prime \prime }}+\nu \stackrel{}{s^{\prime \prime }}+\stackrel{}{f}(t,\xi )$$
(3)
Here $`\stackrel{}{f}(t,\xi )`$ is defined by the correlator which in the one-dimensional $`k`$-space looks like:

$$<\stackrel{}{f}_k(t_1)\stackrel{}{f}_k(t_2)>=\frac{D}{k^y}\delta (t_1t_2)$$
(4)
## 3 The details of the numerical procedure

The technical detail: the figures, mentioned below are available in the individual files. For example the file Fig.1.gif corresponds to Fig.1 etc.

In this chapter the full details of the numerical procedures are described. A reader, not interesting this theme, may omit this section.

A vortex line was presented as a number of points in the coordinate space $`(x,y,z)`$. For example, the ring with radius 1 contains 1200-5000 points, situated on the equal distance ($`\mathrm{\Delta }s`$) from each other. The fourth order Runge-Kutta method was used for modeling the temporal evolution. The time step was defined as $`\mathrm{\Delta }t=(\mathrm{\Delta }s)^2/C`$. Where $`\mathrm{\Delta }s`$ is the minimal for the entire configuration space step, $`C`$ is a constant. In the work  one maintain that this method is stable when $`C>\sqrt{2}`$. The experimental testing in our case gives the result $`C5001000`$ if the vortex line is strongly intricated. The criterion of the explicity was the conservation of the full length and the integral of the curvature square. (This values are some of the LIA invariants.) For example, this integral was equal $`38222.2167961644`$ and it was covserved with the explicity $`\pm 0.0000001`$ during the program ran, the total length was conserved down to 12 decimal signs. Besides to test the explicity, the iterative Crank-Nicolson type numerical scheme was developed. This scheme was similar to the first one, described in the work . But the distinction was the following: the equation for $`\stackrel{}{s}`$ was used, instead of the equation for a tangent vector $`\stackrel{}{l}`$ as in . (There method is really unusable for the required explicity.) The criterion of the iteration convergence was about $`10^{18}`$. Let’s note that the explicity of the floating number representation in a standard computer is about 19 characters (the format ”double”). For example, the same initial configurations (it is shown at the Fig.1 on the (x,y) plane) were calculated with these methods till $`192547`$ and $`1174140`$ time steps, but for the approximately same physical time. One can see from Fig.2 and Fig.3 (the configurations on the $`(x,y)`$-plane, corresponding to these two calculations), the results are absolutely identical. The spectral distributions are also identical.

To create a strongly intricated initial configuration, the following procedure was used. Usual ring, with radius 1 was created at first. After that, sinusoidal disturbances were added to this configuration. The disturbances were added step by step, and after each step the procedure of the rescaling of the length was applied, besides the coordinates of points were corrected to place them uniformly. In this way (step by step) some harmonics were excited.

When the randomized model was studied, the additional $`\stackrel{}{f}(t,\xi )`$ was incorporated into the numerical method. At that this value was modeled in the following way:

$$f^\alpha (t,\xi )=\underset{n=n_{low}}{\overset{n=n_{high}}{}}\frac{c_t^\alpha }{(\frac{2\pi n\xi }{L_t})^{\frac{y}{2}}}sin(\frac{2\pi n\xi }{L_t}+\varphi _t^\alpha )$$

Here $`(\frac{2\pi n}{L})`$ is a one-dimensional wave vector $`k_n`$ with number $`n`$, $`L_t`$ is a total length of a vortex line (recalculated at each time step), $`\xi `$ is a current distance along the line, to clarify $`y`$ see the rel.(4), $`\alpha `$ denotes a spatial component, $`c_t^\alpha `$ is a random amplitude, defined by a random-number generator, $`\varphi _t^\alpha `$ is a random phase. The last two values were generated at each time step.

Besides, a random disturbance also was used to create an initial configuration. At that a bunch of harmonics was excited, and then the random force was switched off.

During the was program running, one was possible to observe the changing of various parameters and the spectral picture on the monitor. To develop the program, Borland C++Builder 5 tool was used.

## 4 The results, discussion and conclusions

Thus, in the previous paper  was shown that the LIA don’t describe a flux of harmonic amplitudes across the spectrum. The aim of the numerical experiment was to confirm this result. The first part of the experiment looks as following. One was created an initial strongly intricated configuration with some spectral structure. After that, the vortex filament evolved accordingly to the eq.(1) (under the condition $`\nu =0`$). The reconecting processes was omitted from the consideration, i.e. parts of a vortex line freely passed across each other when they intersected. The changing of the spectral picture was studed. The initial spectral distribution ploted on the Fig.4 from $`n=1`$ till $`n=150`$ (there and below a red line correspond to an initial state, black line denotes a current distribution). A wave vector $`k_n=2\pi n/L_t`$, where $`L_t`$ is the total length of the filament at the moment $`t`$. In coordinate spase this configuration plotted on the Fig.1 ((x,y)-plane). The successive steps of the spectral evolution are shown at the Fig.5-Fig.8. Finally, the spectral distribution looks as at the Fig.9. At later time a stationary state was formed and any changes were only inside the area, marked by blue line. In the stationary state, the average size of the vortex tangle ceased to grow and only slight pulses occured. Inside of the marked area some peaks appeared, which dassapeared after a time. The general spectral distribution was not changed. The harmonic excitation in the low region (about the first harmonic) was absent just as the excitation in the high area. Thus, after the stationary state was formed, an spectrum changes was ceased. The initial spectral distribution was changed cardinally. One can seem that this point is not in agreement with the result of . But the procedure, used in  implicity implies just stationary state. As well known, a standard situation look as following. A wave packet spreads to infinity when a nonlinear, nondissipative equation is used.

The stability of the harmonic distribution, perhaps connected with the infinite number of the invariants, produced by the LIA. For example, total length $`L=const`$, $`𝑑\xi \stackrel{}{s^{\prime \prime }}^2=const`$, etc. High harmonics can not excitate because it means the growth of the curvature and, hence, of the curvature square, this is impossible. Increasing of low harmonics lead to the growth of a loop size, this is also impossible because $`L`$ is a constant and so the total curvature must decrease in this case. Thus some transfer across the spectrum is certanly absent if a vortex configuration achieved a stationary state.

In the second part of the numerical experiment the temporary evolution of a detached narrow harmonic bunch was studied. The harmonics in the environment of this bunch was equal to zero. The initial configuration presented a ring. Some scores of time steps this ring undergo the random disturbances of the kind: $`f^\alpha (t,\xi )=_{i=n_{low}}^{i=n_{high}}\frac{c_t^\alpha }{(\frac{2\pi n\xi }{L_t})^{\frac{y}{2}}}sin(\frac{2\pi n\xi }{L_t}+\varphi _t^\alpha )`$. The amplitudes $`c_t^\alpha `$ and the phases $`\varphi _t^\alpha `$ possessed random values at each time steps. (Their values was defined by the random number generator.) As the result the initial spectral configuration was obtained (Fig.10). (At the Fig.10 the stripe of the spectrum is plotted from harmonic number $`n=70`$ up to number $`n=120`$.) Later the external force was switched off and the evolution under the LIA nonlinearity only was investigated.

The result is following. There are no any additional harmonics exciting in the system, on condition that they equal to zero initially. This result was checked thoroughly. The boundaries of the excited area was thoroughly examined after about 800 thousands time steps. They were examined at the enlarge scale factor $`10^8`$ (in the vertical direction) relative to the amlitude of the stripe. As it shown at the Fig.11 where the boundaries of the excited harmonic bunch plotted (scale factor $`10^8`$), there are no any disturbances beyong the initial stripe. At the same time the spectral distribution inside of the perturbed area was varied. This change represented the periodical process (the period is about 20 thousands time steps). The steps of the changes is shown on the Fig.12-Fig.18. Obviously, the spectral distribution strictly return to it’s initial state. The period of the return depends on the width of an excited strip and also on current numbers of the harmonics. Note again, in a standard nonlinear task we have a wave packet is smearing and tending to it’s stationary spectral distribution. Note, however, after some millions time steps, distortion of the spectral stripe was observed. Apparently it is produced not strictly zero initial background. It was tested for the smooth ring. Obviously, in the later case only the first harmonic must be present. But other harmonics had amplitudes of order $`10^7`$ (relative to the amplitudes, plotted on the Fig.12 - Fig.18). Because, any nimerical method has inadvertent errors, backgroud amplitudes was conserved within about 1 percent. At the Fig.19 (scale factor $`510^9`$) shown the divergence from zero for the initial (red) and current (black) backgrounds. Besides, the graphic on the Fig.11 is wider the on the Fig.12 - Fig.18 because the greater magnification in the vertical direction allows to see small amplitudes which are lost details.

Also a case when the external force not switched off was considered. The force was during all time when the programm was running. At that, of course, $`\nu 0`$ in the eq.(3). The result is the same. After the quasi-equilibrium between pumping and dissipation was formed (i.e. when the tendency for changing the total length and the total curvature vanished), any new harmonics not appeared in the system.

Thus two conclusions for sure following for the LIA.

1. There are no transfer across the spectrum when a vortex system is in the stationary state.

2. There are only initially excited harmonics remain in a vortex system. New harmonics never excite. All redistributions of hamronic amplitudes take place between excited harmonics.

And the last remark. It can be suppose that the question: ”whether can LIA approximately describe the superflud turbulence” is rather make no sense until we don’t know exactly the role of the processes of reconnections in dynamics of vortex loops. But this problem is too much difficult.

The work was carried out under support of INTAS (grant N 2001-0618) and RFBR (grant N 03-02-16179).