text:
# The Flattened, Rotating Molecular Gas Core of Protostellar Jet HH 212

## 1 Introduction

The youngest protostars are often so deeply embedded in their parental cloud cores that shocks driven by their jets cannot be easily traced in visible shock line emission but are bright in near-infrared molecular hydrogen shock lines. Such is the case for the recently discovered highly collimated jet HH 212, near the Horsehead Nebula region of Orion (Zinnecker, McCaughrean, & Rayner 1998).

HH 212 is remarkable because of its high collimation and symmetric pairs of bow shocks and shock knots on either side of the driving source (Figure 1). Theories of protostellar jet structure often involve a heavy influence of the surrounding medium, with impacts creating shocked globules along the jet, sometimes even redirecting it. But the symmetry of HH 212 suggests that some pulsing mechanism from the source itself is responsible for the nearly periodic bow shock pairs. Somehow the symmetry is maintained even though the jet is moving through a dense gas core. In the innermost region, water masers have been found in motion along the outflow (Claussen et al. 1998).

The inferred position of the jet exciting source coincides with IRAS point source 05413-0104, which was detected at 25, 60, and 100 microns; it is still too cold and embedded to be detected strongly at 12 microns (see Beichman et al. 1986). It is surrounded by a cold (11.5 K), compact ammonia core out of which it presumably formed (Wouterloot, Walmsley, & Henkel 1988). The spectral energy distribution climbs to longer wavelengths. It is detected also at 1.3 mm, with a submillimeter to bolometric luminosity of about 2% (Chini et al. 1997). It has narrow NH3 linewidths compared to emission found at the position of other maser sources. These characteristics suggest that the source is very young.

HH212 is classified as a Class 0 system, the youngest class of protostars with much material to be accreted still in the surrounding molecular gas envelope (André, Ward-Thompson, & Barsony 2000). Class 0 sources have a high rate of submillimeter to bolometric luminosity (over 0.5%) and high rates of infall and outflow (André, Ward-Thompson, & Barsony 1993, Adams, Lada, & Shu 1987). The morphology and kinematics of the surrounding gas at this early stage of the formation of a star are crucial to understand as an important link between prestellar gas cores and the resulting young T Tauri star/disk systems. Our VLA ammonia observations were conducted as a study of the gas density, temperature, and velocity distributions in this pristine, highly embedded young source.

## 2 Observations

Observations were carried out with the Very Large Array (VLA) of the National Radio Astronomy Observatory111The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.. The ammonia (J,K) = (1,1) and (2,2) inversion transitions were observed simultaneously, at rest frequencies of 23.694495 GHz and 23.722733 GHz, respectively. The width of the antenna primary beam response was approximately 2’ FWHP. The array was in its compact, “D” configuration. The synthesized beamsize was $`8.4^{\prime \prime }\times 8.8^{\prime \prime }`$ FWHM, with natural weighting and a 20 k$`\lambda `$ taper applied to the UV data. The bandwidth was 3.125 MHz and the channel separation was 24.4 kHz, corresponding to 0.31 km sec-1.

The absolute flux calibrator was 3C286, with a calculated flux density of 2.41 Jy at 1.3 cm. The bandpass calibrator was 3C84, and a phase closure solution was applied based on interleaved observations of the phase calibrator B0605-085 (1950). The “dirty” beam was deconvolved from the maps using the “Imager” algorithm in the NRAO AIPS image processing software.

## 3 A Flattened, Rotating, and Heated Gas Envelope

The integrated intensity of the NH3 (1,1) line is displayed in the contours of Figure 1, overlaid on the 2.12 $`\mathrm{\mu m}`$ image of HH 212 from Zinnecker, McCaughrean, & Rayner (1998). The gas core is flattened with an axis ratio of about 2:1. The major axis has a FWHM extent of 29”, corresponding to 12000 AU at an assumed distance of 400 pc, with full extent of over 15000 AU; it is centered on the jet source, and it lies perpendicular to the jet which has a position angle of about 25 degrees east of north and only 4$`\mathrm{deg}`$ from the plane of the sky (Claussen et al. 1998).

A velocity gradient is apparent across the flattened gas core. In Figure 2a the integrated intensity (moment 0) is again plotted in contours, while the underlying shades represent the (moment 1) velocity field. Radial velocities change smoothly from blue-shifted gas in the northwest side to redshifted gas in the southeast; the large velocity gradient of about 4-5 km sec-1pc-1is virtually constant. The direction of the gradient within the central FWHM intensity region of the core is perpendicular to the jet axis, with the nodes of constant velocity roughly parallel to the jet.

Along the eastern side of the core, a small amount of blueshifted material is seen to extend to the northeast alongside the northern, blueshifted side of the jet. To a lesser extent, the same effect is also seen to the south with a small protrusion of redshifted gas extending in the same direction as the redshifted side of the jet. Molecular gas is possibly being swept along by the outflow. There is even a hint, at least in the north, that molecular gas is found in a “bowl” surrounding the jet axis, possibly outlining a cavity of a wider-angle flow surrounding the collimated shock emission of the jet (see, e.g., Shu et al. 1995, Shang, Shu,& Glassgold 1998).

Ammonia inversion lines are a convenient probe of both optical depth and temperature. Most ammonia inversion transitions are composed of multiple hyperfine components. The relative intensities of central and satellite component blends of the (1,1) transition can be used to calculate the optical depth. For the main component and total transition of the (1,1) line, we find $`\tau `$ values of 0.8 and 1.6, respectively, at the core center. Using the procedures and assumptions described in Ho & Townes 1983 (see also Wiseman & Ho 1998), we also estimate an NH3 column density of $`1.6\times 10^{14}`$ cm-2 and a total volume density $`n`$ of $`1\times 10^5`$ cm-3. The mass of the core gas detected with the interferometer is approximately 0.2 M, assuming an NH3/H2 abundance ratio of $`2\times 10^8`$ (Harju, Walmsley, & Wouterloot 1993). Within the large uncertainties of this abundance ratio, the mass is the same as that needed to gravitationally bind the rotating system.

The ratio of intensity of the NH3 (2,2) and (1,1) transitions reflects level populations and thus serves as a thermometer for rotational temperature within cold molecular gas cores. In Figure 2b, the contours again show the integrated intensity of (1,1) emission. The shades show the (2,2)/(1,1) ratio, indicating the relatively warmer regions. It is apparent that this is indeed a cold core of gas, with the (2,2) emission showing up only in the central region of the core. The (2,2)/(1,1) ratio peaks near the position of the embedded protostellar source, and extends slightly to the northwest. The peak rotational temperature here is 14 K, as calculated from the line ratio (Ho & Townes 1983). In colder regimes such as this, the kinetic temperature is approximately equal to the NH3 rotational temperature (Danby et al. 1988, Kuiper 1994). The young protostar is evidently still in the center of the gas core, and we may be seeing here the warm gas heated by the protostar and the jet as they interact with the environment.

## 4 Discussion

### 4.1 Rotation in Protostellar Systems

Envelopes of Class 0 protostars carry much (often most) of the mass of the system and can extend to diameters of several thousands of AU. Interferometric maps of such systems tend to show flattened disks perpendicular to the axis of the associated jet (Myers, Evans, & Ohashi 2000). The flattened molecular gas envelope of the HH 212 jet source has a smooth and roughly constant velocity gradient perpendicular to the jet axis (Figure 2a, 2c).

Turbulence in cloud cores can produce velocity gradients that mimic rotation (Burkert & Bodenheimer 2000 (BB)). However, the gradient we find in HH 212 is faster than the most likely effect produced in the BB turbulence models. It should be noted that the BB models were based on more massive cores, so proper scaling of the models would be required to better assess the possibility of a turbulent component to the observed velocity gradient in HH 212. But it is also true that the HH 212 gradient extends systematically across the long axis of the core, and perpendicular to the jet, for over 5 synthesized beamwidths. Thus the combined evidence suggests that the gradient is more likely dominated by rotation.

The more evolved T Tauri systems often show Keplerian rotation signatures in the less plentiful envelope gas, with the radial velocities increasing and diverging toward the central source ($`VR^{1/2}`$). An early report of Keplerian rotation in a young stellar disk was made for 13CO observations of HL Tau (Sargent & Beckwith 1987, 1991), though in this source the kinematics proved to be much more complicated, with evidence for infall (Hayashi et al. 1993) and outflow (Cabrit et al. 1996). The strongest Keplerian velocity curves are seen in circumstellar disks of a few hundred AU in Class II T Tauri systems without much envelope material left, thus having little confusion from large scale infall or outflow motions (Mundy, Looney, & Welch 2000, Dutrey et al. 1998). Keplerian rotation has been detected in the gas disks of GM Aur (Dutrey et al. 1998) and DM Tau (Guilloteau & Dutrey 1998), for example, as well as in multiple systems such as UY Aur and the beautiful ring of GG Tau (Duvert et al. 1998, Guilloteau, Dutrey, & Simon 1999).

Several Class 0 systems such as HH 212 are in dense gas cores that appear to have a more rigid rotation signature (Myers, Evans, & Ohashi 2000); this is consistent with their character of having much of the system mass still distributed in the gas envelope around the embedded accreting protostar(s). Indeed, the 0.2 M we detect with the VLA in the HH 212 core is 10 times the circumstellar gas mass found in some T Tauri systems (e.g. $`0.01`$M for UY Aur (Duvert et al. 1998)). In such early systems, a gradient in the intensity-weighted mean velocity that is constant rather than diverging toward the center can indicate that either a) the non-Keplerian rotation curve is reflecting a mass distribution that is not yet dominated by a central condensation, or b) the system is in Keplerian rotation, but high optical depths result in velocity maps that are dominated by motion in the outer layers along the line of sight toward the center of the core. Models of expected channel map emission can be used to test for Keplerian rotation (e.g. Koerner, Sargent, & Beckwith 1993, Dutrey, Guilloteau, & Simon 1994). Higher spatial and spectral resolution maps would be needed to fully model the HH 212 envelope, although the slightly broader line widths we observe toward the center of the core are possible evidence for some Keplerian component to the motion. In any case, it now appears that the observed kinematics of the protostellar envelope can signify the developmental stage of the embedded YSO. More high resolution kinematical studies are needed to confirm this pattern. The large envelopes such as we observe in HH 212 may feed mass directly to the embedded protostar (if at the earliest stage) or to the inner circumstellar accretion disk just beginning to form. It is important to note that this inner disk ($`<100`$ AU and thus unobservable in our NH3 maps) probably rotates with near-Keplerian velocities (Stahler 2000).

Many rotating cores of these large sizes ($``$12000 AU) do not appear to conserve specific angular momentum, defined as Vrot $`\times `$ Rrot, where Vrot is the velocity of rotation at radius Rrot (Myers, Evans, & Ohashi 2000). Our results show that, with a specific angular momentum of 0.012 km s-1 pc at a radius of 6000 AU, the HH 212 gas envelope would fit well with other cores of similar size in a plot of Ohashi et al. (1997) showing specific angular momentum falling with decreasing source radius. However, Ohashi et al. showed that there is a transition at about this same size scale; smaller dense core regions do show conserved specific angular momentum in regions of both massive infall and rotationally supported disks. Thus this would serve as the size scale for dynamical collapse, where angular momentum is conserved. The HH 212 core seems to be in fast, bound rotation at the size scale of this transition.

Since the velocity gradient and the long axis of the envelope are both perpendicular to the jet axis, it is suggestive that rotation may contribute to the observed flattening. We calculate rotational energy of the system to reach up to about 3-10% of each of the thermal, turbulent, and gravitational energy (Harju et al. 1993, Myers 1983). Preferential contraction along an axis perpendicular to the disk plane, along possible magnetic field lines, could also produce flattened morphology. Another intriguing possibility is that the outflow itself contributes to the morphology by clearing out dense gas along the flow. Indeed, the blue-shifted wisp of gas extending to the northeast alongside the (northward) blue-shifted side of the jet (Figure 2a) seems to form part of a “bowl” around the jet and the associated CO outflow recently studied by Lee et al. (2000); part of the blue lobe of the CO flow shifts to the east corresponding to the NH3 wisp. It has been suggested that molecular outflows associated with jets may broaden over time, as has been seen in B5-IRS1; the broadening flow could contribute to the eventual decline of accretion (Velusamy & Langer 1998). Thus we may be seeing the HH 212 flow in an early stage of broadening and sweeping out envelope material.

### 4.2 Angular Momentum Release Through the Jet

Bipolar outflows appear to be a ubiquitous product of the formation of stars. Strong winds from young stellar sources are collimated into jets which are observed in visible and infrared shock emission. Magnetic fields may be responsible for channeling and collimating the jets (see, e.g., Shu et al. 1995).

Outflows are thought to be a ’release valve’ for the angular momentum carried into the protostellar accretion region by infalling material. Otherwise, matter infalling from even a very slowly rotating circumstellar envelope will have too much conserved angular momentum to fall into the inner accretion region, unless sufficient angular momentum is lost through disk viscosity.

Observational evidence for jets as significant angular momentum carriers has been virtually nonexistent. Recently, however, results of slit spectroscopy on the shock knots in the HH 212 jet have given some evidence of jet rotation (Davis et al. 2000). Five of six knots studied show similar evidence of rotation; in the inner part of the southern flow, for example, a clear velocity gradient is seen. If this gradient represents rotation, the jet rotation is in the same directional sense as that of the rotating molecular gas envelope we observe in ammonia. A rough extrapolation of the large scale angular momentum in the envelope gas fits the observed gradient in the jet, assuming a fraction of infalling mass is channeled into the flow. Thus this may be the first observational evidence for a rotating jet carrying away the angular momentum of material falling in from a rotating envelope. Clearly more detailed observations and evidence of jet rotation and inner disk kinematics are needed.

We thank our referee and also P. T. P. Ho for insightful guidance. Support for J.W. was provided by NASA through Hubble Fellowship grant #HF-01115.01-98A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA under contract NAS 5-26555.


format: text

text:
# Using electronic structure changes to map the 𝐻-𝑇 phase diagram of 𝛼′-NaV2O5

## Abstract

We report polarized optical reflectance studies of $`\alpha `$-NaV2O5 as a function of temperature (4-45 K) and magnetic field (0-60 T). Rung directed electronic structure changes, as measured by near-infrared reflectance ratios $`\mathrm{\Delta }`$R($`H`$)=R($`H`$)/R($`H`$=0 T), are especially sensitive to the phase boundaries. We employ these changes to map out an $`HT`$ phase diagram. Topological highlights include the observation of two phase boundaries slightly below $`T_{SG}`$, enhanced curvature of the 34 K phase boundary above 35 T, and, surprisingly, strong hysteresis effects of both transitions with applied field.

After the discovery of a low-temperature spin-gapped phase in $`\alpha `$-NaV2O5, there was intense interest in this compound as a possible spin Peierls (SP) material . Reexamination of the 300 K structure showed, however, that the situation is more complex . At room temperature, all vanadium ions have a charge of +4.5 and form chains along the $`b`$ axis. Chains are connected by V-O-V rungs. The result is a 1/4-filled 2-leg ladder material with low-dimensional chain character and 1 spin per rung . The spin gap (SG) opens at $`T_{SG}`$ = 34 K, as evidenced by an isotropic drop in the susceptibility . The opening of the spin gap is accompanied by unit cell doubling in the $`a`$ and $`b`$ directions, and quadrupling in the $`c`$ direction . According to Lüdecke et al., the low temperature phase of $`\alpha `$-NaV2O5 also displays alternating modulated and unmodulated ladder chains , although Nakao et. al. report a fully zigzig ordered structure, with variable temperature charge disproportionation . Early NMR experiments pointed out the importance of charge ordering (CO) in this material, as evidenced by both V4+ and V5+ resonances . Zigzag CO is the leading candidate for the low-temperature charge arrangement . More recent NMR experiments have discovered the advent of the charge-ordered phase in $`\alpha `$-NaV2O5 at $`T_{CO}`$ = 37 K, slightly above the 34.7 K lattice distortion and spin gap formation . Early thermal expansion measurements also discerned two transitions . The detailed interaction between charge and magnetic ordering in $`\alpha `$-NaV2O5 is not understood at this time, although, because of the larger energy scale, CO is expected to dominate. Further, the field dependence of $`T_{SG}`$ (on a field vs. temperature ($`HT`$) phase diagram) is anomalous  and does not follow the expected $`H^2`$ behavior for a SP transition, at least at low applied fields . The nearly vertical field dependence of $`T_{SG}`$ and the interplay of CO and magnetic effects have been cited as reasons why $`\alpha `$-NaV2O5 should not be considered a true SP material.

The electronic structure of $`\alpha `$-NaV2O5 has attracted a great deal of attention , as have other physical properties . The optical conductivity is strongly polarized, as expected for such an anisotropic material; excitations in the rung direction ($`a`$) are much stronger than those in the chain direction ($`b`$). Electronic features are observed near 1 eV, between 3 and 4 eV, and above 4 eV in both polarizations. The 1 eV structure in the $`a`$ direction has received the most attention and is currently assigned as charge transfer excitation between V4+ and V5+ on opposite sides of a rung . That changes in electronic structure are sensitive probes of phase boundaries was established by previous work on CuGeO3 and $`\alpha `$-NaV2O5, where integrated intensities of reflectance and transmittance ratio features, as a function of temperature or magnetic field, display “breakpoints” at a phase boundary . These inflection points provide a useful tool to locate a transition in $`HT`$ phase space. In CuGeO3, the phase diagram obtained in this manner is consistent with that obtained with other techniques . At the same time, the detailed spectral features provide microscopic information on the changes which occur at a phase boundary. In the past, the challenge of accessing the high field phase has prevented general extension of this technique to $`\alpha `$-NaV2O5.

In order to provide further information on the electrodynamics of quarter-filled magnetic oxides, we have carried out polarized optical reflectance measurements to probe the high field response of $`\alpha `$-NaV2O5. We use these data to generate an $`HT`$ phase diagram which shows unusual topology and a number of important differences from traditional SP behavior.

High quality single crystals of $`\alpha `$-NaV2O5 were grown in Orsay by the flux method . Typical crystal dimensions were 1.65 $`\times `$ 2.8 $`\times `$ 0.25 mm. The high-field optical reflectance measurements were performed in the 60 T Long-Pulse magnet at the National High Magnetic Field Lab in Los Alamos, NM. This magnet is powered by a 1.4 GVA motor-generator which provides a very long (2 second), user-definable magnet pulse. We elected to use the flat-top field profile, in which the magnet is held at maximum field (60 T) for a full 100 ms. The sample resides in a vacuum jacket in the bore of the magnet, and data were taken in the range 4 to 45 K, concentrating in the regime around the 34 K phase transition. The temperature was held constant for each shot; the sample was thermally cycled to above 40 K after each magnet pulse to avoid possible hysteresis/fatiguing effects. White light was coupled to the sample via a single 600 $`\mu `$m optical fiber, and thin-film, in-situ linear polarizers selected the desired polarization along the $`a`$\- (rung) or $`b`$\- (chain) direction. The reflected light was collected by another optical fiber, dispersed in an optically fast 0.3 m spectrometer, and collected by a high-speed CCD camera. Because the optical spectra of $`\alpha `$-NaV2O5 is modified in broad bands, a 150 line/mm grating was employed. We covered the wavelength range from 350 to 1200 nm. Full spectra were collected every 2.1 ms throughout the entire 2 second magnet pulse, providing for complete field-dependence in a single shot at fixed field . A calibrated coil, mounted on the probe, measured the applied field, which was applied perpendicular to the $`ab`$ plane.

Small field-induced changes in the polarized optical reflectance were studied by calculating reflectance ratios, $`\mathrm{\Delta }`$R($`H`$) = R($`H`$)/R($`H`$=0 T), at a number of fixed temperatures. Reflectance ratios were measured in 1 or 1.5 K steps around $`T_{SG}`$ (34 K) and at larger temperature intervals further away from the transition. We quantified these reflectance ratio changes with applied field by integrating the spectral area. Inflection points in the integrated area were used to identify phase boundaries in $`HT`$ space. Error bars on the location of the high-field phase boundaries vary depending on the sweep; we estimate the largest to be on the order of $`\pm `$ 3 T on the 30 K point near 60 T; other error bars are $`\pm `$ 2 T.

Figure 1 displays the rung-polarized ($`a`$ axis) near-infrared reflectance ratios of $`\alpha `$-NaV2O5 with applied field. As in previous temperature sweep data, the electronic structure is modified in wide bands, with the largest changes occurring in the vicinity of the 1 eV excitation . That the 1 eV excitation is sensitive to an applied magnetic field is related to the charge transfer from vanadium centers across the V-O-V rung, which also involves a shift of the spin to a new site. Chain-polarized and higher energy rung-polarized spectral signatures were also observed at the phase boundaries; we concentrate our analysis and discussion on the rung-directed near-infrared data because the effects are more pronounced.

In order to quantify the spectral changes in $`\alpha `$-NaV2O5 with magnetic field, we calculated the spectral area in the wavelength range of interest and plotted it vs. applied field. Representative data are shown in Fig. 2. The most interesting spectral changes are in the vicinity of the 34 K transition, where there are large hysteresis loops in the integrated area vs. field data below $`T_{SG}`$. That inflection points in the optical response are different on the upsweep compared to the downsweep indicates first-order behavior. Data collected above the 34 K transition temperature does not display hysteresis (Fig. 2). In contrast, CuGeO3 shows second-order character over the majority of $`HT`$ phase space with first order behavior only around the critical field ($`H_c`$=12.5 T) at low temperature, where the spin gap is fully open and the lattice distortion is most pronounced . Previous investigations of $`\alpha `$-NaV2O5 generally report second-order behavior , although these efforts have not explored the $`HT`$ diagram in the field regime reported here. The thermal expansion data are an exception in that two transitions are observed .

Figure 3 displays the $`HT`$ phase diagram for $`\alpha `$-NaV2O5 derived from an analysis of the field dependent optical properties. The field dependence of $`T_{SG}`$ is nearly vertical close to 34 K, softening only above 35 T. Note that previous 28 and 30 T determinations of the transition temperature are in general agreement with this diagram . Two phase boundaries are clearly observed below $`T_{SG}`$ in $`\alpha `$-NaV2O5. That there are indeed two phase boundaries is supported by the spectral area vs. field data in Fig. 2. Here, two distinct inflection points are seen in both up- and down-field sweeps; the position of these breakpoints changes with temperature, indicating that they are characteristic of $`\alpha `$-NaV2O5 and not an experimental artifact . The low-field boundary in Fig. 3 separates two phases directly below $`T_{SG}`$; it saturates near 33 T. This structure is reminiscent of, but different than, the critical field ($`H_c`$) boundary that separates dimerized and incommensurate phases in traditional SP materials . The optical signature of the 33 T boundary weakens rapidly away from $`T_{SG}`$, as indicated on the phase diagram, making the distinction between the two phases short-lived and limiting the range of phase space available for study . In contrast, the optical signature of the high field phase boundary remains strong within the range of our investigation (up to 60 T). The curvature is weaker than that in CuGeO3 , but stronger than that previously reported for $`\alpha `$-NaV2O5 . Based upon the three highest field data points and the work of Bulaevski and Cross  which predicts $`H^2`$ dependence to the phase boundary, we extract $`\alpha `$$``$0.1. The first-order behavior of this (and the 33 T) phase boundary is unexpected and different from CuGeO3 , although recent theoretical work proposes a framework under which hysteresis might be understood .

In summary, we report polarized optical reflectance studies of $`\alpha `$-NaV2O5 between 350 and 1200 nm as a function of temperature (4-45 K) and applied magnetic field (0-60 T). Using inflection points in the $`a`$-axis near-infrared reflectance ratio data, we have identified high-field phase boundaries at a number of different temperatures and mapped out an $`HT`$ phase diagram. Highlights include the observation of two phase boundaries slightly below $`T_{SG}`$, enhanced curvature of the 34 K phase boundary above 35 T, and, surprisingly, strong hysteresis effects of both transitions with applied field. We hope these experiments will spur complementary magnetization and NMR studies of $`\alpha `$-NaV2O5 in the future.

Acknowledgements

We are grateful for financial support from the Materials Science Division, Basic Energy Sciences at the U.S. Department of Energy under Grant # DE-FG02-99ER45741. The National High Magnetic Field Laboratory in Los Alamos, NM, is supported by the National Science Foundation, the State of Florida, and the Department of Energy. We thank Bill Moulton for useful discussions.

FIGURES

Fig. 1. Near-infrared reflectance ratio ($`\mathrm{\Delta }`$R) spectra of $`\alpha `$-NaV2O5 at 31 K polarized in the rung ($`a`$) direction. The $`\mathrm{\Delta }`$R spectra shown here are presented at three representative fields: 60 T/0 T, 45 T/0 T, and 0 T/0 T, illustrating the different optical response in the three regimes of the $`HT`$ phase diagram at this temperature.

Fig. 2. Integrated spectral area of $`\alpha `$-NaV2O5 vs. applied magnetic field at four different temperatures below $`T_{SG}`$ and one temperature above $`T_{SG}`$. Solid arrows indicate breakpoints; dashed arrows denote the field sweep direction.

Fig. 3. Solid symbols: $`HT`$ phase diagram of $`\alpha ^{}`$-NaV2O5, as obtained from an analysis of inflection points in the reflectance ratio data on the upsweep of the magnet. Symbol size is intended to allow visualization of transition strength. Open star: 28 T data as determined by Long et. al. in Ref. ; open squares: data as determined by Bompadre et. al. in Ref. ; open circles: data as determined by Schnelle et. al. in Ref. . Dashed lines guide the eye.


text:
# Notes on computing peaks in k-levels and parametric spanning trees

## 1 Introduction

The $`k`$-level of an arrangement of lines is one of popular geometric objects in computational geometry . The $`k`$-level is the union of $`k`$-th lowest (closed) line-segments of the arrangement, and it can be considered as the trajectory of the $`k`$-th smallest element in a set of $`n`$ data each of which depends on a parameter $`x`$ linearly. Thus, the $`k`$-level is a special case of the locus of the largest element of the minimum base of a parametric matroid  with one parameter; in precise, the $`k`$-level is the locus of the maximum element in the minimum base of a parametric uniform matroid of rank $`k`$. From a different aspect, the $`k`$-level is a dual concept of the $`k`$-set , and the complexity $`g_k(n)`$ of the $`k`$-level of an arrangement of $`n`$ lines is asymptotically same as the number of different $`k`$-sets in a set of $`n`$ points in a plane.

Lovász  first gave a nontrivial $`O(n^{3/2})`$ upper bound for $`g_{n/2}(n)`$, and also introduced Straus’s $`\mathrm{\Omega }(n\mathrm{log}n)`$ lower bound. The current best upper and lower bounds for $`g_k(n)`$ are $`O(k^{1/3}n)`$  and $`2^{\mathrm{\Omega }(\sqrt{\mathrm{log}k})}n`$ , respectively. The upper bound holds for any parametric matroid (with a linear parameter) of rank $`k`$ in $`n`$ elements . Moreover, Cole et al. gave an $`O(g_k(n)\mathrm{log}^2n)`$  time algorithm to compute the $`k`$-level; this time complexity has been improved due to improvement of algoirthms for the dynamic convex hull computation employed as a key subroutine . The current best time complexity (randomized expected time) is $`O(g_k(n)\alpha (g_k(n))\mathrm{log}n)`$ time given by Har-Peled , where $`\alpha ()`$ is the inverse Ackermann function. Thus, the time complexity is very close to $`O(g_k(n)\mathrm{log}n)`$.

However, we often need a compact “outline” of a trajectory of a parametric problem by using a small number of characteristic points on it. Such an outline, generally speaking, will be useful as a compact data to control parametric problems, and possibly utilized in designing kinetic data structures . Local peaks in the trajectory are considered to be natural characteristic points. A key observation to investigate the $`k`$-level is that it has at most $`2k1`$ local peaks (at most $`k`$ maximal peaks and at most $`k1`$ minimal peaks) . One interesting problem is to compute all the local peaks efficiently. This enables us to give a decomposition of the $`k`$-level into monotone chains, and hence create an outline with a size $`O(k)`$ of the $`k`$-level.

We give an algorithm to compute all the local peaks in the $`k`$-level of an arrangement of $`n`$ lines in $`O(n\mathrm{log}n)+\stackrel{~}{O}((kn)^{2/3})`$ time, where $`\stackrel{~}{O}`$ is the big-O notation ignoring polylogarithmic factors. The current estimate for the polylogarithmic factor of the second term is less than $`\mathrm{log}^5n`$; however we do not give it explicitly in this paper, since it is probably loose and will confuse readers. The time complexity is better than $`O(g_k(n)\mathrm{log}n)`$ for some restricted range of $`k`$ even if the current lower bound of $`g_k(n)`$ by Tóth  is tight. If we substitute the current $`O(k^{1/3}n)`$ upper bound to $`g_k(n)`$, the time complexity is better than $`O(g_k(n)\mathrm{log}n)`$ if $`k=O(n/\mathrm{log}^cn)`$, where $`c`$ is a suitable constant.

Another interesting question is how fast we can compute $`\tau `$ largest maximal peaks for $`\tau k`$. If $`\tau =1`$, Roos and Widmayer  gave a neat method to compute the maximum point in the $`k`$-level in $`O(n\mathrm{log}n+(nk)\mathrm{log}^2(nk))`$ time by using an efficient slope selection algorithm. We can compute $`\tau `$ largest peaks in $`O(n\mathrm{log}^2n)+\stackrel{~}{O}((\tau n)^{2/3})`$ time by combining Roos and Widmayer’s technique and the above mentioned method for computing all the peaks.

Finally, we investigate whether we can analogously treat some parametric matroids: Compute peaks in the trajectory of the largest element in the minimum weight base of a parametric matroid. In particular, the graphic matroid is of wide interest : Consider a weighted undirected connected graph $`G(x)`$ with $`k`$ nodes and $`n`$ edges, such that each edge has a parametric weight that is linear in a parameter $`x`$. Here, $`k`$ and $`n`$ become the rank and size of the graphic matroid, respectively. Let $`T(x)`$ be the minimum weight spanning tree of $`G(x)`$ and consider the longest edge $`e(x)`$ in $`T(x)`$. Note that the minimum weight spanning tree becomes a spanning tree that minimizes the length of the longest edge. We call the edge $`e(x)`$ the spanning bottleneck edge(SBE), and write $`SBE(x)`$ and $`w_{SBE}(x)`$ for $`e(x)`$ and its weight, respectively. The naming comes from the fact that $`w_{SBE}(x)`$ is the minimum value of $`w`$ such that the subgraph of $`G(x)`$ consisting of edges with weights less than or equal to $`w`$ is connected.

The following problems are important in sensitivity analysis: (1). Compute the maximum value and the minimum value of $`w_{SBE}(x)`$ for $`xI`$, where $`I`$ is a given interval. (2). Compute all peaks of the trajectory $`y=w_{SBE}(x)`$.

For example, imagine a system represented by the graph $`G`$ where a link represented by an edge with a weight larger than a (controllable) threshold value becomes unreliable, and the edge weight depends on a parameter $`x`$ linearly within an interval $`I`$. For a given subinterval $`JI`$, we want to know the threshold value of the edge weights so that the graph remains connected for every $`xJ`$. This can be reduced to the problem (1). Moreover, if we have computed all peaks in $`I`$ as a preprocessing (problem (2)), we can efficiently query for the threshold value, provided that we have an efficient method (shown in Section 4) to compute $`w_{SBE}(x)`$ at endpoints of $`J`$.

Both of problems (1) and (2) can be solved by computing the whole picture of the transitions of the minimum spanning trees, and the time complexity of the computation is $`O(kn\mathrm{log}n)`$ by using an algoirthm given by Fernandez-Baca et al. . Note that $`O(h_k(n)n^{2/3}\mathrm{log}^{3/4})`$ time output-sensitive algorithm of Agarwal et al. is better for some range of $`k`$, where $`h_k(n)`$ is the number of of transitions of the minimum spanning tree. An $`O(k^{1/3}n)`$ upper bound  and an $`\mathrm{\Omega }(n\alpha (k))`$ lower bound  are known for $`h_k(n).`$

Roos and Widmayer’s method can be directly applied to the first problem. By using dynamic maintenance algorithms  of a minimum spanning tree, the time complexity becomes $`O(\sqrt{k}n\mathrm{log}n)`$. Combined with range searching techniques, we improve the time complexity to $`\stackrel{~}{O}(n^{8/7}k^{1/7}+nk^{1/3})`$. We give some discussion on the second problem, although theoretical improvement on the $`O(kn\mathrm{log}n)`$ time method remains open.

## 2 Preliminaries

### 2.1 Roos and Widmayer’s algorithm

Given a set $``$ of $`n`$ lines in the $`x`$-$`y`$ plane, let $`_k`$ be the $`k`$-level of the arrangement of $``$. Let $`𝒑`$ be a point on $`_k`$ that has the maximum $`y`$-value $`y_{max}`$. Without loss of generality, we assume that such a point is unique. For any given value $`\alpha `$, one can decide whether $`y_{max}\alpha `$ or not in $`O(n\mathrm{log}n)`$ time: We sweep on the line $`h:y=\alpha `$ from the leftmost intersection point to the right to compute the levels of all intersection points on $`h`$ with lines in the arrangement. If all intersection points are above the $`k`$-level, $`\alpha >y_{max}`$; otherwise, $`\alpha y_{max}`$. By using this decision method, a binary search algorithm works to compute $`𝒑`$, and a weakly polynomial time algorithm with a time complexity $`O(n\mathrm{log}n\mathrm{log}\mathrm{\Gamma })`$ can be obtained, if each coefficient of the lines is a quotient number of integers with $`\mathrm{log}\mathrm{\Gamma }`$ bits. Roos and Widmayer applied an efficient slope selection method  to transform the binary search algorithm into strongly polynomial, and gave an $`O(n\mathrm{log}^2n)`$ time algorithm. They further improved the time complexity to $`O(n\mathrm{log}n+k\mathrm{log}^2k)`$ for computing the minimum and $`O(n\mathrm{log}n+(nk)\mathrm{log}^2(nk))`$ for computing the maximum.

### 2.2 Range query and Matoušek’s point set decomposition

We use well-known (although sophisticated) simplex range query data structures : We construct a data structure for a set $`S`$ of $`n`$ points in a plane such that given a query halfplane $`H`$, we can compute the number of points in $`S`$ located in $`H`$ efficiently. If we spend $`O(m)`$ time for constructing the data space for $`n\mathrm{log}nmn^2`$, the query time is $`\stackrel{~}{O}(n/m^{1/2})`$. The query can be done in polylogarithmic time by using $`O(n/m^{1/2})`$ processors. The data structure uses $`\stackrel{~}{O}(m)`$ space, although we do not discuss space complexity in this paper. Moreover, we can query the number of points in the intersection of two (or three) halfplanes in the same query time if we ignore a polylogarithmic factor. We can also do reporting query by spending additional $`O(N)`$ time if the region contains $`N`$ points.

Given a set $``$ of $`n`$ lines in a plane, we consider the set $`𝒟()`$ of their dual points: The dual point of a line $`y=axb`$ is $`(a,b)`$. We construct a range searching data structure for $`𝒟()`$. Given a point $`p=(x_0,y_0)`$, the set of dual points of lines below $`p`$ is the set of points in $`𝒟()`$ located below the line $`Y=x_0Xy_0`$, where $`X`$ and $`Y`$ correspond to coordinates of the dual plane. Thus, we can compute the level of the point $`p`$ in the arrangement of $`n`$ lines by using half-plane range searching. Moreover, we have the highest line below $`p`$ in the same query time. Also, we can query the number of lines which lie below both of a pair of query points.

A main building block for the range query is the point set decomposition structure of Matoušek , which we also need to utilize directly (we only describe its two-dimensional version):

###### Theorem 2.1 (Matoušek)

Given a set $`S`$ of $`n`$ points in the plane, for any given $`r<n`$, we can subdivide $`S`$ into $`r`$ disjoint subsets $`S_i`$ $`(i=1,2,\mathrm{},r)`$ such that $`|S_i|2n/r`$ satisfying the following condition: Each $`S_i`$ is enclosed in a triangle $`\sigma _i`$, and any line in the plane cuts at most $`cr^{1/2}`$ triangles among $`\sigma _1,\sigma _2,\mathrm{},\sigma _r`$ where $`c`$ is a constant independent of $`n`$ and $`r`$. Such a decomposition can be constructed in $`O(n\mathrm{log}n)`$ time.

## 3 Computing all peaks in $`k`$-level

We assume $`kn/2`$ for simplicity from now on; if $`k>n/2`$, replace $`k`$ by $`nk`$ and exchange maximal and minimal in the statements. A key observation for the $`k`$-level is that it is a subset of a union of $`k`$ concave chains such that all concave vertices of the $`k`$-level are vertices of these concave chains ; thus, a $`k`$-level has at most $`k`$ maximal peaks and $`k1`$ minimal peaks. We want to compute all the local peaks in a given interval $`I`$ of the $`x`$-coordinate value.

We apply a version of parametric search paradigm . However, before applying the parametric search, we start with a simpler “$`k`$-branching binary search” method. Without loss of generality, we assume that no line in the arrangement is horizontal nor vertical.

We prepare two key-subroutines: one-shot query and peak counting: Let $`𝒑(x_0)`$ be the point on the $`k`$-level at the $`x`$-coordinate value $`x_0`$. Let $`\mathrm{}_k^{}(x_0)`$ (resp. $`\mathrm{}_k^+(x_0)`$) be the line in the $`k`$-level at the $`x`$-coordinate value $`x_0ϵ`$ (resp. $`x_0+ϵ`$) for infinitesimally small $`ϵ>0`$. If $`x_0`$ is not an $`x`$-coordinate value of a vertex on the $`k`$-level, $`\mathrm{}_k^{}(x_0)=\mathrm{}_k^+(x_0)`$. The above operation to compute the point (together with lines containing the point) on $`_k`$ at a given $`x`$-coordinate value is called one-shot query. One-shot query is an analogue of ray shooting , and thus the following lemma is basically well-known. The complexity $`q(n,m)`$ given in the lemma is called one-shot query-time for the $`k`$-level:

###### Lemma 3.1

If we preprocess the lines in $``$ with $`O(m)`$ time for $`n\mathrm{log}n<m<n^2`$, given an $`x`$-coordinate value $`x_0`$, we can compute $`𝐩(x_0)`$, $`\mathrm{}_k^{}(x_0)`$ and $`\mathrm{}_k^+(x_0)`$ in $`q(n,m)=\stackrel{~}{O}(n/m^{1/2})`$ time, and also in polylogarithmic time by using $`O(n/m^{1/2})`$ processors.

Proof  By using the method given in the preliminary section, we can compute the level of any given point $`(x_0,y_0)`$ in polylogarithmic time by using $`O(n/m^{1/2})`$ processors. We now apply parametric searching to have the sequential time bound to compute the point $`𝒑(x_0)`$. A parametric searching algorithm is usually stated as a sequential algorithm; however, it is naturally a parallel algorithm if we use a parallel decision algorithm and also a parallel sorting algorithm. We remark that we can do it easier without using the parametric searching if we examine the range searching method in precise; however, we omit it in this paper.  $`\mathrm{}`$

The peak-counting is a routine to compute the number of peaks of the $`k`$-level in a given interval $`J=[x_0,x_1]`$ of $`x`$-coordinate values efficiently. The following elementary lemma is essential:

###### Lemma 3.2

Let $`f(x_0)`$ and $`f(x_1)`$ are numbers of positive-slope lines below or on the $`k`$-level at $`x_0`$ and $`x_1`$, respectively. Then, the number of maximal peaks of $`_k`$ in the interval $`J`$ is $`f(x_0)f(x_1)`$.

Proof  At-most-$`k`$-level (the part of the arrangement below $`k+1`$-level) is a union of $`k`$ concave chains such that all concave peaks in the chains appear in the $`k`$-level . If a concave chain among them has a peak in $`J`$, the slope of the chain must be changed from positive to negative. Thus, the number of maximal peaks within $`J`$ is the difference between the numbers of positive slope lines at two endpoints.  $`\mathrm{}`$

We remark that $`f(x_0)f(x_1)`$ equals the number of positive slope lines intersecting the segment between $`𝒑(x_0)`$ and $`𝒑(x_1)`$ if the segment has a nonpositive slope.

###### Lemma 3.3

For a given interval $`J`$ of $`x`$-coordinate value, the number $`\kappa (J)`$ of peaks of $`_k`$ in $`J`$ can be computed in $`O(q(n,m))`$ time if we preprocess the lines with $`O(m)`$ time. Also, the number of maximal peaks can be computed in
$`O(q(n,m))`$ time.

Proof  If we construct the dual of range search data structure for the set of lines with positive slopes, the number of positive slope lines below a given point $`(x_0,y_0)`$ can be computed in $`O(q(n,m))`$ time. Hence, $`f(x_0)`$ and $`f(x_1)`$ can be computed in $`O(q(n,m))`$ time, and the number of maximal peaks can be computed by using Lemma 3.2. The number of minimal peaks is easily computed from the number of maximal peaks and slopes of the $`k`$-level at endpoints.  $`\mathrm{}`$

Now, we can apply a binary search paradigm to design a weakly-polynomial time algorithm. First, for the input interval $`I`$, we compute the number of peaks $`\kappa =\kappa (I)`$ within the interval. The time complexity for this initialization is negligible, and obviously $`\kappa 2k1`$. Next, we construct a data structure for the one-shot query in $`O(m)`$ time, where the choice of $`m`$ will be explained later. Suppose that coefficients of the equations of lines are quotient numbers of $`\mathrm{log}\mathrm{\Gamma }`$ bit integers. We apply $`\kappa `$-branching binary search to find all peaks; At each stage of the binary search, we have at most $`\kappa `$ subintervals which has at least one local peak of $`_k`$ (such a subinterval is called an active interval), and we recursively search in active subintervals. Thus, after examining $`\kappa \mathrm{log}\mathrm{\Gamma }`$ candidates of $`x`$-coordinate values, we can find all of the peaks.

###### Proposition 3.4

All the local peaks of $`_k`$ in the interval $`I`$ can be computed in $`O(\kappa q(n,m)\mathrm{log}\mathrm{\Gamma })`$ time.

To make the complexity into strongly polynomial, we apply parametric search by using the parallel algorithm for the one-shot query given in Lemma 3.1 as its guide algorithm. We run the guide algorithm without inputting the parameter value (in our case, an $`x`$-coordinate value), and decide the $`x`$-coordinate values of the peaks by using the sequential one-shot query and the counting algorithm of Lemma 3.3 as decision algorithms. The counting algorithm needs two values of the parameter, which corresponds to endpoints of each of intervals obtained by splitting active intervals in the parametric searching process. Usually, parametric searching method solves optimization problems on monotone or convex functions. Here, $`k`$-level is neither monotone nor convex, but it consists of $`\kappa `$ monotone fragments. Thus, while running the guide algorithm, there are at most $`\kappa `$ different critical parameter values to determine all the comparisons in the current parallel step that are necessary to proceed into the next parallel step. In precise, the number of different choices is the number of active intervals obtained breaking $`I`$ by the critical parameter values found so far in the guide algorithm. We make a clone of the guide algorithm for each active interval. If the current interval is split into $`f`$ active subintervals, $`f1`$ new clones are created. Naturally, we create at most $`\kappa `$ clones in our process. There is only one critical parameter value to determine the comparisons in a “usual” parametric search, and such a value can be found if we run the decision algorithm $`O(\mathrm{log}N)`$ times if the guide algorithm is a parallel algorithm on $`N`$ processors. In our case, we run the decision algorithm $`O(\kappa \mathrm{log}N)`$ times at each level. Thus, we obtain the following theorem:

###### Theorem 3.5

All the peaks on $`_k`$ within an interval $`I`$ can be computed in $`O(n\mathrm{log}n)+\stackrel{~}{O}((\kappa n)^{2/3})`$ time if $`I`$ has $`\kappa `$ local peaks.

Proof  The parametric searching method gives $`\stackrel{~}{O}(\kappa q(n,m))`$ time complexity apart from the $`O(m)`$ preprocessing time. We balance $`\kappa q(n,m)=\stackrel{~}{O}(\kappa n/m^{1/2})`$ and $`m`$ to have $`m=(\kappa n)^{2/3}`$. If $`(\kappa n)^{2/3}<n`$, we instead use $`m=n\mathrm{log}n`$. This gives the time complexity.  $`\mathrm{}`$

Since $`\kappa 2k`$, we have the following:

###### Corollary 3.6

All the peaks on $`_k`$ can be computed in $`O(n\mathrm{log}n)+\stackrel{~}{O}((kn)^{2/3})`$ time.

### 3.1 Computing selected peaks

When $`\kappa `$ is large, it may be too expensive to compute all the local peaks. Suppose that we want to compute $`\tau `$ largest maximal peaks in the input interval $`I`$ for $`\tau \kappa `$ more efficiently than computing all the peaks. This can be done by combining Roos and Widmayer’s algorithm and the algorithm given above.

We first run a binary search process with respect to $`y`$-coordinate values similarly to Roos and Widmayer’s algorithm. At the intersection of the arrangement with a horizontal line $`y=y_0`$, we compute intervals on the line which are below the $`k`$-level (we call them semi-active intervals) in $`O(n\mathrm{log}n)`$ time. We next compute the sum $`s(y_0)`$ of numbers of local peaks in the semi-active intervals. We could apply our counting method of Lemma 3.3 for each semi-active intervals to compute the sum of maximal peaks in the intervals by using $`O(\tau q(n,m))`$ time. More simply, we can compute it in $`O(n\mathrm{log}n)`$ time by counting the number of intersecting positive slope lines with the horizontal line during the sweep. In precise, we also need information of the arrangement at endpoints of the input interval $`I`$ if one (or both) of them is below $`k`$-level (we omit details).

If $`s(y_0)`$ is not between $`\tau `$ and $`2\tau `$, we continue binary search on $`y_0`$: If $`s(y_0)>2\tau `$, we increase $`y_0`$ while if $`s(y_0)<\tau `$, we decrease $`y_0`$. Thus, we can eventually find a value $`y_0`$ such that $`\tau s(y_0)2\tau `$. We have spent $`O(n\mathrm{log}^2n)`$ time so far. Now, we search all peaks in the union of active intervals by using the method given in the previous section. The following theorem is easy to see:

###### Theorem 3.7

We can compute $`\tau `$ largest maximal peaks of $`_k`$ in an interval $`I`$ in $`O(n\mathrm{log}^2n)+\stackrel{~}{O}((\tau n)^{2/3})`$ time. We can also compute $`\tau `$ largest local peaks (including both maximal and minimal peaks) in the same time complexity.

Note that if we only use Roos and Widmayer’s algorithm in a naive fashion to find $`\tau `$ largest peaks, it would cost $`O(\tau n\mathrm{log}^2n)`$ time. Analogously, we can compute the $`\tau `$-smallest minimal peaks.

###### Theorem 3.8

We can compute $`\tau `$ smallest minimal peaks of $`_k`$ in an interval $`I`$ in $`O(n\mathrm{log}^2n)+\stackrel{~}{O}((\tau n)^{2/3})`$ time. We can also compute $`\tau `$ smallest local peaks (including both maximal and minimal peaks) in the same time complexity.

For this $`\tau `$ smallest peak finding problem, we can modify Roos and Widmayer’s method  to compute them in $`O(n\mathrm{log}^2n+k\tau \mathrm{log}^2k)`$ time. Although $`n+(\tau n)^{2/3}2(n+k\tau )`$ always holds, the time complexity is better than the one in Theorem 3.8 by a polylogarithmic factor if $`n^{1/2}\mathrm{log}^cn<k<n^{1/2}\mathrm{log}^cn`$ for some constant $`c`$. Moreover, the algorithm does not use range search data structure, and hence much simpler and uses smaller data space.

###### Proposition 3.9

The $`\tau `$ smallest local peaks in an interval $`I`$ can be found in $`O(n\mathrm{log}^2n+k\tau \mathrm{log}^2k)`$ time.

Proof  First we search for a horizontal line $`h`$ such that it intersects $`k`$-level, and the number of peaks below it is not less than $`\tau `$. Such a line $`h`$ can be found in $`O(n\mathrm{log}^2n)`$ time. Next, let $`v`$ and $`w`$ be the leftmost intersection and the rightmost intersection with the $`k`$-level on $`h`$, and let $`J`$ be the interval between them. If the $`k`$-level at an endpoint of the input interval $`I`$ is below $`h`$, we connect the point on $`k`$-level at the $`x`$-value of the endpoint to $`J`$ with a segment to form a chain $`C`$ with at most three segments. Let $`_0`$ be the set of lines in the arrangement intersecting with the chain $`C`$. The cardinality of $`_0`$ is at most $`2k`$ because of Lemma 3.2. Finally, we find all peaks below $`h`$ by using the $`\tau `$-branching binary search. Here, by using the windowing method of , we should only take care of lines in $`_0`$ together with lines below endpoints of the chain $`C`$. There are at most $`4k`$ such lines. Hence, this second step can be done in $`O(k\tau \mathrm{log}^2k)`$ time.  $`\mathrm{}`$

## 4 Bottleneck edge length in a <br>parametric spanning tree

Next, we consider the parametric spanning tree problem. Consider a connected graph $`G=(V,E)`$ with $`k`$ nodes and $`n`$ edges. Because of the connectivity, $`k1nk(k+1)/2`$. For each edge $`eE`$, we associate a weight function $`w_e(x)`$, which is linear on a parameter $`x`$. We assume that the arrangement generated by lines $`y=w_e(x):eE`$ is simple, i.e., no three lines intersect at a point. We can remove this assumption by giving a symbolic perturbation. $`G`$ is denoted by $`G(x)`$ if it is considered as a weighted graph with parametric weights. For a given value $`x`$, we consider the minimum spanning tree $`T(x)`$ of $`G(x)`$.

It is known that the number of transitions of the structure of the minimum spanning tree $`T(x)`$ is $`O(k^{1/3}n)`$, and all the transitions can be computed in $`O(kn\mathrm{log}n)`$ time . Moreover, the average edge weight in the minimum spanning tree is a concave function in $`x`$, and the value of $`x`$ maximizing the average edge weight of $`T(x)`$ can be computed in $`O(n\mathrm{log}n)`$ time .

As parametric matroid problems, the average edge weight is a counterpart of the average of $`y`$-values of $`k`$ lines below (or on) the $`k`$-level. A natural counterpart of the $`k`$-level itself in the minimum spanning tree is the longest (i.e. maximum weight) edge in the minimum spanning tree. The edge is also called the spanning bottleneck edge at $`x`$ ($`SBE(x)`$ in short), and its weight is denoted by $`w_{SBE}(x)`$. It is easily observed that $`W_{SBE}(x)`$ is the minimum value of $`w`$ such that the subgraph of $`G(x)`$ constructed from the set of edges whose weights are less than or equal to $`w`$ is connected.

It is natural and important problem in sensitivity analysis  to trace the trajectory $`y=w_{SBE}(x)`$ of the weight of $`SBE(x)`$. Analogously to the $`k`$-level, there are at most $`k`$ maximal peaks in the trajectory $`y=w_{SBE}(x)`$. We want to compute peaks in the trajectory.

### 4.1 One-shot query for the longest edge in MST

We first consider efficient query for $`SBE(x_0)`$ at any given value of $`x_0`$ of the parameter. This query is called one-shot query for the $`SBE`$. A naive method is the following: First construct $`T(x_0)`$ in $`O(n)`$ time, and select its longest edge. Instead, we use the Matoušek’s set partition. In the dual space, the dual points of $`n`$ weight functions of the edges are partitioned into $`r`$ subsets of size $`O(n/r)`$. Each subset is contained in a triangle, and $`O(r^{1/2})`$ triangles are cut by any query line.

Accordingly, we partition the set of $`n`$ edges into $`r`$ subsets each has $`O(n/r)`$ edges. For each subset, we compute a spanning forest (irrelevant to edge weights) and store the connected components except singletons. Thus, each component has a forest with $`O(\mathrm{min}\{k,n/r\})`$ edges. This computation can be done in $`O(n)`$ additional time.

If we are given a parameter value $`x_0`$, we sort $`O(r)`$ vertices of the triangles with respect to the inner product of them with the vector $`(x_0,1)`$. We do binary search on this sorting list. We guess a vertex $`v`$, and consider a line $`\mathrm{}:Y=x_0X+c`$ which goes through $`v`$. We recognize the triangles which are below $`\mathrm{}`$; thus, the edges in the subsets associated with the triangles has weights less than $`c`$. We construct a spanning forest $`F`$ of the union of forests in these subsets: since they have $`O(rk)`$ edges, this can be done in $`O(rk)`$ time. If the forest $`F`$ is a spanning tree, we decide that $`v`$ is too large in the sorting list, and continue the binary search.

Otherwise, we consider the subsets associated with the triangles cut by $`\mathrm{}`$. They contain $`O(n/r^{1/2})`$ points in total. We sort them with respect to the weights, and greedily insert them into $`F`$ until we have a spanning tree. If we do not have a spanning tree, we decide $`v`$ is too small, and continue the binary search. If we have a spanning tree, we decide $`v`$ is a candidate, but it may be too large, and continue to search for the lowest vertex $`v`$ satisfying the above condition, and return the longest edge in the tree for that $`v`$. Note that the spanning tree is not a minimum spanning tree in general; however, we correctly recognize the longest edge in a minimum spanning tree.

This process needs $`O(rk+n/r^{1/2})`$ time, and we do this process $`O(\mathrm{log}r)`$ times during the binary search. Thus, if we set $`r=(n/k)^{2/3}`$, the time complexity is $`O(n^{2/3}k^{1/3}\mathrm{log}(n/k))`$, which is slightly better than $`O(n)`$ if $`k=o(n)`$. By applying a hierarchical subdivision, we can further improve it: We fist start $`r=r_1`$, and decompose the subset of size $`O(r_1n)`$ into $`r_2`$ smaller subsets, where $`r_2=r_1^{1/2}`$, and we further continue the refinement for $`r_i=r_{i1}^{1/2}`$ until $`r_i`$ becomes below a constant. The query time becomes
$`k(r_1+r_1^{1/2}r_2+\mathrm{}+(r_1r_2\mathrm{}r_{i1})^{1/2}r_i)+n/(r_1r_2\mathrm{}r_i)^{1/2}`$. Setting $`r_1=(n/k)^{1/2}`$, this enables $`\stackrel{~}{O}((nk)^{1/2})`$ time computation for $`SBE(x_0)`$. Similarly to the case of halfspace range searching, we can combine hierarchical cutting  of the arrangement to have a preprocess-query trade-off (we omit details in this version). Indeed, we have the following proposition:

###### Proposition 4.1

If we spend $`\stackrel{~}{O}(m)`$ preprocessing time for $`nm`$, we can do the one-shot query for $`SBE`$ in
$`\stackrel{~}{O}(n/(m/k)^{1/2}+k)`$ time.

Moreover, we will later use the following two-shot reporting query for a spanning forest, which reports a spanning forest consisting of edges whose weight functions are below both of given two query points $`(x_0,y_0)`$ and $`(x_1,y_1)`$. This can be done similarly to one-shot query (this is a counterpart of the simplex range searching if the one-shot query is a counterpart of the halfplane range searching).

###### Proposition 4.2

If we spend $`\stackrel{~}{O}(m)`$ preprocessing time for $`nm`$, we can do the two-shot reporting query in
$`\stackrel{~}{O}(n/(m/k)^{1/2}+k)`$ time.

### 4.2 Computing the maximum peak

Let us consider the problem of computing the maximum peak in $`I`$. First, we straightforwardly apply Roos-Widmayer’s algorithm. For a given $`y`$-value $`y_0`$, we want to decide whether $`Max_{xI}w_{SBE}(x)y_0`$ or not. We dynamically update the spanning forest associated with edges with weight below $`y_0`$ from $`x=x_0`$ to $`x=x_1`$ if $`I=[x_0,x_1]`$. If we find a value $`xI`$ such that the spanning forest becomes a spanning tree, we know $`Max_{xI}w_{SBE}(x)y_0`$. It costs $`O(k^{1/2})`$ time to update a minimum spanning forest for insertion and deletion of edges. Suppose that we sweep on the line $`y=y_0`$ updating the minimum spanning forest. The line $`y=y_0`$ has at most $`n`$ intersections with lines associated with weight functions, and hence the method needs $`O(nk^{1/2})`$ time for the decision. Thus, the maximum peak can be found in $`O(nk^{1/2}\mathrm{log}n)`$ time.

We try to improve the above time complexity. We subdivide the line $`y=y_0`$ into $`n/s`$ intervals such that each interval contains at most $`s`$ intersection points. For each interval $`I_i=[x_i,x_{i+1}]`$, we perform the two-shot reporting query at $`(x_i,y_0)`$ and $`(x_{i+1},y_0)`$. The reported forest $`F`$ is constructed from edges whose weight is less than $`y_0`$ both at $`x=x_i`$ and $`x=x_{i+1}`$. If the forest has more than $`s+1`$ connected components, it is impossible that $`w_{SBE}(t)y_0`$ for a $`tI_i`$. Otherwise, we dynamically maintain the spanning tree, where we contract nodes into at most $`s+1`$ super nodes each of which associate with a connected component of the forest $`F`$. Our graph has only $`s`$ edges, and hence the update can be done $`O(\sqrt{s})`$ time per intersection.

Hence, total time complexity becomes
$`\stackrel{~}{O}(n\sqrt{s}+(n/s)[n/(m/k)^{1/2}+k]+m)`$.
If we optimize this, we have $`\stackrel{~}{O}(n^{8/7}k^{1/7}+nk^{1/3})`$. This is an improvement over $`O(nk^{1/2})`$, since $`nk(k+1)/2`$.

The minimum of $`w_{SBE}(x)`$ can be analogously computed. Hence, we have the following theorem:

###### Theorem 4.3

For a given interval $`I`$ of the parameter value $`x`$, we can compute both the maximum and the minimum of $`w_{SBE}(x)`$ for $`xI`$ in $`\stackrel{~}{O}(n^{8/7}k^{1/7}+nk^{1/3})`$ time.

We can generalize the above theorem for the truncated matroid of the graphic matroid to obtain the following proposition (we omit the proof):

###### Proposition 4.4

For a constant $`c`$, let $`w_{SBEc}(x)`$ be the minimum value of $`w`$ such that the set of edges in $`G(x)`$ with weights less than or equal to $`w`$ has at most $`c`$ connected components. We can compute both the maximum and the minimum of $`w_{SBEc}(x)`$ in $`\stackrel{~}{O}(n^{8/7}k^{1/7}+nk^{1/3})`$ time.

### 4.3 Computing all local peaks for $`SBE`$

It is desired to apply the method of computing all peaks in the $`k`$-level to $`SBE`$ in a parametric graph. It is known that there are at most $`k`$ maximal peaks in $`y=w_{SBE}(x)`$. Unfortunately, we have no good method to know the number of peaks within an interval $`I=[x_0,x_1]`$ exactly, since we do not have a property that is a counterpart of Lemma 3.2. The only known method for the authors is to compute the minimum spanning trees $`T(x_0)`$ and $`T(x_1)`$ explicitly, and compute the difference $`d(I)`$ between the number of edges whose weight functions are positive slopes. For any disjoint set of intervals, the sum of $`d(I)`$ over the intervals is at most $`k`$, and $`d(I)`$ gives an upper bound of maximal peaks of $`y=w_{SBE}(x)`$ within $`I`$. However, it is often an overestimate, since $`d(I)`$ gives the number of maximal peaks of $`k`$ trajectories of weights of all edges (not only maximum one) in the parametric minimum spanning tree, where we only include peaks where an edge in the spanning tree is replaced by an edge outside the spanning tree.

By using the above observation, we have an $`\stackrel{~}{O}(kf(n,k))`$ time algorithm if we have an algorithm to compute $`T(x_0)`$ for a given $`x_0`$ in $`O(f(n,k))`$ time. Unfortunately, we only have $`f(n,k)=O(n)`$, which leads to an $`\stackrel{~}{O}(kn)`$ time algorithm, which is inferior to a known algorithm to compute all transitions of the parametric minimum spanning tree. We can compute $`\tau `$ largest peaks in transitions of edge weights in MST in $`\stackrel{~}{O}(\tau n)`$ time, if we include all peaks of all edges in MST; however, the number of peaks appeared at the transitions of the longest edge among them may be much smaller than $`\tau `$.

Although the above method is not attractive for the minimum spanning tree, the method is applicable to any parametric matroid, and hence it is useful if we do not have a dynamic algorithm to maintain a minimum weight base. The current $`O(kn\mathrm{log}n)`$ time algorithm to compute the transitions of parametric minimum spanning tree needs $`O(k^{2/3})`$ time method (indeed, it can be done in $`O(k^{1/2})`$ time) to update a minimum spanning tree. Thus, the analogue of the algorithm needs $`\stackrel{~}{O}(k^{1/3}nq+kn)`$ time to compute all the transitions of a parametric matroid of rank $`k`$, where $`q`$ is the time complexity to update its minimum weight base.

## 5 Concluding remarks

The number of maximal peaks in a $`k`$-level is known to be at most $`{}_{k}{}^{}C_{d}^{}`$ (number of combinations choosing $`d`$ elements from $`k`$ elements) if we have $`d`$ dimensions . Hence, this is much smaller than the complexity of whole arrangement, especially if $`k`$ is much smaller than $`n`$. However, to the author’s knowledge, the problem of computing peaks in the $`k`$-level for a higher dimensional arrangement is open. Although the algorithm of Roos and Widmayer can be applied to $`3`$-dimensional case, it needs $`\stackrel{~}{O}(n^2)`$ time to compute the largest peak (i.e. global maximum) if we naively implement it. One necessary constituent is to develop a counterpart of Lemma 3.3: Given an arrangement of $`n`$ hyperplanes in the three-dimensional space, preprocess it, and for any given three points $`A`$, $`B`$, and $`C`$ in the plane $`z=0`$, decide whether the triangle $`ABC`$ contains (a projection) of a peak in the $`k`$-level or not efficiently. For the purpose, we probably need a counterpart of Lemma 3.2: Give a criterion of the existence of a peak from the information of the set of hyperplanes below $`k`$-level at each of $`A`$, $`B`$, and $`C`$. In two-dimensional space, the lines are classified into positive slope lines and negative slope lines, while this kind of natural discrete classification of planes in the space does not exist. Moreover, it is difficult to find a counterpart of concave chain decomposition for three dimensional $`k`$-level . These lacks make the problem difficult, although the authors think it is quite attractive.

Another interesting problem is an an extension of parametric $`SBE`$ problem to the $`c`$-edge-connectedness for $`c2`$. Here, we hope we can develop efficient solutions by combining geomertic methods and graph theoretical methods.

### Acknowledgement

The authors gratefully acknowledge referees of the SOCG conference for informing several references on recent improved algorithms for computing $`k`$-level and parametric minimum spanning trees.


text:
# A one-dimensional toy model of globular clusters

## 1 Introduction

Globular clusters are gravitationally bound concentrations of large numbers of stars, spherically distributed in space. They orbit around a galaxy spending most of the time in the galactic halo . The most important elements governing globular clusters structure are two-body relaxation and truncation due to tidal forces. Different dynamical models considering these specific phenomena, have been investigated both analytically  and numerically . For an extended discussion, see  and Refs. therein.

Dynamical evolution causes stars to escape as an effect of the gravitational interaction with the nearby galaxies. This evaporation process drives the cluster towards a configuration with a high density core and the velocity dispersion of stars in the bulk can increase without limit. This phenomenon is known as gravothermal catastrophe and its study goes back to Antonov  and to Lynden-Bell and Wood .

Referring to the pioneering work of Chandrasekhar , it is possible to calculate the perturbations induced by stellar encounters on star motion. This is done by means of a diffusion model, which allows to derive a quantitative description of changes of star velocities in terms of single encounters. Considering weak encounters, i.e. solving the diffusion model in the Fokker-Planck approximation, King  found the following expression for the velocity profile

$$\left\{\begin{array}{cc}f_0(v)=A\left(\mathrm{exp}\left(\frac{v^2}{2\sigma ^2}\right)\mathrm{exp}\left(\frac{v_c^2}{2\sigma ^2}\right)\right)\hfill & \text{for }vv_c\hfill \\ & \\ f_0=0\hfill & \text{for }v>v_c\text{,}\hfill \end{array}\right\}$$
(1)
where $`v_c`$ is a cutoff velocity of the stars, $`\sigma `$ is the one-dimensional velocity dispersion and $`A`$ a normalization constant. King’s models have been shown to be in agreement with the observed brightness surface profiles of globular clusters . Further developments of King’s model  consider the same functional form (truncated Gaussian) applied to the gravitational energy, and hence propose a general form for the distribution function $`f(𝐱,𝐯)`$. We will not consider here these extensions.

In this paper we discuss a simplified one-dimensional $`N`$-body model which reproduces King’s distribution. We let the particles, all of equal mass $`m`$, interact through the one-dimensional gravitational potential $`V=2\pi Gm^2|x|`$, where $`G`$ is the gravitational constant. Bearing in mind the comparison with globular clusters, we imitate the effect of galactic tidal forces by introducing a finite cut-off in positions. Thus, the evaporation of stars from the system is the only “dissipative” effect we consider. This is enough to drive the system towards an asymptotic non stationary regime, which we analyze in detail, and which reveals striking similarity with King’s model.

The main difference between the simplified one-dimensional model considered here and a more realistic three-dimensional one, is the lack of any singularity of the potential at the origin. The presence of a finite lower bound for the $`1D`$ potential makes less energy available to support the evaporation process as the system cools down. In the $`3D`$ case an infinite amount of energy can indeed be extracted from the singular pair-wise gravitational interaction, which is the main origin of the gravothermal catastrophe. This is the reason why the model we discuss cannot reproduce the core collapse corresponding to the gravothermal catastrophe.

In the next Section we present the model and the results concerning velocity distribution and density profiles. Section 3 is devoted to the discussion of power-laws for the number of particles in the cluster. Finally, in Section 4 we draw some conclusions.

## 2 One-dimensional model

Let us consider a one-dimensional classical Newtonian self-gravitating system of $`N`$ particles with equal mass $`m`$, with Hamiltonian

$$H=\frac{1}{2}m\underset{i}{\overset{N}{}}v_i^2+2\pi Gm^2\underset{j<i}{}|x_ix_j|,$$
(2)
where $`x_i`$ is the position and $`v_i`$ the velocity of $`i`$-th star. $`G`$ is the universal gravitational constant. We choose in the following $`m=1`$ and $`2\pi G=1`$. This system has been recently the subject of intensive investigations . Particle accelerations are constant inbetween two collisions and are proportional to the net difference of particles respectively on the right and on the left. When a collision occurs particles cross each other, or, equivalently, collisions are elastic. This particle approach is known to corresponds in the continuum limit ($`N\mathrm{}`$) to the Vlasov-Poisson equations for the distribution function $`f(x,v)`$

$`{\displaystyle \frac{f}{t}}+v{\displaystyle \frac{f}{x}}{\displaystyle \frac{V}{x}}{\displaystyle \frac{f}{v}}=0`$
$`{\displaystyle \frac{^2V}{x^2}}=4\pi Gm{\displaystyle f(x,v)𝑑v},`$ (3)
where $`V(x)`$ is the self-consistent gravitational potential.

We add the following features to the model

* Particles are confined in a box of size $`L`$, i.e. $`x_i[L/2,L/2]`$
* The effect of tidal forces induced by the parent galaxy is imitated by requiring that each time a particle reaches the boundary of the box with a finite velocity, it drops out of the system, which therefore experiences a mass and energy loss (“evaporation”). This last feature implies that system (3) is solved with absorbing boundary conditions.

The numerical implementation is based on an “event-driven” scheme, first introduced in plasma physics , which is adapted to the present case as follows. The algorithm looks for the particles which collide the first and for the time when the event occurs, $`t_{coll}`$. Then it computes the first “evaporation time”, $`t_{evap}`$, and makes the system evolve until the minimum, $`t_{min}`$, between the collision and the evaporation time is reached. Once the system experiences evaporation, the total mass is reduced and the escaping particle stops interacting with the residual bulk. By rescaling the position and the velocity of each particle, i.e. introducing a local dissipation, we maintain the position of the center of mass fixed and its velocity to zero. This means we simply translate velocity and position of the remaining particles to keep the system centered in position and momentum space. A particle can escape from the system as a result of this rescaling: this possibility is taken into account even if it is has a low probability.

Evaporation is a singular event, which, in fact, marks the transition between two self-gravitating system having a different number of particles and energy. We remark that the integration scheme is “exact”. Time $`t`$ elapsed from the initial configuration is obtained by summing all values of $`t_{min}`$ up to the last event.

In all numerical experiments, $`N(0)`$ particles are initially uniformly distributed in the box and the initial velocities are Gaussian i.i.d random variables, with the temperature $`T_0`$ given by twice the average kinetic energy.

In Fig. 1 we show the time evolution of the number of particles $`N(t)`$ which remain inside the box up to time $`t`$. After an abrupt decrease of $`N(t)`$, which strongly depends on the initial condition, the system reaches a state where rare evaporations are present, making $`N(t)`$ decrease much slowlier. Our numerical experiments show that such a quasi-stationary state lives indefinitely, although we cannot exclude that, finally, $`N(t)`$ relaxes to an asymptotic value $`N_{as}`$. As the best approximation for this value, we take the one the system reaches in the longest computer runs.

Given $`L`$, for large enough $`T_0`$ the system approaches a state characterized by a single cluster, which adapts itself to the size of the box. In Fig. 2 we show the phase-space portrait for a system in its late stage evolution, when the most energetic particles have dropped out from the box and the system has relaxed to an asymptotic plateau, as the ones reported in Fig. 1. The particles are almost uniformly distributed within a bounded region of the phase plane. Both this fact and the shape of the contour suggest a possible connection with the so called water-bag (WB) distribution . This is a stationary solution of the Vlasov-Poisson system (3), $`f^{WB}(x,v)`$, which is constant in a simply connected domain $`\mathrm{\Omega }`$ of the phase-plane and strictly zero outside.

Adopting the notation of Ref. , we call respectively $`x_s`$ and $`v_s`$ the maximum position and velocity of the WB. The potential $`V(x)`$ for such a distribution is then implicitly specified by the following integral equation

$$x=\frac{ϵ^{\frac{3}{4}}}{N_{as}}_0^{V(x)}\left[ϵ^{\frac{3}{2}}(ϵ\zeta )^{\frac{3}{2}}\right]^{\frac{1}{2}}𝑑\zeta .$$
(4)
The maximal energy of the water-bag, $`ϵ`$, is such that, if we express $`f(x,v)`$ in terms of the energy $`u=v^2/2+V(x)`$, $`f(x,v)=F(u)=0`$ if $`u>ϵ`$. The energy $`ϵ`$ is related to $`x_s`$ and $`v_s`$ by $`ϵ=v_s^2/2`$, $`V(x_s)=ϵ`$ and the zero energy level is fixed by requiring that $`V(0)=0`$. The density profile, $`\rho ^{WB}(x)`$, is expressed as function of the potential

$$\rho ^{WB}(x)=\frac{3N_{as}^2}{8ϵ^{3/2}}(ϵV)^{1/2}.$$
(5)
Since the distribution function $`f(x,v)`$ is constant over $`\mathrm{\Omega }`$, it follows immediately

$$\rho ^{WB}(x)=_{v^+(x)}^{v^+(x)}f(x,v)𝑑v=2cv^+(x)$$
(6)
where $`v^+(x)`$ represents the profile of the upper branch of the WB contour, which we assumed to be symmetric, and $`c=3N_{as}^2/(16\sqrt{2}ϵ^{3/2})`$. Using Eq. (5), this implies

$$v^+(x)=\sqrt{2(ϵV)}.$$
(7)
To compute the velocity contour $`v^+(x)`$, we need to know $`V(x)`$, which we do by solving Eq. (4) by an adaptive recursive Newton-Cotes 8 panel rule with tolerance $`10^7`$. This velocity contour is drawn in Fig. 2 and, as predicted, encloses all phase points. In deriving $`V(x)`$ we have taken the value $`v_s`$ from the cluster phase plot; this is the only phenomenological input in this calculation and the agreement with the data has to be considered quite satisfactory.

Moreover, we can compare the theoretically derived potential $`V(x)`$ with the one computed directly from the asymptotic positions

$$V(x_i)=\underset{j}{}|x_ix_j|.$$
(8)
We need, of course, to perform a vertical shift to fix the zero in the origin. The result of formula (8) is reported in Fig. 3 together with the theoretically derived potential. The agreement is really good. We are thus led to conclude that our asymptotic state is well described by a water-bag. However, this latter is a stationary solution of the self-gravitating 1D system, while in our simulations we continue to observe some particle evaporations even at very long times. This is why we have called our asymptotic state quasi-stationary and its description in terms of a water-bag distribution can only be approximate.

An alternative treatment of the asymptotic state is based on King’s formula (1). In this case one does not try to reproduce the full distribution function, but just its projections along the $`x`$ and $`v`$ axis: $`\rho (x)=f𝑑v`$ and $`f_0(v)=f𝑑x`$. Following the standard derivation of the equilibrium isothermal distribution , we are led to introduce an analytical ansätz for the density profile

$$\left\{\begin{array}{cc}\rho (x)=A(cosh^2(Bx)cosh^2(B\frac{L}{2}))\hfill & \text{for }|x|<\frac{L}{2}\hfill \\ & \\ & \\ \rho (x)=0\hfill & \text{for }|x|>\frac{L}{2}\text{ ,}\hfill \end{array}\right\}$$
(9)
where the normalization $`A`$ is fixed by $`_{L/2}^{L/2}\rho (x)𝑑x=N_{as}`$. For the velocities we take King’s distribution, specified in Eq. (1). Assuming $`\rho (x)`$ as in Eq. (9) we can derive a close analytical expression for the potential. For a one-dimensional system, the following relation holds in the continuum limit

$$V(x)=_{\frac{L}{2}}^{\frac{L}{2}}|yx|\rho (y)𝑑y.$$
(10)
Inserting Eq. (9) into Eq. (10) and performing the integral, we get

$$V(x)=V_0\left[\frac{L}{B}\mathrm{tanh}(B\frac{L}{2})(x^2+\frac{L^2}{4})\mathrm{cosh}^2(B\frac{L}{2})+\frac{2}{B^2}\mathrm{ln}(\frac{\mathrm{cosh}(Bx)}{\mathrm{cosh}(B\frac{L}{2})})\right],$$
(11)
with

$$V_0=\frac{N}{\frac{2}{B}\mathrm{tanh}(B\frac{L}{2})L\mathrm{cosh}^2(B\frac{L}{2})}.$$
(12)
$`V(x)`$ is quadratic for small $`x`$. To verify the reliability of our guess we re-analyze the data previously discussed in connection with the water-bag distribution. In Fig. 4 we plot the potential calculated numerically from Eq. (8) together with a one-parameter fit, using Eq. (11). Again the agreement with the data is very good and even apparently superior to the one obtained using the WB picture. This is simply due to the fact that here we perform a one parameter fit, while, in the previous discussion, $`v_s`$ was arbitrarily deduced from the phase space analysis. As a cross-check, we introduce in Eq. (9) the coefficient $`B`$ determined from the fit of the potential. The resulting density profile is plotted in Fig. 5 and it agrees with the normalized histogram of particles position. Finally, an histogram of the velocity is represented in Fig. 6. The reverse-cup shape due to the cut-off of the tails is evident. The solid line in Fig. 6 is a numerical fit which uses the expression of Eq. (1) with $`v_c`$ and $`\sigma `$ as free parameters.

As a side remark we observe that, coherently with the observed form of the potential, each particle oscillates almost harmonically inside the box. This can be seen by looking at the asymptotic orbit of a single particle (Fig. 7). The slight diffusion of the orbit is the signature of the interaction with the other particles, which induces a weak chaoticity.

A further aspect that we have tested is the dependence of the dynamics on the initial temperature, $`T_0`$. Indeed, for small values of $`T_0`$, the system shows a pronounced collapse, which leads to a massive central core, as it is clearly displayed in the main plot of Fig. 8. This phase-space distribution significantly differs from the one in Fig. 2 and cannot be represented by a water-bag. The histograms of positions and velocities are computed and plotted in the right and left inserts, respectively. Both the density and the velocity profiles are very well reproduced by a numerical fit based on our ansätz (9) and on King’s distribution (1).

In conclusion, the ansätz we have introduced shows a good agreement with numerical data for all values of $`T_0`$ we have simulated, while the water-bag distribution fails to reproduce the velocity and density profiles at very low initial temperature. However, in the high temperature range, the water-bag treatment is superior, because it leads to an accurate description of the full distribution $`f(x,v)`$.

## 3 Scaling laws

In this Section we discuss some numerically found scaling laws which do not have presently a theoretical justification but which are an important signature of the presence of a finite box.

We follow the system until it reaches the asymptotic state. Being $`N_{coll}`$ the number of collisions, we define the “average collision time”, $`\tau =\frac{t}{N_{coll}}`$.

In the main plot of Fig 9 we represent $`N_{as}`$ as a function of $`\tau `$. Each point refers to a different value of the initial temperature $`T_0`$, varying form $`0.2`$ to $`7`$, while $`N_{coll}`$ is maintained constant for each realization. The initial temperature controls the rate of evaporation at a very early stage of the evolution. Larger values of $`T_0`$, produce higher mass loss, inducing the system to relax to a quasi-stationary state characterized by less residual particles $`N_{as}`$ (see Fig. 1). This process has consequences at the dynamical level, determining a larger mean free path and consequently a larger value of $`\tau `$. The curve in Fig 9 is consistent with this qualitative picture, showing a power-law decay with exponents $`a=0.4`$ which is valid over more than two decades. In the insert of Fig. 9, we plot $`\tau `$ vs. $`T_0`$.

We have also checked the dependence on $`N`$ of the collision time $`t_{coll}`$, defined as an average of all the collision times corresponding to each fixed value of $`N`$, i.e. inbetween two successive evaporations. In Fig. 10 we plot the results of different numerical experiments where we vary the initial temperature $`T_0`$. The dependence of $`t_{coll}`$ on $`N`$ is again a power-law with exponent $`b=2.5`$. Note that $`b1/a`$, hence Fig. 9 can be thought as a macroscopic averaged image of the microscopic properties shown in Fig. 10.

## 4 Conclusions

We have introduced a one-dimensional toy model of globular clusters with an emphasis on the evaporation process. With this in mind, we have discussed the effect of introducing a finite size box in a classical one dimensional self-gravitating medium. The dynamics of the system has been investigated for a special class of initial conditions. We pointed out the appearance of two distinct, non stationary, asymptotic regimes which occur depending on the temperature of the initial realization. For small values of $`T_0`$, similarities with the isothermal solution are found while for larger temperatures the density and velocity profiles are well reproduced also assuming a water-bag distribution.

We propose a form of the density profile, with a cut-off in the tails, which fits well numerical data in all the explored regimes, allowing to derive a close analytical expression of the gravitational potential. Moreover, a King-like velocity profile is shown to be in good agreement with numerical data. The asymptotic truncated profiles are thus a direct consequence of the evaporation from the finite box.

Finally we have also given numerical evidence of some scaling laws, which remain to be theoretically explained, but which are strongly related to the escaping process.

In the future, we plan to extend this study to the system of concentric spherical mass shells  by introducing an external absorbing boundary in the configuration space as done here.

Acknowledgements We thank an anonymous referee for his helpful comments, which led to a substantial improvement of this paper. We acknowledge useful discussions with E. Aurell, M.C. Firpo, M. Henon and A. Noullez. This work is part of the MURST-COFIN00 grant on Chaos and localization in quantum and classical mechanics. D.F. warmly acknowledges the research group Dynamics of Complex Systems in Florence for the kind hospitality.