text:
# Beam Simulation Tools for GEANT4 (and Neutrino Source Applications)

## 1 Introduction

Geant4 is a tool kit developed by a collaboration of physicists and computer professionals in the High Energy Physics (HEP) field for simulation of the passage of particles through matter. The motivation for the development of the Beam Tools is to extend the Geant4 applications to accelerator physics. The Beam Tools are a set of C++ classes designed to facilitate the simulation of accelerator elements such as r.f. cavities, magnets, absorbers. These elements are constructed from the standard Geant4 solid volumes such as boxes, tubes, trapezoids, or spheres.

A variety of visualization packages are available within the Geant4 framework to produce an image of the simulated apparatus. The pictures shown in this article were created with Open Inventor , which allows direct manipulation of the objects on the screen, plus perspective rendering via the use of light.

Although there are many computer programs for beam physics simulations, Geant4 is ideal to model a beam through a material or to integrate a beam line with a complex detector. There are many such examples in the current international High Energy Physics programs.

## 2 A Brief Introduction to Geant4

Geant4 is the object oriented C++ version of the Geant3 tool kit for detector simulation developed at CERN. It is currently being used in many fields, such us HEP, space exploration, and medicine.

As a tool kit, Geant4 provides a set of libraries, a `main` function, and a family of initialization and action classes to be implemented by the user. These classes are singlets, and their associated objects are constructed in `main`. The objects contain the information related to the geometry of the apparatus, the fields, the beam, and actions taken by the user at different times during the simulation. The Geant4 library classes start with the `G4` prefix. The example described in this section, called `MuCool`, uses only some of the many available user classes.

### 2.1 Detector and Field Construction

The detector and field geometry, properties, and location are implemented in the constructor and methods of the `MuCoolConstruct` user class, which inherits from `G4VUserDetectorConstruction`. In the `Construct()` method the user does the initialization of the electromagnetic field and the equation of motion. There are a variety of Runge-Kutta steppers to select from, which perform the integration to different levels of accuracy. Next comes the detector description, which involves the construction of solid, logical, and physical volume objects. They contain information about the detector geometry, properties, and position, respectively. Many solid types, or shapes, are available. For example, cubic (box) or cylindric shapes (tube), are constructed as:

```
G4Box(const G4String& pName, G4double pX, G4double pY, G4double pZ);
G4Tubs(const G4String& pName, G4double pRMin, G4double pRMax,
       G4double pDz, G4double pSPhi, G4double pDPhi);

```

where a name and half side lengths are provided for the box. Inner, outer radii, half length, and azimuthal coverage are the arguments of a cylinder (tube). A logical volume is constructed from a pointer to a solid, and a given material:

```
G4LogicalVolume(G4VSolid* pSolid, G4Material* pMaterial,
                const G4String& name)

```

The physical volume, or placed version of the detector is constructed as:

```
G4PVPlacement(G4RotationMatrix *pRot, const G4ThreeVector &tlate,
              const G4String& pName, G4LogicalVolume *pLogical,
              G4VPhysicalVolume *pMother, G4bool pMany, G4int pCopyNo);

```

where the rotation and translation are performed with respect to the center of its “mother” volume (container). Pointers to the associated logical volume, and the copy number complete the list of arguments.

### 2.2 Physics Processes

Geant4 allows the user to select among a variety of physics processes which may occur during the interaction of the incident particles with the material of the simulated apparatus. There are electromagnetic, hadronic and other interactions available like: “electromagnetic”, “hadronic”, “transportation”, “decay”, “optical”, “photolepton\_hadron”, “parameterisation”. The different types of particles and processes are created in the constructor and methods of the `MuCoolPhysicsList` user class, which inherits from `G4VUserPhysicsList`.

### 2.3 Incident Particles

The user constructs incident particles, interaction verteces, or a beam by typing code in the constructor and methods of the `MuCoolPrimaryGeneratorAction` user class, which inherits from
`G4VUserPrimaryGeneratorAction`.

### 2.4 Stepping Actions

The `MuCoolSteppingAction` user action class inherits from `G4UserSteppingAction`. It allows to perform actions at the end of each step during the integration of the equation of motion. Actions may include killing a particle under certain conditions, retrieving information for diagnostics, and others.

### 2.5 Tracking Actions

The `MuCoolTrackingAction` user action class inherits from `G4UserTrackingAction`. For example, particles may be killed here based on their dynamic or kinematic properties.

### 2.6 Event Actions

The `MuCoolEventAction` user action class inherits from `G4UserEventAction`. It includes actions performed at the beginning or the end of an event, that is immediately before or after a particle is processed through the simulated apparatus.

## 3 Description of the Beam Tools Classes

This Section is devoted to explain how to simulate accelerator elements using the Beam Tools. Brief descriptions of each class and constructor are included.

### 3.1 Solenoids

The Beam Tools provide a set of classes to simulate realistic solenoids. These are BTSheet, BTSolenoid, BTSolenoidLogicVol and BTSolenoidPhysVol.

* The `BTSheet` class inherits from `G4MagneticField`. The class objects are field maps produced by an infinitesimally thin solenoidal current sheet. The class data members are all the parameters necessary to generate analytically a magnetic field in $`r`$-$`z`$ space (there is $`\phi `$ symmetry). No geometric volumes or materials are associated with the `BTSheet` objects. `GetFieldValue` is a concrete method of `BTSheet` inherited from `G4Field`, through `G4MagneticField`. It returns the field value at a given point in space and time.
* The `BTSolenoid` class inherits from `G4MagneticField`. The class objects are field maps in the form of a grid in $`r`$-$`z`$ space, which are generated by a set of `BTSheet`. The sheets and the `BTSpline1D` objects, containing the spline fits of $`B_z`$ and $`B_r`$ versus $`z`$ for each $`r`$ in the field grid, are data members of `BTSolenoid`. No geometric volumes or materials are associated with `BTSolenoid`. The field at a point in space and time is accessed through a `GetFieldValue` method, which performs a linear interpolation in $`r`$ of the spline fit objects.
* The `BTSolenoidLogicVol` class defines the material and physical size of the coil system which is represented by the set of current sheets. A `BTSolenoid` must first be constructed from a list of current `BTSheets`. The `BTSolenoid` object is a data member of `BTSolenoidLogicVol`. The `BTSolenoidLogicVol` class constructor creates `G4Tubs` solid volumes and associated logical volumes for the coil system, the shielding, and the empty cylindric regions inside them. Only the logical volumes are constructed here. No physical placement of a magnet object is done.
* The `BTSolenoidPhysVol` class is the placed version of the
`BTSolenoidLogicVol`. It contains the associated `BTSolenoid` object as a data member, as well as the pointers to the physical volumes of its logical constituents.

Figure 1 shows a group of four solenoidal copper coil systems modeled with four infinitesimally thin sheets equally spaced in radius.

### 3.2 Magnetic Field Maps

The Beam Tools also allow to simulate generic field maps using the `BTMagFieldMap` and `BTMagFieldMapPlacement` classes.

* `BTMagFieldMap` class inherits from `G4MagneticField`. The constructor reads the map information from an `ASCII` file containing the value of the field at a set of nodes of a grid. No geometric objects are associated with the field. The field at a point in space and time is accessed through a `GetFieldValue` method, as in the case of the solenoid.
* The `BTMagFieldMapPlacement` class is a placed `BTMagFieldMap` object. Only the field is placed because there is no coil or support system associated with it.

### 3.3 r.f. Systems: Pill Box Cavities and Field Maps

This section explains how to simulate realistic r.f. systems using Pill Box cavities. The Beam Tools package provides the classes: `BTAccelDevice`, `BTPillBox`, `BTrfCavityLogicVol`, `BTrfWindowLogicVol`, and `BTLinacPhysVol`.

* `BTAccelDevice.hh` class is abstract. All accelerator device classes are derived from this class, which inherits from G4ElectroMagneticField.
* The `BTPillBox` class inherits from `BTAccelDevice` and represents single $`\pi /2`$ Pill Box field objects. No solid is associated with `BTPillBox`. The time dependent electric field is computed using a simple Bessel function. It is accessed through a `GetFieldValue` method. The field is given by:

$`E_z=V_pJ_0\left({\displaystyle \frac{2\pi \nu }{c}}r\right)\mathrm{s}in(\varphi _s+2\pi \nu t)`$ (1)
$`B_\phi ={\displaystyle \frac{V_p}{c}}J_1\left({\displaystyle \frac{2\pi \nu }{c}}r\right)\mathrm{c}os(\varphi _s+2\pi \nu t)`$ (2)
where $`V_p`$ is the cavity peak voltage, $`\nu `$ the wave frequency, $`\varphi _s`$ the synchronous phase, and $`J_{0,1}`$ the Bessel functions evaluated at $`\left(\frac{2\pi \nu }{c}r\right)`$.
* The `BTrfMap` class also inherits from BTAccelDevice. The class objects are electromagnetic field maps which represent an r.f. cavity. In this way, complex r.f. fields can be measured or generated and later included in the simulation. The field map, in the form of a grid, is read in the `BTrfMap` constructor from an `ASCII` file. The `BTrfMap` object is a field, with no associated solid. A `GetFieldValue` method retrieves the field value at a point in space and time.
* The `BTrfCavityLogicVol` class constructor creates solid and logical volumes associated with the r.f. field classes. In the case of a map, a vacuum cylinder ring represents its limits. In addition to geometric and material parameters of the cavity, the class contains field and accelerator device information.
* The `BTrfWindowLogicVol` class is used with `BTCavityLogicVol` to create the geometry and logical volume of r.f. cavity windows, including the support structure, which may be placed to close the cavity iris at the end cups.
* The `BTLinacPhysVol` class is a placed linac object. A linac is a set of contiguous r.f. cavities, including the field, the support and conductor material, and windows. The `BTLinacPhysVol` constructor is overloaded. One version places a linac of Pill Box cavities and the other places field maps.

Fig. 2 shows a Pill Box cavity (in red) with windows. It also shows a cooling channel where solenoids are embedded in large low frequency cavities. Since the beam circulates inside the solenoid, the cavity is represented by a field map (in red) restricted to a cylindric volume with radius slightly smaller than the inner radii of the magnets.

### 3.4 Tuning the r.f. Cavity Phases

One of the critical elements of an accelerator simulation is the “r.f. tuning”. Each cavity must be operated at the selected synchronous phase at an instant coincident with the passage of the beam. The r.f. wave must be therefore synchronized with the beam, more specifically, with the region of beam phase space that the user needs to manipulate. For this, there is the concept of a reference particle, defined as the particle with velocity equal to the phase velocity of the r.f. wave. If the kinematic and dynamic variables of the reference particle are set to values which are coincident with the mean values of the corresponding variables for the beam, the r.f. system should affect the mean beam properties in a similar way it affects the reference particle.

The Beam Tools allow the use of a “reference particle” to tune the r.f. system before processing the beam. The time instants the particle goes through the phase center of each cavity are calculated and used to adjust each cavity phase to provide the proper kick, at the selected synchronous phase.

### 3.5 Absorbers

The Beam Tools provide a set of classes to simulate blocks of material in the path of the beam. The constructors create the solid, logical, and physical volumes in a single step. They are all derived from the abstract class of absorber objects `BTAbsObj`.

* `BTCylindricVessel` is a system with a central cylindric rim, and two end cup rims with thin windows of radius equal to the inner radius of the vessel. The material is the same for the vessel walls and windows, and the window thickness is constant. The vessel is filled with an absorber material.
* Two classes are available to simulate absorber lenses:
`BTParabolicLense` and `BTCylindricLense`. The first one is a class of parabolic objects with uniform density, and the second a cylinder object with the density decreasing parabolically as a function of radius. From the point of view of the physics effect on the beam, both objects are almost equivalent. The `BTParabolicLense` is built as a set of short cylinders. The radius is maximum for the central cylinder and reduces symmetrically following a parabolic equation for the others in both sides. The `BTCylindricLense` object is built from concentric cylinder rings of the same length, different radius, and different densities to mimic a real lens.

The gray cylinder in Fig. 3 is a schematic representation of a liquid hydrogen vessel with aluminum walls and windows. Figure 3 also shows a set of six parabolic lenses in the center of a complex magnetic system. The lenses are placed to mitigate the effect of the decrease in $`p_z`$ at large radii in a magnetic field flip region, using an emittance exchange mechanism.

Wedge absorbers are also useful in some cases. They can be easily constructed using the Geant4 trapezoid shape `G4Trap`.

## 4 Applications to Neutrino Factory Feasibility Studies

The neutrino beam in a Neutrino Factory would be the product of the decay of a low emittance muon beam. Muons would be the result of pion decay, and pions would be the product of the interaction of an intense proton beam with a carbon or mercury target. Thus the challenge in the design and construction of a Neutrino Source is the muon cooling section, aimed to reduce the transverse phase space by a factor of ten, to a transverse emittance of approximately $`\epsilon _x`$1 cm.

The ionization cooling technique uses a combination of linacs and light absorbers to reduce the transverse emittance of the beam, while keeping the longitudinal motion under control. There are two competing terms contributing to the change of transverse emittance $`\epsilon _x`$ along the channel. One is a cooling term, associated with the process of energy loss, and the other is a heating term related to multiple scattering.

### 4.1 The Double Flip Cooling Channel

The double flip cooling channel is a system consisting of three homogeneous solenoids with two field-flip sections. The first flip occurs at a relatively small magnetic field, B=3 T, to keep the longitudinal motion under control. The field is then increased adiabatically from -3 to -7 T, and a second field flip performed at B=7 T. Figure 4 shows a side view of a lattice unit cell, consisting of a six 201 MHz Pill Box cavities linac and one liquid hydrogen absorber, inside a solenoid. Details on the design and performance of this channel are available in Ref. .

### 4.2 The Helical Channel

The helical channel cools both in the transverse and longitudinal directions. The lattice is based on a long solenoid with the addition of a rotating transverse dipole field, lithium hydride wedge absorbers, and 201 MHz r.f. cavities. Figure 4.4 shows a side view of the helical channel, including the wedge absorbers, idealistic (thin) r.f. cavities, and the trajectory of the reference particle. The design details and performance of this channel are described in Ref. .

### 4.3 The Low Frequency Channel

This is a design based on 44/88 MHz r.f. technology. A unit cell is composed of four solenoids embedded in four r.f. cavities, followed by a liquid hydrogen absorber. Figure 2 shows a unit cell of the low frequency channel, including the solenoids, the absorber, and the relevant section of the r.f. field map (inside the magnets). More information about this channel may be found in Ref. .

### 4.4 Other Systems

Among other simulations performed with the Beam Tools for Geant4 we may cite: the Alternate Solenoid Channel (sFoFo) , and a High Frequency Buncher/Phase Rotator scheme for the neutrino factory .

## 5 Summary

The Beam Physics Tools for Geant4 are used in numerous accelerator studies, reported in conference proceedings and proposals. Geant4 is especially suited to systems where accelerators, shielding, and detectors must be studied jointly with a simulation. The Beam Tool libraries, a software reference manual, and a user’s guide, are available from the Fermilab Geant4 web page .

We thank Mark Fishler, for contributing the data cards and and spline fit classes, Jeff Kallenbach for helping with visualization issues, and Walter Brown for providing C++ consultancy. We are also grateful to the Geant4 Collaboration for answering our questions. In particular, we thank J. Apostolakis, M. Asai, G. Cosmo, M. Maire, L. Urban, and V. Grichine.


format: text

text:
# A Bounded Rational Driver Model

## I Model description

Assume that at a certain instant of time $`t`$ the driver has decided to correct the car motion and chosen the acceleration $`a(t)`$ (Fig. 2). As discussed above, the optimal path $`\left\{𝔥_{\text{opt}}(𝔱,t)\right\}`$ of the further motion ($`𝔱>t`$) is too complex for her to compute and to follow it. So, she regards the path $`\{𝔥_a(𝔱,t):𝔞(𝔱,t)=a(t)\}`$ characterized by the constant acceleration as the optimal one.

A certain time interval $`\tau _a`$ later, the driver has to correct the car motion again. This can be done by shifting the current acceleration $`a(t+\tau _a)`$ towards the desired optimal value $`a_{\text{opt}}(t+\tau _a)=_𝔱^2𝔥_{\text{opt}}(𝔱,t+\tau _a)|_{𝔱=t+\tau _a}`$ known to her approximately:

$$a(t+\tau _a)a(t)=C\left(a_{\text{opt}}(t+\tau _a)a(t)\right)+a_{\text{rnd}}(t+\tau _a),$$

where $`C1`$ is a constant about unity and the random term $`a_{\text{rnd}}(t+\tau _a)`$ allows for the uncertainty in the driver evaluation of the optimal acceleration at the current time. Its mean amplitude $`a_c`$ characterizes physiological properties of drivers and can be considered constant. Thereby, $`a_{\text{rnd}}(t)a_{\text{rnd}}(t^{})=a_c^2\delta _{t,t^{}}`$, where $`\delta _{t,t^{}}`$ is Kronecker’s delta.

This discrete representation of the car motion correction is converted to a continuous description based on stochastic differential equations. Namely, the above discrete governing equation is reduced to

$$\frac{da}{dt}=\frac{1}{\tau _a}\left(aa_{\text{opt}}(h,v,V)\right)+\eta \xi (t).$$
(1)
Here, $`a_{\text{opt}}(h,v,V)`$ is the optimal acceleration specified by the current values of headway, car velocity, and leading car velocity. The term $`\xi (t)`$ is white noise of unit amplitude which models the uncertainty in the driver evaluation of the optimal motion.

The acceleration increment $`\delta a`$ caused by the random force $`\eta \xi (t)`$ acting during the time $`\tau _a`$ is actually the random component $`a_{\text{rnd}}(t)`$ entering the discrete governing equation. Thus, it follows from the estimate $`(\delta a)^2\eta ^2\tau _a`$ that

$$\eta =\frac{a_c}{\sqrt{\tau _a}}.$$
(2)
The time scale $`\tau _a`$ of the driver control over the car motion depends on the state $`(h,v,V,a)`$. Thus, the stochastic differential equation (1) contains multiplicative noise. So its type with respect to the corresponding Fokker-Planck equation has to be specified. The adopted driving strategy (Fig. 2) implies that all the characteristics of correcting the car motion are determined by its state at the “terminal” point $`t+\tau _a`$ rather than at the “initial” point $`t`$. Therefore, it is reasonable for Eq. (1) to be of Klimontovich type or, according to the classification in MH , to describe a “postpoint” random process.

To complete the model, $`a_{\text{opt}}(h,v,V)`$ and $`\tau _a(h,v,V,a)`$ have to be specified. The simple ansatz

$$a_{\text{opt}}(h,v,V)=\frac{1}{\tau }\left[\left(vV\right)\frac{1}{\tau }g_h\left(hh_V\right)\right]$$
(3)
is well justified, at least, near the stationary state of the car motion, $`v=V`$ and $`h=h_V`$. It should be noted that similar ideas about $`a_{\text{opt}}(h,v,V)`$ and a dependence of $`\tau _a`$ on the motion state had been discussed already in Ref. Hel . (See also D3  for a discussion.)

Here, $`\tau `$ is the characteristic time of the velocity variations and the constant $`g_h1`$. The limit $`g_h1`$ deserves special attention because it is just the condition that a driver, at first, prefers to eliminate the velocity difference $`vV`$ between her car and the car ahead and only then optimizes the headway. In this case the optimal dynamics of car motion, i.e., the car dynamics governed by the relation $`a=a_{\text{opt}}(h,v,V)`$ is a pure fading relaxation towards the stationary state. Conversely, the model under consideration predicts complex oscillations in the car motion. Note, that the adopted assumption about the value of the coefficient $`g_h`$ can be justified by applying to the general principles of the car motion we2 .

If the car motion state is far from equilibrium the necessity for correcting the velocity and headway distance is obvious. In this case it is natural to suppose that the characteristic time interval $`\tau _a`$ between sequential attempts to correct the car motion should be comparable to $`\tau `$ which characterizes the velocity variations, i.e., $`\tau _a\tau /g_v`$. Here, $`g_v1`$ is an additional model parameter. When the car motion comes close to the equilibrium and the inequality $`|a_{\text{opt}}(h,v,V)|a_c`$ is fulfilled the uncertainty $`a_{\text{rnd}}(t+\tau _a)`$ in evaluating the optimal acceleration becomes significant. Under such conditions there is no reason for the driver to affect the car motion and she may not correct it at all. It means that the car motion control is depressed and, correspondingly, the correction time interval $`\tau _a`$ grows dramatically inside a domain $`_u`$ of the phase plane $`\{h,v,V\}`$ where the inequality $`|a_{\text{opt}}(h,v,V)|a_c`$ holds.

To compute the function $`\tau _a(h,v,V,a)`$, the boundary of the domain $`_u`$ has to be analyzed. Note, that the acceleration itself enters the driver’s perception of motion quality: without any reason, a driver prefers not to accelerate at all. When the car motion control is active the estimate $`\dot{a}a/\tau _a`$ by virtue of Eq. (1) can be adopted. So, the boundary of the domain $`_u`$ is specified by $`a_{\text{opt}}^2(h,v,V)+\mu ^2a^2a_c^2`$, where $`\mu 1`$ is a certain coefficient about unity. Assuming the variables $`h`$, $`v`$, $`a`$ to be independent of one another inside $`_u`$ and averaging the latter expression over $`_u`$ its boundary $`\mathrm{\Phi }(h,v,V,a)1`$ can be derived:

$$\mathrm{\Phi }(h,v,V,a)=\frac{\left(vV\right)^2}{a_c^2\tau ^2}+g_h^2\frac{\left(hh_V\right)^2}{a_c^2\tau ^4}+\mu ^2\frac{a^2}{a_c^2}.$$
(4)
If $`\mathrm{\Phi }(h,v,V,a)1`$ the driver activity in correcting the car motion is depressed completely. Otherwise, $`\mathrm{\Phi }(h,v,a)1`$, the driver controls the car motion actively. This is described by the dependence of the correction time interval $`\tau _a`$ on the car motion state,

$$\frac{1}{\tau _a}=g_v\mathrm{\Omega }\left(\mathrm{\Phi }(h,v,V,a)\right)\frac{1}{\tau }.$$
(5)
The form of the function $`\mathrm{\Omega }\left(x\right)`$ is illustrated in Fig. 3. Equation (1) together with expressions (2), (4), and (5) form the proposed car following model with bounded rational drivers.

When $`g_h>g_v\mathrm{\Omega }(0)`$ the stationary motion with $`v=V`$ and $`h=h_V`$ is unstable, leading to non-damped but bounded oscillations in the headway and velocity of the following car. The particular form of the function $`\mathrm{\Omega }\left(x\right)`$ is of minor importance, it is only necessary that its value inside $`_u`$ to be small in comparison with the ratio $`g_h/g_v`$. When analyzing the model numerically the following ansatz

$$\mathrm{\Omega }(x)=\mathrm{exp}[(x1)/\mathrm{\Delta }]/(\mathrm{exp}[(x1)/\mathrm{\Delta }]+1)$$

is used, with the parameter $`\mathrm{\Delta }0.2`$. Below, numerical results will be presented that demonstrate the characteristic properties of the developed model.

Figure 4 displays an example of this dynamic in the $`hv`$-phase plane for the dimensionless headway $`x=(hh_V)/(a_c\tau ^2)`$ and the relative car velocity $`u=(vV)/(a_c\tau )`$. As seen in Fig. 4, the behavior of this model is qualitatively similar to the empirical data in Fig. 1.

Preliminary results have shown that, first, the quasi-period of these oscillations in the car velocity is equal to $`\tau `$ times a numerical factor (about ten) depending weakly on the model parameters. For $`\tau 1`$ s this period is similar to the observed quasi-period. Second, the amplitude of velocity oscillations does not change substantially as the model parameters vary and is about $`a_c\tau `$. By visually comparing Figs. 1, 4 the estimate $`a_c0.3`$ m/s2 is obtained. It should be noted that the amplitude of the acceleration oscillations exceeds $`a_c`$ by a numerical factor about three. Third, the amplitude of the headway super-oscillations, in contrast, depends essentially on the parameter $`g_v`$, enabling one to fix this parameter based on experimental data.

## II Summary

A model regarding the bounded rational behavior of car drivers has been supposed in this contribution. It takes into account that drivers, although having detailed ideas about their preferred driving strategy, are not able to control this driving strategy sufficiently precisely. Namely, drivers introduce three main sources of error into the optimal driving strategy: instead of keeping track of the changes in acceleration they simply choose a constant one, that additionally is not the optimal one but blurred by noise. This noise models the inability of drivers to evaluate exactly the very complex integrations leading to an optimal driving strategy. Therefore, the need to correct the motion from time to time arises, with the correction time intervals distributed randomly but inversely proportional to the deviation from the desired optimal acceleration.

It is shown, that these ideas can be captured in a simple model for the car-following dynamics, however at the cost of introducing a non-Newtonian term, the jerk (change in acceleration). The benefit of doing so is that the resulting model has smooth trajectories in headway, velocity and acceleration but still being a stochastic one. This discerns the approach proposed here from almost all models of car-following introduced so far.

Although the trajectories generated by this model have some similarities with real car-following data, the approach proposed here still needs thorough testing with empirical data. This will be done in the near future and will be reported soon.

These investigations were supported in part by RFBR Grants 01-01-00389, 00439 and INTAS Grant 00-0847.


text:
# Listening to the coefficient of restitution and the gravitational acceleration of a bouncing ball

## Abstract

We show that a well known method for measuring the coefficient of restitution of a bouncing ball can also be used to obtain the gravitational acceleration.

Three contributions to this journal have described how to measure the coefficient of restitution between a ball and a flat surface using the sound made by the collision of the ball with the surface Bernstein ; Smith ; Stensgaard . The procedure reported in these papers is to drop the ball vertically on a horizontal surface, allow it to bounce several times, while recording the sound produced by the impacts. Analysis of the recording gives the time intervals between successive rebounds, and from these the coefficient of restitution is obtained.

The evolution of the techniques described in these papers is a nice example of how the development of microcomputers has changed the science teaching laboratory. In 1977, Bernstein Bernstein  detected the sound with a microphone, amplified and filtered the signal, and fed it to a pen recorder. Smith, Spencer and Jones Smith , in 1981, connected the microphone to a microcomputer via a homemade data collection and interface circuit, and then uploaded the resulting data to a larger computer for analysis and graphical display. In 2001, Stensgaard and Lægsgaard Stensgaard  used the microphone input of a PC sound card to make the recording, reducing the experimental equipment to basically a standard microcomputer.

To see how the coefficient of restitution $`ϵ`$ is related to the time between bounces, note that if $`ϵ`$ is constant (independent of velocity), and air resistance is negligible, the velocity of the ball just after the $`n`$th bounce on the fixed surface is given by

$$v_n=v_0ϵ^n$$
(1)
where $`v_0`$ is the velocity just before the first impact. The time-of-flight $`T_n`$ between the $`n`$th and $`(n+1)`$th collisions is proportional to $`v_n`$,

$$T_n=\frac{2v_n}{g},n=1,2,\mathrm{}$$
(2)
where $`g`$ is the gravitational acceleration. Thus

$$T_n=T_0ϵ^n,$$
(3)
where we have defined $`T_02v_0/g`$ . Taking the logarithm of both sides of Eq. (3) we obtain

$$\mathrm{log}T_n=n\mathrm{log}ϵ+\mathrm{log}T_0,$$
(4)
so that the plot of $`\mathrm{log}T_n`$ vs $`n`$ is a straight line of slope $`\mathrm{log}ϵ`$ and intercept $`\mathrm{log}T_0`$. Thus, as long as it is independent of velocity, the coefficient of restitution can be obtained by fitting the straight-line of Eq. (4) to the time-of-flight data.

The purpose of this note is to point out that this straight-line fit can also be used to determine another physical quantity of interest: the gravitational acceleration $`g`$. The (rather simple) observation is that, if the ball is released from a known height $`h`$, then $`T_0=(8h/g)^{1/2}`$, and

$$g=\frac{8h}{T_{0}^{}{}_{}{}^{2}}.$$
(5)
Thus, just as the slope parameter of Eq. (4) fixes the coefficient of restitution, the intercept parameter determines the acceleration of gravity (if the easily measured initial height $`h`$ is known).

In order to check how this works in practice, we have dropped a “superball” from a measured height onto a smooth stone surface and recorded the sound produced by the successive impacts. The recording was made with the microphone and sound card of a PC running Windows, using the sound recorder program that comes with the operating system. The sampling frequency was 22,050 Hz, resulting in a time resolution of 45 $`\mu `$s. The audio file, stored in the binary WAV format, was converted to ASCII text format with the shareware program awave audio awave . The recorded signal is plotted in Fig. 1, where the pulses corresponding to individual impacts are easily recognized (only the first six collisions are shown). We have used 8-bit resolution in the recording, so that data values can only go from 0 to 255. The no-signal value corresponds to 128.

The time intervals $`T_n`$ between collisions $`n`$ and $`n+1`$ were obtained directly through inspection of the ASCII sound file. They are plotted in Fig. 2 (in logarithmic scale) as a function of $`n`$. The least-squares fit of the $`T_n`$ data set to Eq. (4) gives

$`ϵ=0.9544\pm 0.0002,`$ (6a)
$`T_0=0.804\pm 0.001\text{ s}.`$ (6b)
The best-fit line is also shown in Fig. 2.

The ball was released from a height $`h=79.4\pm 0.1`$ cm above the surface. Taking this and the adjusted $`T_0`$ into Eq. (5), we obtain for the gravitational acceleration

$$g=982\pm 3\text{ cm/s}^2.$$

For comparison, the value of $`g`$ in Rio de Janeiro is 978.8 cm/s2.

The applicability of the method described above depends on $`ϵ`$ being constant over the range of impact velocities involved in the experiment. That this condition is satisfied in the present case is seen in Fig. 3, where the coefficient of restitution for an impact at velocity $`v_n`$, $`ϵ=v_{n+1}/v_n=T_{n+1}/T_n`$, is plotted as a function of $`T_n`$ (recall that $`v_nT_n`$, see Eq. (2)). The coefficients of restitution for the different impacts are all very close to the least-squares value given in Eq. (6), indicated by the dashed line in Fig. 3.

A case in which the coefficient of restitution depends on the velocity is shown in Fig. 4, where we display the times of flight of a superball dropped from $`h=27.5\pm 0.1`$ cm onto a wood surface. A plot of $`ϵ`$ at each collision, shown in Fig. 5, reveals a clear dependence of the coefficient of restitution on the time-of-flight (or impact velocity). Assuming a linear relation between $`ϵ`$ and $`T`$, as suggested by Fig. 5,

$$ϵ=ϵ_0(1+\alpha T),$$
(7)
we obtain an extension of Eq. (3)

$$T_n=T_0ϵ_{0}^{}{}_{}{}^{n}\underset{i=0}{\overset{n1}{}}(1+\alpha T_i).$$
(8)
The least-squares fit of Eq. (8) to the data shown in Fig. 4 gives

$`ϵ_0=0.921\pm 0.001,`$ (9a)
$`\alpha =0.078\pm 0.003\text{ s}^1,`$ (9b)
$`T_0=0.4752\pm 0.0005\text{ s}.`$ (9c)
The curves corresponding to these parameters are also shown in Figs. 4 and 5. The above value for $`T_0`$ yields

$$g=974\pm 5\text{ cm/s}^2,$$

again a very reasonable value. Consideration of the velocity dependence of the coefficient of restitution was important in order to get an accurate result; had we assumed a constant $`ϵ`$, we would obtain $`g=935\pm 10`$ cm/s2.

To summarize, we have seen that the value of the gravitational acceleration is a useful by-product of experiments devised to “hear” the coefficient of restitution of a bouncing ball. The measurement of $`g`$ is particularly simple if the coefficient of restitution is independent of the impact velocity, but more complicated cases can also be handled.

After this work was completed we learned of a recent paper by Cavalcante *et al.* Cavalcante , in which $`g`$ was measured using the sound of a bouncing ball. The analysis presented in the paper is, however, somewhat different from ours. Another related reference is the article by Guercio and Zanetti Guercio  in this journal.


text:
# Is KT Mon a Classical Nova?

## 1 Introduction

KT Mon (Nova Mon 1942) was discovered by A. N. Vyssotsky on a Harvard plate taken on 1943 January 2. The object was first recorded at mpg = 10.3 on 1942 December 30. Gaposchkin (1954) presented a light curve and a description of a spectrum taken by Vyssotsky on the discovery plate. According to Gaposchkin (1954), the object rather rapidly faded (2 mag in 27 d). The object was fainter than mpg = 13.35, 13.9 d prior to the first positive record on 1942 December 30. Gaposchkin (1954) described the low-resolution spectrum by Vyssotsky as “broad bright lines H$`\beta `$, H$`\gamma `$, H$`\delta `$, 4686 HeII and 4640 NIII”. From these findings, Gaposchkin (1954) concluded that KT Mon is a nova. Duerbeck (1987) classified the object as an established moderately fast nova.

However, several authors have noted the peculiarity of KT Mon. Payne-Gaposchkin (1977) deduced the distances of past novae, and found three stars (KT Mon, CG CMa and WX Cet) are more than 15 kpc distant from the Galactic Center. The two objects in this list, CG CMa (Kato et al., 1999; Duerbeck et al., 1999) and WX Cet (Bailey, 1979; O’Donoghue et al., 1991; Kato et al., 2001) are now established to be SU UMa-type dwarf novae (For a recent review of dwarf novae and SU UMa-type dwarf novae, see Osaki (1996) and Warner (1995), respectively.) Regarding KT Mon, Payne-Gaposchkin (1977) described that “KT Mon: whose spectrum observed by Vyssotsky, and whose light curve obtained by Gaposchkin (1954) place it unequivocally as a nova” and suggested that “further study might find KT Mon to be recurrent” based on this unusual distance estimate. Similar long-distance estimates have been given in Warner (1987); Shafter (1997). No quiescent counterpart has been reported down to mag 20 (Duerbeck, 1987).

## 2 Overlooked Unusual Features of KT Mon

In addition to the unusually long distance estimate (Payne-Gaposchkin, 1977), which originally questioned the classical nova-nature of KT Mon, we further noted a previously overlooked peculiarity from the available material.

The early presence of high excitation lines, particularly HeII $`\lambda `$4686 and NIII $`\lambda `$4640, is unusual for typical development of nova spectra (Payne-Gaposchkin, 1957; Bode & Evans, 1989; Williams, 1992), since strong appearance of NIII $`\lambda `$4640 Bowen brend requires high-energy photons (van Paradijs & McClintock, 1995). Such a condition is hard to achieve during the early stage of a classical nova outburst when thermal emissions from a lower ($`10^4`$ K) temperature expanded atmosphere of the white dwarf dominates. High-excitation emission lines can appear when the contracting nova (pseudo)photosphere reaches $`T_{\mathrm{eff}}2.5\times 10^4`$ K (see e.g. Bode & Evans (1989) Chap. 5), which then produces sufficient high-energy UV photons. This condition is expected to be achieved at $``$2.5 mag below the optical maximum, which agrees with the general observational records (Williams et al., 1994). From an examination of photographic records of past nova spectra, these lines were reported to appear at around 4 mag below the optical maximum. Only a few very fast novae (e.g. V2487 Oph, Filippenko et al. (1998)) and fast recurrent novae (e.g. U Sco, Sect. 3.1) are known to exhibit these high-excitation lines soon after their maxima, which are presumably caused by an unusually thin ejecta. Since the decline rate of KT Mon was not particularly fast, this object apparently did not achieve this condition at the time of Vyssotsky’s observation.

A relatively large ($``$1.5) color index (blue-red) inferred by Gaposchkin (1954) also needs to be reexamined. Taking into the recent observation of H$`\alpha `$ flux of post-outburst novae (Ciardullo et al., 1990) into account, the contribution of the H$`\alpha `$ line to the $`R`$ band is estimated to $``$1.0 mag. Although much of the originally reported color index was tentativity attributed to interstellar reddening in interpreting KT Mon as being a classical nova (Payne-Gaposchkin, 1977; Shafter, 1997), a more recent estimate of the maximum reddening in this direction $`E(BV)`$ = 0.545 or $`A_V`$ = 1.806 (Schlegel et al., 1998) suggests that the observed color more strongly reflects the contribution of the H$`\alpha `$ emission line. If the ejecta was indeed unusually thin as in recurrent novae, the contribution of the H$`\alpha `$ line, however, might be insufficient (cf. Sekiguchi et al. (1988)) to explain the overall color index.

## 3 New Interpretations

In recent years, more categories of eruptive objects have been identified to show high-excitation emission lines during their early stage of outbursts: U Sco-type recurrent novae (Rosino & Iijima, 1988; Sekiguchi et al., 1988; Iijima, 2002), WZ Sge-type dwarf novae (Baba et al., 2002; Kuulkers et al., 2002) and X-ray transients (Tanaka & Shibazaki, 1996; van Paradijs & McClintock, 1995). Considering the strong expected impact if KT Mon indeed belongs to any of these objects, we further examine these possibilities.

First of all, we have performed the prenova search with available plate scans (the DSS1 red, the DSS2 red, blue, infrared) and available catalogs in the VizieR service. Gaposchkin (1954) gave a position of nova as 6h 19m 58s.8, +5 29’ 46” (equinox 1900), which precesses to 6h 25m 18s.9, +5 26’ 28” (J2000.0). Duerbeck (1987) re-examined the Harvard plates and give the position as 6h 25m 18s.46, +5 26’ 31”.7. Duerbeck (1987) further indicated that the identification by Khatisov (1971) (who used the crude GCVS position for identification) with a GSC-cataloged star having position end figures 19s.69, 32”.88 (GSC 2.2.1) is incorrect. Taking this indication into consideration, the position of Duerbeck (1987) must have an accuracy better than 10”. We confirmed that there is no promising optical, IR, and X-ray counterpart within 10” from the Duerbeck (1987) position. We have thus confirmed the safe upper limit of the quiescent KT Mon to be $`V`$ = 20.

### 3.1 Recurrent Nova?

This possibility was originally considered by Payne-Gaposchkin (1977), but has been overlooked to date. Since the outburst amplitude of KT Mon is larger than 10, we can safely exclude the possibility of a recurrent nova with a giant secondary (Anupama & Mikolajewska, 1999). The remaining possibilities are a U Sco-like object or an IM Nor/CI Aql-like object (Kiss et al., 2001; Matsumoto et al., 2001; Kato et al., 2002b). The former class represents extremely rapidly evolving novae and the latter slower ones. The reported light curve of KT Mon (Gaposchkin, 1954) (decline by 2 mag in 27 d) more suggests the latter class. However, the spectroscopic feature of KT Mon (presence of the HeII emission line) does not resemble the early spectra of IM Nor and CI Aql (see Kiss et al. (2001); Matsumoto et al. (2001)), which did not show high excitation lines. Since the difference between these classes can be primarily attributed to the mass of the white dwarf (Hachisu et al., 2000a; Hachisu & Kato, 2001), there still remains the possibility that KT Mon is essentially a U Sco-type object, with a slightly less massive white dwarf.

Based on modern calculations (Hachisu et al., 2000a, b), the optical maximum of U Sco corresponds to $`M_\mathrm{V}`$ = $``$6.4$`\pm `$0.4. By assuming $`0A_V1.8`$, and the same peak maximum $`M_\mathrm{V}`$ as in U Sco, the range of acceptable distance of KT Mon becomes 6$`d`$26 kpc. This range can only slightly reduce the Galactocentric distance questioned in Payne-Gaposchkin (1977). Even if the distance is acceptable, the upper limit quiescent $`M_\mathrm{V}`$ = +3.3$`\pm `$0.4 is hard to accept for a recurrent nova which requires a high accretion rate (Hachisu et al., 2000b). If KT Mon is indeed a recurrent nova, there is a need for a special unidentified mechanism to reduce the quiescent luminosity.

### 3.2 WZ Sge-Type Star?

It has been now widely demonstrated that some of WZ Sge-type dwarf novae show strong HeII and CIII/NIII emission lines during their outbursts. There properties are not inconsistent with the existing observation of KT Mon. Since KT Mon is 3.3 mag fainter than WZ Sge at maximum, a reasonable upper limit of $`d`$ 200 pc can be derived from the recent distance estimate of WZ Sge (45 pc, J. Thorstensen (2001), cited in Steeghs et al. (2001)). If this interpretation is correct, the quiescent upper limit of KT Mon corresponds to $`M_\mathrm{V}=+13.5`$. This value is not inconsistent with a recent example of a very faint WZ Sge-type object van Teeseling et al. (1999). Although an argument exists regarding the distance of V592 Her itself (Kato et al., 2002a), $`M_\mathrm{V}=+13.5`$ is not incompatible with a combination of a binary with a cool white dwarf and a brown dwarf van Teeseling et al. (1999). If KT Mon belongs to WZ Sge-type dwarf novae, this object is a strong candidate for a close binary containing a brown dwarf.

This interpretation may seem to show a difficulty in interpreting the reported large color index (Gaposchkin, 1954) in its late-stage light curve. In recent years, however, WZ Sge-type dwarf novae can become exceptionally red ($`BR+1.0`$ has been reported, see e.g. Patterson et al. (1998)) during the late stage of their outbursts. This red color may explain some part of the color index reported by Gaposchkin (1954). WZ Sge-type dwarf novae are found to frequently show “rebrightenings”. In WZ Sge itself, such rebrightenings with amplitudes of 1–2 mag (Ishioka et al., 2002). Since these rebrightenings show a very rapid rise and fall, a scatter in the fragmentary light curve in Gaposchkin (1954) may reflect such a phenomenon.

We also note that V4338 Sgr (Nova Sgr 1990), with similar spectroscopic characteristics to those of KT Mon, has been classified as a possible WZ Sge-type dwarf nova (Wagner et al., 1990).

### 3.3 X-ray Transient?

The presence of strong Balmer and HeII and CIII/NIII emission lines is also very characteristic to X-ray transients (mostly black-hole candidates; see comprehensive reviews van Paradijs & McClintock (1995); Tanaka & Lewin (1995); White et al. (1995)). It will be noteworthy that the X-ray transient V404 Cyg = GS 2023+338 (Wagner et al., 1991; Casares et al., 1991) had been considered to be a classical nova based on its 1938 outburst observation (Duerbeck, 1987). KT Mon therefore can be a X-ray transient. In particular, the description of the spectrum closely agrees to that of V404 Cyg Wagner et al. (1991). The decline rate (0.07 mag d-1, see Gaposchkin (1954)) is not inconsistent with the statistics of optical properties of X-ray transients (Chen et al., 1997). Since the optical maximum of KT Mon is 0.9 mag brighter than that of the closest “classical” X-ray transient V616 Mon ($`d`$ = 1.4 kpc, Esin et al. (2000)), this interpretation would bring KT Mon to $`d`$1 kpc. This possibility makes KT Mon a candidate for a closest black-hole binary. Although a different distance can be acceptable considering wide diversity of properties of X-ray transients (Chen et al., 1997; Mirabel et al., 2001), the strong presence of HeII and CIII/NIII emission lines is more consistent with a high-luminosity X-ray transient. The distance $`d`$1 kpc corresponds to a quiescent upper limit of $`M_\mathrm{V}=+10`$. This value is not inconsistent with a short-period system accreting at a very low quiescent accretion rate (cf. van Paradijs & McClintock (1994, 1995)). The reported color index is mildly consistent with the reported reddening of V616 Mon $`E(BV)=0.350.9`$, Marsh et al. (1994)) combined with the redder color during the late stage of outbursts of X-ray transients (e.g. King et al. (1996)).

## 4 Summary

KT Mon has been considered as a classical nova in 1942 based on its light curve and spectroscopy. However, we noticed a previously overlooked unusual feature: the early presence of high excitation emission lines such as HeII. We examine three new possibilities (recurrent nova, WZ Sge-type dwarf nova, X-ray transient). We have found that the possibilities of a WZ Sge-type dwarf nova or an X-ray transient are more consistent with the modern knowledge, although the possibility of a recurrent nova (or a nova with an unusual spectroscopic evolution) is not completely excluded. If KT Mon is a WZ Sge-type dwarf nova, the system should have a brown dwarf. If KT Mon is an X-ray transient, the system is a good candidate for the nearest black hole. In all possibilities, we can expect a recurrent outburst in meaningfully near future (10–100 yr depending on the classification). Detailed observations during the possible next outburst, as well as deep quiescent observations, for this obscure object are strongly recommended.

This work is partly supported by a grant-in aid \[13640239 (TK), 14740131 (HY)\] from the Japanese Ministry of Education, Culture, Sports, Science and Technology. This research has made use of the Digitized Sky Survey producted by STScI, the ESO Skycat tool, the VizieR catalogue access tool.


text:
# Tunneling time scale of under-the-barrier forerunners

## Abstract

Time-dependent analytical solutions to Schrödinger’s equation with quantum shutter initial conditions are used to investigate the issue of the tunneling time of forerunners in rectangular potential barriers. By using a time-frequency analysis, we find the existence of a regime characterized by the opacity of the barrier, where the maximum peak of a forerunner measured at the barrier transmission edge $`x=L`$ corresponds to a genuine tunneling process. The corresponding time scale represents the tunneling time of the forerunner through the classically forbidden region.

Tunneling, that describes the possibility that a particle traverses through a classically forbidden region, is one of the paradigms of quantum mechanics. In the energy domain one solves Schrödinger’s equation at a fixed energy $`E`$ to obtain the probability of transmission through a barrier region, a subject discussed in every quantum mechanics textbook. In the time domain, however, there are still aspects open to scrutiny. A problem that has remained controversial over the years is the tunneling time problem, that may be stated by the question : How long does it take to a particle to traverse a classically forbidden region? Different authors have proposed and defended different views in answering the above questiontraversal . Some authors have argued that it might not be a unique definition of the tunneling time since different arrangements for the tunneling process may lead to different relevant time scales. An interesting idea that arises from the above considerations is that the relevance of the different tunneling time definitions may depend on the context where they may become physically significant quantitiesgcvdm02 .

In recent work we have investigated the effect of the transient solutions to the time-dependent Schrödinger’s equation for cutoff wave initial conditions (quantum shutter) on the tunneling processgcr97 ; gcv01 ; gcv02 . In particular for tunneling through a barrier we found that just across the tunneling barrier, the probability density as a function of time may exhibit a sharp peak maximum that we called time domain resonance. The peak value $`t_{max}`$ of that forerunner represents the largest probability of finding the particle at the barrier width $`L`$. We found also that the behavior of $`t_{max}`$ as a function of the barrier width $`L`$ exhibits for small values of $`L`$ a basin region, followed by a region where $`t_{max}`$ grows linearly as $`L`$ increasesgcv01 . More recently, in collaboration with Delgado and Mugagcvdm02 , we have investigated the time scale for forerunners preceeding the main tunneling signal of a wave created by a source with a sharp onsetmb00 ; vrs02  and by a quantum shuttergcv01  for systems with very large or infinite barrier width. For the particular case of opaque finite barriers gcvdm02  we found a basin regime occurring at fixed positions $`0<x\kappa _0^1`$ along the internal region, where forerunners are dominated by under-the-barrier frequency components. That is, in opaque barrier systems the tunneling forerunners are observed at distances of the order or smaller than the penetration length $`\kappa _0^1`$, where $`\kappa _0=[2m(VE)]^{1/2}/\mathrm{}`$, with $`E`$ and $`V`$ corresponding to the incidence energy and the barrier height, respectively. However, at distances $`\kappa _0^1<xL`$ the non-tunneling components eventually dominate the time evolution process, and hence it is not possible to speak of a genuine tunneling time scale of forerunners at $`x=L`$. We believe that in order to gain more insight on the properties of transient tunneling structures and their corresponding time scales, the existence of tunneling forerunners at the barrier edge $`x=L`$ needs to be investigated.

The aim of this work is to show the existence of a regime along the basin region, where the maximum $`t_{max}`$ of the time domain resonance measured at the barrier width $`L`$ corresponds to under-the-barrier tunneling and consequently provides a genuine time scale for tunneling.

Our approach to the tunneling time problem is based on a model that deals with an explicit solutiongcr97  of time-dependent Schrödinger equation for an arbitrary potential $`V(x)`$ ($`0xL`$) that vanishes outside the internal region. We consider the problem of the time evolution of a cutoff plane wave $`\mathrm{\Psi }(x,k;t=0)=\mathrm{\Theta }(x)(e^{ikx}e^{ikx})`$ following the instantaneous opening at $`t=0`$ of a quantum shutter at $`x=0`$. Along the tunneling region the solution reads,

$`\mathrm{\Psi }^i(x,k,t)=`$ $`\mathrm{\Phi }_k(x)M(y_k)\mathrm{\Phi }_k(x)M(y_k)`$ (1)
$`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}\mathrm{\Phi }_n(x)M(y_{k_n}),(0xL)`$
where the $`\mathrm{\Phi }_{\pm k}(x)`$’s refer to the stationary solutions of the problem and $`\mathrm{\Phi }_n(x)=2iku_n(0)u_n(x)/(k^2k_n^2)`$. Similarly, the solution $`\mathrm{\Psi }^e(x,k;t)`$ for the external or transmitted region ($`xL`$), is given by gcv01 ,

$`\mathrm{\Psi }^e(x,k;t)=`$ $`T_kM(y_k)T_kM(y_k)`$ (2)
$`i{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}T_nM(y_{k_n}),(xL)`$
where the $`T_{\pm k}`$’s refer to the transmission amplitudes, and the factor $`T_n=2iku_n(0)u_n(L)\mathrm{exp}(ik_nL)/(k^2k_n^2)`$. In Eqs. (1) and (2) the coefficients $`\mathrm{\Phi }_n`$ and $`T_n`$ are given in terms of the resonant eigenfunctions $`\{u_n(x)\}`$ with complex energy eigenvalues $`E_n=\mathrm{}^2k_n^2/2m`$, with $`k_n=a_nib_n`$ ($`a_n,b_n>0`$). The resonant sums in Eqs. (1) and(2) run over the full set of complex poles $`\{k_n\}`$. The $`M^{}s`$ are defined as gcr97 ,

$$M(y_q)=\frac{1}{2}e^{imx^2/2\mathrm{}t}w(iy_q)$$
(3)
where $`w`$ is the complex error function abrwtz defined as $`w(z)=\mathrm{exp}(z^2)erfc(iz)`$, with arguments $`y_q(x,t)=e^{i\pi /4}(m/2\mathrm{}t)^{1/2}[x\mathrm{}qt/m]`$, where $`q=\pm k,`$ $`k_{\pm n}`$.

We shall explore the issue of the tunneling time scale in typical one dimensional potential barriers mendez  of height $`V`$ and thickness $`L`$, defined along the interval $`0xL`$. Let us analyze the behavior of the maximum of the time domain resonance, $`t_{max}`$, as a function of the potential barrier width, $`L`$. Let us recall that this transient structure corresponds to the first maximum of the probability density $`|\mathrm{\Psi }^e|^2`$, measured at the fixed position, $`x=L`$. This can be appreciated in Fig. 1 where we plot the $`L`$dependence of $`t_{max}`$ (full dot), corresponding to potential barrier systems with parameters: $`V=0.3`$ eV, incidence energy $`E=\mathrm{}^2k^2/2m=0.001`$ eV, and effective mass for the electron $`m=0.067m_e`$. Here we can clearly observe a basin corresponding to a range of values of the barrier width. We can also appreciate that if $`L`$ is further increased, $`t_p`$ grows linearly with $`L`$. Such a linear regime occurs at large barrier widths. We have argued gcv01  that in this case, the tunneling process is inhibited and that the particle goes mainly over the barrier.

We have recently suggested gcv01  that for small values of the barrier width $`L`$, the basin exhibited by $`t_{max}`$ , is a result of a subtle interplay between tunneling and top-barrier resonant processes. In what follows we shall investigate under what conditions the time scales associated to the basin, are in fact related to a tunneling process. We begin our analysis by choosing the case $`L=4.0`$ nm depicted in Fig. 1, which is located around the minimum of the basin. In Fig. 2 we plot the normalized probability density $`|\mathrm{\Psi }^e|^2`$ (solid line), normalized to the transmission coefficient $`|T_k|^2`$, as a function of time $`t`$ for a fixed value of position $`x=L`$. We can clearly appreciate a time domain resonance peaked at the value $`t_{max}=5.17`$ fs.

In order to determine if this structure is related to a genuine tunneling process we analyze the frequency content of $`|\mathrm{\Psi }^e|^2`$. We compute the local average frequency $`\omega _{av}`$ mb00 ; cohen ,

$$\omega _{av}=\mathrm{Im}\left[\frac{1}{\mathrm{\Psi }^s}\frac{d}{dt}\mathrm{\Psi }^s\right],$$
(4)
and the instantaneous bandwidth $`\sigma `$cohen ,

$$\sigma =\left|\mathrm{Re}\left[\frac{1}{\mathrm{\Psi }^s}\frac{d}{dt}\mathrm{\Psi }^s\right]\right|,$$
(5)
where $`s=i,e`$ according to consider, respectively, the internal or external solutions.

In Fig. 3 we plot the relative average local frequency (relative frequency for short) $`\omega _{av}/\omega _V`$, where $`\omega _V=V/\mathrm{}`$ is the cut-off frequency, along the relevant time interval, discussed in Fig. 2. We can appreciate that in the vicinity of the maximum of the time-domain resonance, $`t_{max}`$, the probability density is composed entirely by under-the-barrier frequency components i.e $`\omega _{av}/\omega _V<1`$. This also occurs at the exact value $`t_{max}`$, also indicated in the figure by a solid square. In the inset of Fig. 3 we plot the instantaneous bandwidth $`\sigma `$ of the spectrogram. Notice the absence of a frequency dispersion around the maximum $`t_{max}`$, i.e., $`\sigma (t_{max})=0`$. The above result indicates that in this case the peak of the time domain resonance, and the values close to it, refer to a genuine tunneling event.

In Fig. 4 we plot the relative frequency $`\omega _{av}/\omega _V`$ associated to different values of the maximum $`t_{max}`$ of the time-domain resonance, measured at different positions along both the internal and external regions of a potential barrier with parameters: $`V=0.3`$ eV, $`L=4.13`$ nm, and $`\alpha =3.0`$. In this case we choose the following values of the incidence energy: $`E=0.001`$ eV (solid dot), and $`E=0.01`$ eV (hollow dot). In the inset of Fig. 4 we show, for the particular case of $`E=0.01`$ eV, the values of $`t_{max}`$ (solid square) at the different values of position considered in the main graph. As can be clearly appreciated in that figure, the tunneling process along the whole internal region is governed by under-barrier-frequency components, i.e., $`\omega _{av}/\omega _V<1`$. The above results should be contrasted with those obtained for very large or infinite barrier widths recently studied in Ref. gcvdm02 , where under-the-barrier tunneling components arise only at distances smaller than a characteristic length. Here we have found a particular combination of potential parameters characterizing a genuine tunneling process. The existence of under-the-barrier tunneling forerunners is not only restricted to values of the position along the internal region of the potential. In fact, we can see in Fig.4, that for distances up to $`x2L`$ along the external region, we can still observe frequency components below the cutoff-frequency, $`\omega _V`$. As $`x/L`$ increases further, $`\omega _{av}/\omega _V>1`$. This behavior follows because the corresponding values of $`t_{max}`$ increase linearly, indicating the prevalence of non-tunneling components (see the inset in Fig. 4).

We have found that the existence of tunneling forerunners at the barrier edge $`x=L`$ may be described more generally by referring to the opacity $`\alpha `$ of the system, defined as,

$$\alpha =\frac{\left[2mV\right]^{1/2}}{\mathrm{}}L,$$
(6)
and by the dimensionless parameter $`u`$, the ratio between the potential barrier height and the incidence energy,

$$u=\frac{V}{E}.$$
(7)
There is a regime, that we shall refer to as tunneling regime, characterized by the opacity $`\alpha `$ and the parameter $`u`$, where the relative frequencies associated to the time domain resonance are below the cutoff-frequency $`\omega _V`$. In order to show this, in Fig. 5 we plot the relative frequency $`\omega _{av}/\omega _V`$ as function of the opacity $`\alpha `$, for three different values of the parameter $`u`$, $`u=5`$ (solid dot), $`u=10`$ (solid triangle), and $`u=300`$ (solid square). Although in this case we have chosen a value of $`V=0.3`$ eV, we have checked that the relative frequency $`\omega _{av}/\omega _V`$, exhibits a striking regularity when plotted as a function of $`\alpha `$. For a given value of $`\alpha `$ all the systems with the same parameter $`u`$, share the same curve. That is, all systems characterized by the same parameters $`\alpha `$ and $`u`$ yield the same relative frequency, $`\omega _{av}/\omega _V`$. In view of this observed regularity, we can completely characterize the regime associated with under-the-barrier frequency components. Although in Fig. 5 we have considered values of the parameter $`u`$ such that $`5u300`$, the case corresponding to very large values of $`u`$ ($`u\mathrm{}`$), almost overlaps with the case $`u=300`$. Thus, for very large values of $`u`$ there is an upper-bound for the opacity $`\alpha _u3.3`$. This result, and the fact that the lower-bound for the opacity notelb  is given by $`\alpha _c=2.065`$, recently reported in Ref. gcv02 , allows to characterize an opacity “window” given by $`\alpha _c\alpha \alpha _u`$ where the relative frequencies are always below the cutoff-frequency $`\omega _V/\omega _V=1`$, irrespective of the value of the parameter $`u`$, namely of the value of the incidence energy. Note that the values of $`\alpha `$ within the opacity “window” may be obtained using typical parameters of semiconductor heterostructuresmendez . From the figure one can appreciate the existence of another regime of $`\alpha `$’s where the existence of under-the-barrier processes associated to the time-domain resonance may depend on the value of the parameter $`u`$.

It is of interest to remark that the time scale given by the peak maximum $`t_{max}`$ of the forerunner at the barrier edge $`x=L`$, i.e., the time domain resonance, differs from both the semi-classical Büttiker-Landauer and Büttiker traversal times, which exhibit a linear dependence with $`L`$, and also differs from the phase-timegcv01 ; gcv02 . In Ref. gcv02  we have discussed an unexpected relationship between the delay time and the existence of time domain resonances.

To conclude we remark that the analytical solution to the time-dependent Schrödinger’s equation with quantum shutter initial conditions applies in general to arbitrary potentials provided they vanish beyond a distance and can also be extended to deal with cutoff pulses as discussed in Ref. gcv02 . Our results refer to novel transient effects in time-dependent tunneling. To test our results experimentally would require to consider the detection of tunneling particles in time domain at distances close to the interaction region.

###### Acknowledgements.

The authors thank J. G. Muga for a critical reading of the manuscript and acknowledge financial support of DGAPA-UNAM under grant No. IN101301.


text:
# Compact Floor-Planning via Orderly Spanning Trees

## 1 Introduction

In VLSI chip design, floor-planning  refers to the process of, given a graph whose nodes (respectively, edges) representing functional entities (respectively, interconnections), partitioning a rectangular chip area into a set of non-overlapping rectilinear polygonal modules (each of which describes a functional entity) in such a way that the modules of adjacent nodes share a common boundary. For example, Figure 1(b) is a floor-plan of the graph in Figure 1(a).

Early stage of the floor-planning research focused on using rectangular modules as the underlying building blocks. A floor-plan using only rectangles to represent nodes is called a rectangular dual. It was shown in  that a plane triangulation $`G`$ admits a rectangular dual if and only if $`G`$ has four exterior nodes, and $`G`$ has no separating triangles. (A separating triangle, which is also known as complex triangle , is a cycle of three edges enclosing some nodes in its interior.) As for floor-planning general plane graphs, Yeap and Sarrafzadeh  showed that rectilinear modules with at most two concave corners are sufficient and necessary.

In a subsequent study of floor-planning, He  measured the complexity of a module in terms of the number of its constituent rectangles, as opposed to the number of concave corners. A module that is a union of $`k`$ or fewer disjoint rectangles is called a $`k`$-rectangular module. Since any rectilinear module with at most two concave corners can be constructed by three rectangular modules, the result of Yeap and Sarrafzadeh  implies the feasibility of floor-planning plane graphs using 3-rectangular modules. He  presented a linear-time algorithm to construct a floor-plan of a plane triangulation using only 2-rectangular modules. He’s floor-planning algorithms consists of three phases: The first phase utilizes the canonical ordering  to assign nodes on separating triangles. The second phase involves the so-called vertex expansion operation to break all separating triangles. The third phase adapts rectangular-dual algorithms  to finalize the drawing of the floor-plan. Figure 2 depicts the shapes of the 2-rectangular modules required by He’s algorithm. For convenience, these four shapes are referred to as I-module, L-module, T-module, and Z-module throughout the rest of this paper.

In this paper, we provide a “simpler” linear-time algorithm that computes “compact” floor-plans for plane triangulations. The “compactness” of the output floor-plans is an important advantage of our algorithm. Although previous work  reveals no area information, one can verify that a floor-plan using only $`O(1)`$-rectangular modules may require area $`\mathrm{\Omega }(n)\times \mathrm{\Omega }(n)`$. The output of our algorithm for an $`n`$-node plane triangulation has area no more than $`(n1)\times \frac{2n+1}{3}`$, which can be shown to be almost worst-case optimal. What “simplicity” means is two-fold:

* First, as opposed to the multiple-phase approach of , our algorithm is based upon a recent development of orderly spanning trees , which provides an extension of canonical ordering  to plane graphs not required to be triconnected and an extension for realizer  to plane graphs not required to be triangulated. Our approach bypasses the somewhat complicated rectangular-dual phase. Aside from the two applications of orderly spanning trees reported in  (namely, succinct encodings for planar graphs with efficient query support  and 2-visibility drawings for planar graphs ), our investigation here finds another interesting application of orderly spanning trees. (A similar concept called ordered stratification and its application in constructing 2-visibility drawing were independently studied by Bonichon, Le Saëc, and Mosbah .)
* Second, the floor-plan design of our algorithm is “simpler” (in comparison with ) in its own right, in the sense that I-modules, L-modules, and T-modules suffice. (Recall that Z-modules are needed by He’s algorithm .) Our result is worst-case optimal, since there is a plane triangulation that does not admit any floor-plan consisting of only I-modules and L-modules \[22, Figure 4\].

The remainder of this paper is organized as follows. Section 2 reviews the definition and property of orderly spanning tree for plane graph. Section 3 presents our linear-time floor-planning algorithm as well as its correctness proof. Section 4 provides a lower bound for the required area for floor-planning plane triangulations. Section 5 concludes the paper.

## 2 Orderly spanning tree

A plane graph is a planar graph equipped with a fixed planar embedding. The embedding of a plane graph divides the plane into a number of connected regions, each of which is called a face. The unbounded face of $`G`$ is called the exterior face, whereas the remaining faces are interior faces. $`G`$ is a plane triangulation if $`G`$ has at least three nodes and the boundary of each face, including the exterior face, of $`G`$ is a triangle. Let $`T`$ be a rooted spanning tree of a plane graph $`G`$. Two nodes are unrelated in $`T`$ if they are distinct and neither of them is an ancestor of the other in $`T`$. An edge of $`G`$ is unrelated with respect to $`T`$ if its endpoints are unrelated in $`T`$. Let $`v_1,v_2,\mathrm{},v_n`$ be the counterclockwise preordering of the nodes in $`T`$. A node $`v_i`$ is orderly in $`G`$ with respect to $`T`$ if the neighbors of $`v_i`$ in $`G`$ form the following four blocks in counterclockwise order around $`v_i`$:

1. the parent of $`v_i`$,
2. the unrelated neighbors $`v_j`$ of $`v_i`$ with $`j<i`$,
3. the children of $`v_i`$, and
4. the unrelated neighbors $`v_j`$ of $`v_i`$ with $`j>i`$,

where each block could be empty. $`T`$ is an orderly spanning tree of $`G`$ if $`v_1`$ is on the boundary of $`G`$’s exterior face, and each $`v_i,1in`$, is orderly in $`G`$ with respect to $`T`$. It is not difficult to see that if $`G`$ is a plane triangulation, then $`B_2(v_i)`$ (respectively, $`B_4(v_i)`$) is nonempty for each $`i=3,4,\mathrm{},n`$ (respectively, $`i=2,3,\mathrm{},n1`$). For each $`i=2,3,\mathrm{},n`$, let $`p(i)`$ be the index of the parent of $`v_i`$ in $`T`$. Let $`w(i)`$ denote the number of leaves in the subtree of $`T`$ rooted at $`v_i`$. Let $`\mathrm{}(i)`$ and $`r(i)`$ be the functions such that $`v_{\mathrm{}(i)}`$ (respectively, $`v_{r(i)}`$) is the last (respectively, first) neighbor of $`v_i`$ in $`B_2(v_i)`$ (respectively, $`B_4(v_i)`$) in counterclockwise order around $`v_i`$. For example, in the example shown in Figure 1(a), one can easily verify that node 3 is indeed orderly with respect to $`T`$, where $`B_1(3)=\left\{1\right\}`$, $`B_2(3)=\left\{2\right\}`$, $`B_3(3)=\{4,5\}`$, $`B_4(3)=\{6,9\}`$, $`p(3)=1`$, $`w(3)=2`$, $`\mathrm{}(3)=2`$, and $`r(3)=9`$. When $`G`$ is a plane triangulation, it is known  that for each edge $`(v_i,v_j)`$ of $`GT`$ with $`i<j`$, at least one of $`i=\mathrm{}(j)`$ and $`j=r(i)`$ holds. To be more specific, if $`i=2`$ and $`j=n`$, then both $`2=\mathrm{}(n)`$ and $`n=r(2)`$ hold; otherwise, precisely one of $`i=\mathrm{}(j)`$ and $`j=r(i)`$ holds.

The concept of orderly spanning tree for plane graphs  extends that of canonical ordering  for plane graphs not required to be triconnected and that of realizer  for plane graphs not required to be triangulated. Specifically, when $`G`$ is a plane triangulation, (i) if $`T`$ is an orderly spanning tree of $`G`$, then the counterclockwise preordering of the nodes of $`T`$ is always a canonical ordering of $`G`$, and (ii) if $`(T_1,T_2,T_n)`$ is a realizer of $`G`$, where $`T_i`$ is rooted at $`v_i`$ for each $`i=1,2,n`$, then each $`T_i`$ plus both external edges of $`G`$ incident to $`v_i`$ is an orderly spanning tree of $`G`$. Our floor-planning algorithm is based upon the following lemma.

###### Lemma 1 (see )

Given an $`n`$-node plane triangulation $`G`$, an orderly spanning tree $`T`$ of $`G`$ with at most $`\frac{2n+1}{3}`$ leaves is obtainable in $`O(n)`$ time.

## 3 Our floor-planning algorithm

A floor-plan $`F`$ of $`G`$ is a partition of a rectangle into $`n`$ non-overlapping rectangular modules $`r_1,r_2,\mathrm{},r_n`$ such that $`v_i`$ and $`v_j`$ are adjacent in $`G`$ if and only if the boundaries of $`r_i`$ and $`r_j`$ share at least one non-degenerated line segment. The size of $`F`$ is the area of the rectangle being partitioned by $`F`$ with the convention that the corners of all modules are placed on integral grid points. For example, the size of the floor-plan shown in Figure 1(b) is $`9\times 8`$. This section proves the following main theorem of the paper.

###### Theorem 1

Given an $`n`$-node plane triangulation $`G`$ with $`n3`$, a floor-plan $`F`$ of $`G`$ can be constructed in $`O(n)`$ time such that

1. $`F`$ consists of I-modules, L-modules, and T-modules only, and
2. the size of $`F`$ is bounded by $`(n1)\times \frac{2n+1}{3}`$.

Let $`T`$ be an orderly spanning tree of $`G`$, where $`v_1,v_2,\mathrm{},v_n`$ is the counterclockwise preordering of $`T`$. Our floor-planning algorithm is described as follows. Although the first two steps of our algorithm follow how Chiang et al.  obtained their 2-visibility drawing of $`G`$ with respect to $`T`$, we list them this way to make the presentation of our algorithm more self-contained.

Algorithm $`\text{FloorPlan}(G,T)`$

Produce a (vertical) visibility drawing of $`T`$ as follows: For each $`i=1,2,\mathrm{},n`$, if $`v_i`$ is a leaf of $`T`$, then draw $`v_i`$ as a unit square; otherwise, draw $`v_i`$ as a $`1\times w(i)`$ rectangle. Place each node beneath its parent such that the children of each node is placed in the same order as in $`T`$.

Turn the above visibility drawing of $`T`$ into a 2-visibility drawing of $`G`$ by stretching the nodes downward in the least necessary amount such that $`v_i`$ and $`v_j`$ are horizontally visible to each other if and only if $`(v_i,v_j)`$ is an unrelated edge of $`G`$ with respect to $`T`$. Specifically, for each $`i=3,4,\mathrm{},n`$, the $`i`$-th iteration of this step ensures the horizontal visibility between $`v_i`$ and each node in $`B_2(v_i)`$.

First, grow a horizontal branch for $`v_n`$ from boundary of $`v_n`$ visible to $`v_2`$ such that the left boundary of the horizontal branch touches $`v_2`$. Second, for each $`i=3,4,\mathrm{},n1`$, grow horizontal branches for $`v_i`$ from the boundaries of $`v_i`$ visible to $`v_{\mathrm{}(i)}`$ and $`v_{r(i)}`$ such that the left (respectively, right) boundary of the horizontal branch touches $`v_{\mathrm{}(i)}`$ (respectively, $`v_{r(i)}`$). Furthermore, when extending the boundary of $`v_i`$, we also extend the boundaries of the descendants of $`v_i`$ to maintain the property that the bottom boundary of each internal node of $`T`$ is completely occupied by the top boundaries of its children. Note that some former extended modules might be covered by latter extending.

For each $`i=n1,n2,\mathrm{},3`$, if $`v_i`$ has a horizontal branch with height greater than one, then reduce the height of the thick branch down to one.

Pictures of intermediate steps are shown to illustrate how our algorithm obtains the floor-plan in Figure 1(b) for the plane graph $`G`$ with respect to the orderly spanning tree $`T`$ shown in Figure 1(a). Figure 3 shows how Step 1 obtains the visibility drawing for $`T`$. Figure 4 shows how Step 2 obtains the resulting 2-visibility drawing for $`G`$. Observe that the resulting drawing satisfies the property that the bottom boundary of each internal node of $`T`$ is completely occupied by the top boundaries of its children. Figure 5 illustrates how Step 3 obtains the resulting drawing for $`G`$. Note that when the horizontal branch of node 3 is extended to the right by one unit to touch the left boundary of node 9, the right boundary of node 5 is also extended to the right by the same amount. To see the necessity of Step 4, one can verify that the module for node 10 in Figure 5(d) has a thick horizontal branch. The height of this thick branch can be reduced by moving down the top boundary of the thick branch that is adjacent to the bottom boundary of node 11. The resulting floor-plan consists of only I-modules, L-modules, and T-modules. Moreover, each horizontal branch of the L-modules and T-modules has height exactly one.

###### Lemma 2

The following statements hold for our algorithm FloorPlan.

1. The algorithm can be implemented to run in $`O(n)`$ time.
2. The output is a floor-plan of $`G`$ of size no more than $`(n1)\times w(v_1)`$.
3. The resulting floor-plan consists of I-modules, L-modules, and T-modules, where the height of each horizontal branch of L-modules and T-modules is one.

Proof.  Statement 1. One can verify that our algorithm is implementable to run in linear time as follows.

Since $`w(v_1),w(v_2),\mathrm{},w(v_n)`$ can be computed from $`T`$ in $`O(n)`$ time, the described (vertical) visibility drawing of $`T`$ can easily be computed in $`O(n)`$ time.

Note that we have to ensure that $`v_i`$ and $`v_j`$ are horizontally visible to each other if and only if $`v_jB_2(v_i)`$ at the end of the stretch-down iteration for $`v_i`$. Therefore, when the boundaries of $`v_i`$ and the nodes in $`B_2(v_i)`$ are stretched down, the boundaries of some other nodes might require being stretched down as well. For example, when we obtain Figure 4(c) from Figure 4(b) by stretching down the boundary of node 6 to ensure that nodes 6 and 8 are horizontally visible to each other, we also have to increase the the heights of nodes 2 and 3 by one. Thus, a naive implementation of this step may require $`\mathrm{\Omega }(n^2)`$ time. However, this step can be implemented by directly computing the position $`y(i)`$ of the bottom boundary of $`v_i`$ for each $`i=1,2,\mathrm{},n`$ and the position $`y(i,j)`$ of the bottom boundary of each unrelated edge $`(v_i,v_j)`$ with $`i<j`$ according to the following recurrence relation:

$`y(i)`$ $`=`$ $`\{\begin{array}{cc}1\hfill & \text{if }i=1;\hfill \\ \mathrm{max}\{y(\mathrm{}(i),i),y(i,r(i))\}\hfill & \text{otherwise};\hfill \end{array}`$
$`y(i,j)`$ $`=`$ $`1+\mathrm{max}\{y_{\mathrm{}}(i,j),y_r(i,j)\},`$
where $`y_{\mathrm{}}(i,j)`$ and $`y_r(i,j)`$ are defined as follows. Let $`v_j^{}`$ be the neighbor of $`v_i`$ that immediately succeeds $`v_j`$ in counterclockwise order around $`v_i`$. Let $`v_i^{}`$ be the neighbor of $`v_j`$ that immediately precedes $`v_i`$ in counterclockwise order around $`v_j`$. By $`i<j`$, one can easily see that either $`i^{}=p(j)`$ or $`v_i^{}B_2(v_j)`$ holds. Similarly, either $`j^{}=p(i)`$ or $`v_j^{}B_4(v_i)`$ holds. Let

$`y_{\mathrm{}}(i,j)`$ $`=`$ $`\{\begin{array}{cc}y(j^{})\hfill & \text{if }j^{}=p(i);\hfill \\ y(i,j^{})\hfill & \text{otherwise};\hfill \end{array}`$
$`y_r(i,j)`$ $`=`$ $`\{\begin{array}{cc}y(i^{})\hfill & \text{if }i^{}=p(j);\hfill \\ y(i^{},j)\hfill & \text{otherwise}.\hfill \end{array}`$
Clearly, the bottom positions $`y(i)`$ of all nodes $`v_i`$ can be obtained in $`O(n)`$ time by dynamic programming. Since the top position of $`v_i`$ is simply $`y(p(i))`$, the resulting 2-visibility drawing of $`G`$ can be obtained in $`O(n)`$ time.

On the one hand, a naive implementation of this step may require $`\mathrm{\Omega }(n^2)`$ time, since growing the horizontal branches for a node may cause boundary extension for its descendants. On the other hand, although in the $`i`$-th iteration we are supposed to extend the boundary of some descendants $`v_j`$ of $`v_i`$, we do not need to actually extend the boundaries of $`v_j`$ until the beginning of the $`j`$-th iteration. Note that how far should the boundary of $`v_j`$ be extended can be determined directly from the boundary of $`v_{p(j)}`$ in the $`j`$-th iteration. Clearly, the above “lazy” strategy reduces the time complexity of this step to $`O(n)`$. Since the unrelated edge $`(v_i,v_j)`$ of $`GT`$ with $`i<j`$ and $`(v_i,v_j)(v_2,v_n)`$ satisfies exactly one equality of $`i=\mathrm{}(j)`$ and $`j=r(i)`$, the resulting drawing is a partition of a rectangle into $`n`$ rectilinear regions. (That is, there is no gap among modules in the rectangle.) To prove that the resulting drawing is indeed a floor-plan of $`G`$, it suffices to show that growing a horizontal branch of $`v_i`$ is to reach the boundary of $`v_j`$ does not result in new adjacency among these rectilinear modules. Suppose $`v_k`$ is a node whose bottom boundary touches the top bottom of the horizontal branch of $`v_i`$. Assume for a contradiction that $`v_k`$ is not adjacent to $`v_i`$ in $`G`$. Since the resulting drawing of the previous step is a 2-visibility drawing of $`G`$, there must be a node $`v_k^{}`$ lies between $`v_i`$ and $`v_k`$ preventing their horizontal visibility to each other. It follows that there is a face of $`G`$ containing at least four nodes $`v_i,v_j,v_k,v_k^{}`$, contradicting the fact that $`G`$ is triangulated.

Since $`T`$ is an orderly spanning tree of $`G`$ and $`G`$ is a plane triangulation, one can see that if $`v_i`$ grows a horizontal branch to reach $`v_j`$, then there must be a unique node $`v_k`$ whose bottom boundary touches the top boundary of that horizontal branch of $`v_i`$. It is also not difficult to verify that both $`(v_i,v_k)`$ and $`(v_j,v_k)`$ are unrelated edges $`G`$ with respect to $`T`$. Thus, in the resulting drawing of the previous step, the left and right boundaries of $`v_k`$ have to touch $`v_i`$ and $`v_j`$. Therefore, the height of that horizontal branch of $`v_i`$ can be reduced to one by moving downward the bottom boundary of $`v_k`$, which is also the top boundary of that horizontal branch, without changing the adjacency of $`v_k`$ to other nodes in the floor-plan. Clearly, each height-reducing operation takes $`O(1)`$ time by adapting lazy strategy, so this step runs in $`O(n)`$ time. Since the for-loop of this step proceeds from $`i=n1`$ down to $`3`$, each horizontal branch has height exactly one at the end of this step.

Statement 2. Since Steps 3 and 4 do not affect the adjacency among the rectilinear modules, it suffices to ensure that the 2-visibility drawing obtained in Step 2 has size no more than $`(n1)\times w(v_1)`$. By the definition of Steps 1 and 2, it is straightforward to see that the width of the resulting drawing is precisely $`w(v_1)`$. It remains to show that $`y(2,n)`$, which is exactly the height of the resulting 2-visibility drawing, is no more than $`n1`$ as follows. Assume for a contradiction that $`y(2,n)n`$. It follows that there is a sequence of unrelated edges $`(v_{s_1},v_{t_1}),(v_{s_2},v_{t_2}),\mathrm{},(v_{s_n},v_{t_n})`$ with

$$2=s_ns_{n1}s_{n2}\mathrm{}s_1<t_1t_2\mathrm{}t_n=n$$

such that at least one of $`s_is_{i+1}`$ and $`t_it_{i+1}`$ holds for each $`i=1,2,\mathrm{},n1`$. It follows that the set $`\{s_1,s_2,\mathrm{},s_n,t_1,t_2,\mathrm{},t_n\}`$ contains at least $`n`$ distinct integers, thereby, contradicting the assumption $`2s_i,t_in`$.

Statement 3. By the definition of Step 3, one can easily verify that the resulting floor-plan consists of I-modules, L-modules, and T-modules. By the height-reducing operation performed on the horizontal branches in Step 4, the statement is proved.

We are ready to prove the main theorem as follows.

Proof. \[Proof for Theorem 1\] Straightforward by Lemmas 1 and 2.

## 4 Lower bounds on the worst-case area of floor-plan

This section shows the near optimality of the output of our algorithm.

###### Lemma 3

For each $`n3`$, there is an $`n`$-node plane triangulation graph $`G_n`$ such that any $`h_n\times w_n`$ floor-plan of $`G_n`$ satisfies $`\mathrm{min}\{h_n,w_n\}\frac{2n+1}{3}`$ and $`h_n+w_n\frac{4n}{3}`$.

Proof.  The lower-bound examples are constructed inductively: For each $`n4`$, $`G_n`$ is obtained from $`G_{n3}`$ by adding an external triangle and arbitrarily triangulating the face between the external triangle of $`G_n`$ and the external boundary of $`G_{n3}`$. As for the base cases, let $`G_n`$ be an arbitrary $`n`$-node plane triangulation for each $`n=3,4,5`$. Now we show that the required inequalities hold for each $`n3`$. As for the inductive basis, one can verify $`\mathrm{min}\{h_3,w_3\}2`$, $`h_3+w_34`$, $`\mathrm{min}\{h_4,w_4\}3`$, $`h_4+w_46`$, $`\mathrm{min}\{h_5,w_5\}3`$, and $`h_5+w_57`$. Therefore the inequalities hold for the base cases. It remains to ensure the induction step as follows.

$`\mathrm{min}\{h_n,w_n\}`$ $``$ $`\mathrm{min}\{h_{n3},w_{n3}\}+2`$
$``$ $`{\displaystyle \frac{2(n3)+1}{3}}+2`$
$`=`$ $`{\displaystyle \frac{2n+1}{3}};`$
$`h_n+w_n`$ $``$ $`h_{n3}+w_{n3}+4`$
$``$ $`{\displaystyle \frac{4(n3)}{3}}+4`$
$`=`$ $`{\displaystyle \frac{4n}{3}}.`$
## 5 Conclusion

A linear-time algorithm for producing compact floor-plans for plane triangulations has been designed. Our algorithm is based upon a newly developed technique of orderly spanning trees with bounded number of leaves . In comparison with previous work on floor-planning plane triangulations , our algorithm is simpler in the algorithm itself as well as in the resulting floor-plan in the sense that the Z-modules required by  is not needed in our design. Another important feature of our algorithm is the upper bound $`(n1)\times \frac{2n+1}{3}`$ on the area of the output floor-plan. Previous work  does not provide any area bounds on their outputs. Investigating whether the $`(n1)\times \frac{2n+1}{3}`$ area is worst-case optimal is an interesting future research direction.

## Acknowledgment

We thank the anonymous referees for their helpful comments, which significantly improve the presentation of the paper. We also thank Ho-Lin Chen for his comments regarding an early version of this work.