text:
# Untitled Document

An Investigation of Time in Relativity

Sidi Cherkawi Benzahra111sbenzahra@wittenberg.edu or benzahra@physics.umn.edu

Wittenberg University

Springfield, Ohio 45501

## Abstract

I added an imaginary term by hand to the proper time, $`\mathrm{\Delta }T`$, of special theory of relativity and wanted to see what will happen to it if I tried to get the time dilation relationship between $`\mathrm{\Delta }T^{}`$ and $`\mathrm{\Delta }T`$. I found that this added imaginary term of time does not have to be zero.

If $`\mathrm{\Delta }T`$ is the real, proper time of travel of a light beam from one point to the other,

$$\mathrm{\Delta }t=\mathrm{\Delta }T+i\mathrm{\Delta }\xi $$
(1)
$`\mathrm{\Delta }t`$ will be the time I choose and work with in order to see what will happen to the imaginary part $`\mathrm{\Delta }\xi `$ once I apply Einstein’s special relativity. $`i\mathrm{\Delta }\xi `$ here is an imaginary term which I added by hand. Consider two observers $`O`$ and $`O^{}`$ . $`O`$ fires a beam of light at a mirror a distance L away and measures the time interval 2 $`\mathrm{\Delta }T`$ for the beam to be reflected from the mirror and to return to $`O`$. Observer $`O^{}`$ is moving at a constant velocity $`u`$. As seen from the point of view of $`O`$, the beam is sent and received from the same point, while $`O^{}`$ moves off in a perpendicular direction. The beam is sent from some point $`A`$ and received at some point $`B`$ a time 2 $`\mathrm{\Delta }T^{}`$ later, according to $`O^{}`$. The distance AB is just $`2u\mathrm{\Delta }T^{}`$. According to $`O`$, the light beam travels a distance $`2L`$ in a time $`2\mathrm{\Delta }T`$. According to $`O^{}`$, the light beam travels a distance of $`2\sqrt{L^2+(u\mathrm{\Delta }T^{})^2}`$ in a time $`2\mathrm{\Delta }T^{}`$. According to Galilean relativity, $`\mathrm{\Delta }T=\mathrm{\Delta }T^{}`$, and $`O`$ measures a speed $`c`$ while $`O^{}`$ measures a speed $`\sqrt{c^2+u^2}`$. According to Einstein’s second postulate, this is not possible–both $`O`$ and $`O^{}`$ must measure the speed $`c`$. Therefore $`\mathrm{\Delta }T`$ and $`\mathrm{\Delta }T^{}`$ must be different. I can find a relationship between $`\mathrm{\Delta }T`$ and $`\mathrm{\Delta }T^{}`$ by setting the two speeds equal to $`c`$. According to $`O`$, $`c=2L/2\mathrm{\Delta }T`$, so $`L=c\mathrm{\Delta }T`$. According to $`O^{}`$, $`c=2\sqrt{L^2+(u\mathrm{\Delta }T^{})^2}/2\mathrm{\Delta }T^{}`$, so $`c\mathrm{\Delta }T^{}=\sqrt{L^2+(u\mathrm{\Delta }T^{})^2}`$. Combining these, I find

$$c\mathrm{\Delta }T^{}=\sqrt{(c\mathrm{\Delta }T)^2+(u\mathrm{\Delta }T^{})^2}$$
(2)
and, solving for $`\mathrm{\Delta }T^{}`$,

$$\mathrm{\Delta }T^{}=\frac{\mathrm{\Delta }T}{\sqrt{1u^2/c^2}}.$$
(3)
This relationship summarizes the effect known as time dilation. But in my calculation, where $`\mathrm{\Delta }t=\mathrm{\Delta }T+i\mathrm{\Delta }\xi `$, I have according to $`O`$,

$$c=\frac{L}{\mathrm{\Delta }T+i\mathrm{\Delta }\xi }.$$
(4)
I can see here in the equation just above that since I added an imaginary part to time, length in particular or space in general picked up an imaginary part also. According to $`O^{}`$

$$c=\frac{\sqrt{L^2+[u(\mathrm{\Delta }T^{}+i\mathrm{\Delta }\xi ^{})]^2}}{\mathrm{\Delta }T^{}+i\mathrm{\Delta }\xi ^{}}.$$
(5)
Taking the expression of $`L`$ from equation (4) and inserting it in equation (5) and squaring, I get

$$c^2(\mathrm{\Delta }T^{}+i\mathrm{\Delta }\xi ^{})^2=c^2(\mathrm{\Delta }T+i\mathrm{\Delta }\xi )^2+u^2(\mathrm{\Delta }T^{}+i\mathrm{\Delta }\xi ^{})^2.$$
(6)
Separating the real terms from the imaginary terms of equation (6) I get

$$[\mathrm{\Delta }T_{}^{}{}_{}{}^{2}(1u^2/c^2)\mathrm{\Delta }T^2][\mathrm{\Delta }\xi _{}^{}{}_{}{}^{2}(1u^2/c^2)\mathrm{\Delta }\xi ^2]+2i[\mathrm{\Delta }T^{}\mathrm{\Delta }\xi ^{}(1u^2/c^2)\mathrm{\Delta }T\mathrm{\Delta }\xi ]=0$$

I can see that time dilation relationship is given by setting the first part of the equation above equal to zero

$$\mathrm{\Delta }T^{}=\frac{\mathrm{\Delta }T}{\sqrt{1u^2/c^2}}.$$
(7)
I can also see that there is an effect of time dilation in the imaginary part, $`\mathrm{\Delta }\xi `$, given by

$$\mathrm{\Delta }\xi ^{}=\frac{\mathrm{\Delta }\xi }{\sqrt{1u^2/c^2}}.$$
(8)
Since a complex number vanishes only if its real and imaginary parts vanish too, I get

$$\mathrm{\Delta }T^{}\mathrm{\Delta }\xi ^{}(1u^2/c^2)\mathrm{\Delta }T\mathrm{\Delta }\xi =0.$$
(9)
Taking the right-hand side of equation (8) and inserting it in equation (9), I get

$$[\mathrm{\Delta }T^{}\sqrt{1u^2/c^2}\mathrm{\Delta }T]\mathrm{\Delta }\xi =0.$$
(10)
Equation (10) is the product of two terms. Since Einstein’s relativity is considered to be true, the first part vanishes, and leaves us with the condition that $`\mathrm{\Delta }\xi `$ can be different from zero. And if $`\mathrm{\Delta }\xi `$ is different from zero that makes $`\mathrm{\Delta }t`$ imaginary since $`\mathrm{\Delta }t`$ is equal to $`\mathrm{\Delta }T+i\mathrm{\Delta }\xi `$ as indicated in equation (1). I have mathematically found that relativity allows time to have an imaginary part.

References

 Kenneth Krane, Modern Physics (John Wiley and Sons, 1983).


format: text

text:
# Lee-Yang zeros and phase transitions in nonequilibrium steady states

(18th April, 2002; revised 19th June, 2002)

## Abstract

We consider how the Lee-Yang description of phase transitions in terms of partition function zeros applies to nonequilibrium systems. Here one does not have a partition function, instead we consider the zeros of a steady-state normalization factor in the complex plane of the transition rates. We obtain the exact distribution of zeros in the thermodynamic limit for a specific model, the boundary-driven asymmetric simple exclusion process. We show that the distributions of zeros at the first and second order nonequilibrium phase transitions of this model follow the patterns known in the Lee-Yang equilibrium theory.

With equilibrium statistical physics now firmly established as a successful general theory for predicting the macroscopic behavior of many-body systems, attention has more recently turned to the the wider class of nonequilibrium systems. Examples of nonequilibrium systems are those relaxing to thermal equilibrium or those driven into a steady state, far from thermal equilibrium. A general framework for understanding these systems is proving elusive. However, over the last two decades the detailed study of specific models has revealed the wide range of phenomena that may emerge.

It is by now well known that models of driven diffusive systems can exhibit phase transitions in their nonequilibrium steady states Schmittmann and Zia (1995). Typically such models comprise one or more species of particles being driven by an external field along with prescribed interactions between the boundaries and the outside world. Applications of these models are diverse and include the kinetics of biopolymerization MacDonald et al. (1968), transport across a membrane Chou and Lohse (1999) and traffic flow Chowdhury et al. (2000). From a more fundamental viewpoint, the interest lies in the richness of collective phenomena displayed, including jamming O’Loan et al. (1998) and spontaneous symmetry breaking Evans et al. (1995), even in one dimension.

Half a century ago Lee and Yang provided a theory of equilibrium phase transitions based around the zeros of the partition function. In this work we focus on nonequilibrium steady states and how the phase transitions they admit may be placed in the context of the Lee-Yang theory. In the original exposition Yang and Lee (1952) the fugacity of a gas was generalized to the complex plane and the zeros of the grand-canonical partition function in this plane were considered. In particular, it was shown that the free energy is analytic in any region of the complex-fugacity plane devoid of any zeros. Conversely, there are nonanalyticities at points where the density of partition function zeros accumulate in the thermodynamic limit. These accumulation points correspond to physically observable phase transitions if they lie on the positive real axis. A similar phenomenon occurs in the complex plane of other intensive fugacity-like variables. For example, Lee and Yang also showed Lee and Yang (1952) that the zeros of the Ising model partition in the complex plane of activity (defined as $`\mathrm{exp}(h/T)`$ where $`h`$ is magnetic field and $`T`$ temperature) lie on the unit circle. Similarly, partition function zeros in the complex temperature plane also accumulate at phase transition points Fisher (1965); Grossmann and Rosenhauer (1969). The applicability of the Lee-Yang theory to equilibrium transitions driven by intensive field-like variables can be understood from the fact that, mathematically, these variables play similar roles in the partition function.

The difficulty in applying the Lee-Yang theory to, e.g., the steady state of a driven diffusive system is that one does not have to hand a partition function defined in terms of thermodynamic state variables. However, one can always define a quantity $`Z`$ as a sum of the steady state weights (unnormalized probabilities) $`f(𝒞)`$:

$$Z=\underset{𝒞}{}f(𝒞).$$
(1)
Thus the probability of a configuration is given by $`P(𝒞)=f(𝒞)/Z`$. In this work, we treat the normalization $`Z`$ as a nonequilibrium analog of the partition function.

In order to obtain $`Z`$ one still has to calculate the steady state weights. Up to a multiplicative factor they are implied by the transition rates $`W(𝒞𝒞^{})`$ that define the model through the requirement that the total inflow of probability into a given configuration $`𝒞`$ must be balanced by the outflow, i.e.,

$$\underset{𝒞^{}𝒞}{}\left[f(𝒞^{})W(𝒞^{}𝒞)f(𝒞)W(𝒞𝒞^{})\right]=0.$$
(2)
The solution of this set of equations for the steady-state weights $`f(𝒞)`$ can, in principle, be obtained using standard methods, such as Gaussian elimination, Cramer’s rule or graphically Evans and Blythe (2001). If one applies such a method, the weights can always be written as polynomials of the elementary transition rates. Thus, for a finite system, $`Z`$ is also a polynomial of the transition rates.

In this work we generalize the transition rates to the complex plane and in doing so consider the zeros of the partition function. In analogy with the Lee-Yang picture for equilibrium systems, for a phase transition to arise we expect the zeros to pinch the real axis at a point which corresponds to a (real) critical transition rate.

As an aside, we note that the general methods of solving (2) mentioned above are rarely tractable when the number of configurations becomes large. One exception is when the system satisfies detailed balance. In this case, each term within the square brackets in (2) is identically zero and it is easy to construct the set of steady-state weights as Boltzmann weights with respect to an energy function. The resulting $`Z`$ is then an equilibrium partition function. A simple example of such a system is the Ising model evolving under Glauber dynamics wherein studying the zeros of the partition function as a function of the transition rate would be equivalent to studying the complex-temperature zeros.

In the present work, however, we are interested in models which do not satisfy detailed balance and hence for which the normalization $`Z`$ is not a priori known. This makes it difficult to make general statements about the significance of the normalization $`Z`$. We have argued that a nonequilibrium phase transition will be signaled by an accumulation of the zeros of $`Z`$ in the complex plane of transition rates towards the real axis. We demonstrate this in the case of a particular nonequilibrium model for which an expression for $`Z`$ is exactly known and we can calcuate the distribution of zeros analytically. We find that these zeros lie on circles in the complex plane of an appropriately chosen function of a transition rate that plays the role of a fugacity. The characteristic properties of the distributions of partition function zeros at first and second order equilibrium transitions Yang and Lee (1952); Grossmann and Rosenhauer (1969); Grossmann and Lehmann (1969) are also exhibited in this nonequilibrium case. In short, we show how the the Lee-Yang theory of phase transitions generalizes to a nonequilibrium steady state.

Before discussing the model and our results in greater detail, we wish to remark how our study differs from a related work due to Arndt Arndt (2000). The aim of that investigation was to locate a phase transition induced by varying the relative numbers of different particle species in a particle-conserving driven diffusive system. To this end, a fugacity-like quantity was introduced in an ad hoc way. Physically, this corresponds to placing the driven diffusive system in contact with a particle reservoir thus introducing equilibrium fluctuations in the particle densities. We stress that in the present work, it is elementary transition rates, which need not have any connection to intensive field (fugacity-like) variables, that are generalized to the complex plane.

We now introduce the model we consider, namely the one-dimensional asymmetric simple exclusion process (ASEP) with open boundaries. Since its introduction MacDonald et al. (1968); Krug (1991) this model has received a great deal of attention—see, e.g., Schütz (2000); Evans and Blythe (2001) for further details and references. The model comprises a lattice of $`N`$ sites, each of which may be occupied or empty. In an infinitesimal time interval $`\mathrm{d}t`$ one of the following transitions may occur: a particle may hop one site to the right with probability $`\mathrm{d}t`$ (i.e. at unit rate) subject to the receiving site being vacant; a particle may be inserted onto the leftmost site (if vacant) with probability $`\alpha \mathrm{d}t`$ or a particle occupying the rightmost site (if present) may be removed with probability $`\beta \mathrm{d}t`$. Fig. 1 illustrates the microscopic dynamics.

As outlined above, there is a current (defined as the flux of particles between neighbouring sites per unit time) in the steady state which in the thermodynamic limit $`N\mathrm{}`$ exhibits nonanalyticities as the boundary rates $`\alpha `$ and $`\beta `$ are varied. These are associated with phase transitions, and the phase diagram is given in Fig. 2.

Phase (i) is a high-density phase with a density profile that decays exponentially towards the right boundary. In this phase the current $`J=\beta (1\beta )`$. Since the ASEP is invariant under particle-hole exchange and the swap $`\alpha \beta `$, it follows that phase (ii) is a low-density phase that has an exponential decay in the density profile from the left boundary and a current $`J=\alpha (1\alpha )`$. At the first order transition line $`\alpha =\beta <\frac{1}{2}`$ the current exhibits a discontinuity in its first derivative. Here the system exhibits a shock front separating regions of high and low density; this is an example of phase coexistence at a nonequilibrium phase transition.

In phase (iii), the current assumes a constant value $`J=\frac{1}{4}`$ which is the largest possible current admitted for any combination of $`\alpha `$ and $`\beta `$. Hence phase (iii) is called the maximal current phase. This phase has the interesting feature that the density profile decays as a power-law from both boundaries and hence at the transition from either phase (i) or (ii) to phase (iii) a correlation length diverges. Additionally, $`J`$ has a discontinuity in its second derivative and so this is a second-order transition.

These results are known mainly from the exact solution Derrida et al. (1993); Schütz and Domany (1993). The two key results we call upon here are that the normalization for an $`N`$-site system $`Z_N`$ is given by

$$Z_N=\underset{p=1}{\overset{N}{}}\frac{p(2N1p)!}{N!(Np)!}\frac{(1/\beta )^{p+1}(1/\alpha )^{p+1}}{(1/\beta )(1/\alpha )}$$
(3)
and the current by $`J=Z_{N1}/Z_N`$. In order later to construct the distribution of the zeros of $`Z_N`$ for the ASEP we also need the large-$`N`$ forms of $`Z_N`$. From the exact solution it can be shown that for all $`\alpha ,\beta `$, $`Z_NAJ^NN^\gamma `$ where $`J`$ is the current and $`A,\gamma `$ depend only on the phase being considered. We see $`\mathrm{ln}J`$ is the extensive part of $`\mathrm{ln}Z_N`$ and thus $`\mathrm{ln}J`$ plays the role of the free energy and $`J`$ the fugacity for this nonequilibrium system.

We now consider the zeros of $`Z_N`$ given by (3) in the complex-$`\alpha `$ plane at fixed $`\beta `$ Derrida (2002). (The symmetry in $`\alpha `$ and $`\beta `$ of $`Z_N`$ implies that one would obtain the same pattern of zeros in the complex-$`\beta `$ plane at fixed $`\alpha `$). Using Mathematica, we obtained numerical estimates of the zeros of $`Z_N`$ for system size (and number of zeros) $`N=300`$ and $`\beta =1,\frac{1}{3}`$. The results are shown in Fig. 3.

For the case $`\beta =\frac{1}{3}`$, the zeros appear to approach the positive real $`\alpha `$ axis at a value $`\alpha _c=\frac{1}{3}`$, coincident with the first-order phase transition point in this regime. The plot for the case $`\beta =1`$ is suggestive of a slow accumulation of the zeros to the positive real $`\alpha `$ axis at $`\alpha _c=\frac{1}{2}`$, i.e. the second-order phase transition point. It appears also in this case that the zeros approach the real axis at an angle $`\frac{\pi }{4}`$ as opposed to an angle $`\frac{\pi }{2}`$ at the first-order transition. We will shortly demonstrate that the above assertions—consistent with the patterns of zeros at first- and second-order phase transitions known from the equilibrium Lee-Yang theory Yang and Lee (1952); Grossmann and Rosenhauer (1969); Grossmann and Lehmann (1969)—are correct.

It is useful to make a change of variable from $`\alpha `$ to $`\xi =\alpha (1\alpha )`$ (note that in the low density phase $`\xi =J`$) and we plot the numerically-obtained zeros in the complex-$`\xi `$ plane in Fig. 4. For both $`\beta =1,\frac{1}{3}`$, it appears that, as $`N\mathrm{}`$, the zeros become uniformly distributed around a circle of radius $`\xi _c=\alpha _c(1\alpha _c)`$ (these circles are shown as solid lines in the figure). We shall now show how to obtain this result analytically.

We first note (for a proof see Derrida (1991)) that for any polynomial $`p_N(\xi )`$ with $`N`$ zeros, its density of zeros $`\rho (x,y)`$ in the complex $`\xi =x+iy`$ plane is given by $`\rho (x,y)=2\pi ^2\mathrm{ln}|p_N(\xi )|`$ where $`^2\frac{^2}{x^2}+\frac{^2}{y^2}`$. In the thermodynamic limit ($`N\mathrm{}`$) we expect a continuous distribution of zeros and wish to normalize $`\rho `$ so that its integral over all space is unity. We thus define the key quantity

$$\varphi =\underset{N\mathrm{}}{lim}\frac{\mathrm{ln}|Z_N|}{N}$$
(4)
so that $`\rho =2\pi ^2\varphi `$. In equilibrium statistical physics $`\varphi `$ would be interpreted as the extensive part of the free energy; as noted above $`\varphi =\mathrm{ln}|J|`$ for the ASEP.

It is useful to regard $`\varphi `$ as an ‘electrostatic’ potential in a two-dimensional space. We introduce the field $`\stackrel{}{E}=\varphi `$ and consider the discontinuities that arise as a line of charge (zeros) is crossed. The appropriate forms of Gauss’s and Ampere’s laws imply that the component of $`\stackrel{}{E}`$ parallel to a line of zeros separating two phases must be continuous across it whereas the perpendicular component differs by a value $`2\pi g(x,y)`$ where $`g(x,y)`$ is the line density of charge (zeros) along the boundary.

The current itself has only two analytic forms. In the interior region of Fig. 4 (i.e., the region that includes $`\xi =0`$) $`J=\alpha (1\alpha )=\xi `$ and hence $`\varphi =\mathrm{ln}|\xi |`$. Over the exterior region, the current is constant: $`J=\alpha _c(1\alpha _c)=\xi _c`$ and so $`\varphi =\mathrm{ln}|\xi _c|`$.

The field that results from the potential $`\varphi `$ is zero in the exterior region and directed radially towards the origin in the interior region; specifically $`\stackrel{}{E}=\stackrel{}{r}/r^2`$ where $`\stackrel{}{r}=(\mathrm{Re}\xi ,\mathrm{Im}\xi )`$. Now, since the component of the electric field parallel to the boundary between the two phases varies continuously and is zero in the exterior region, the phase boundary must lie perpendicular to the field lines in the interior region. In other words, the phase boundary is an equipotential $`\varphi =\mathrm{ln}|\xi |=\text{const}|\xi |=\xi _c`$ since the boundary must pass through the transition point at $`\xi =\xi _c`$. The line-density of zeros $`g(\theta )`$ along the circle $`\xi =\xi _ce^{i\theta }`$ can be calculated by noting that the perpendicular component of the field at the boundary has a magnitude $`1/|\xi _c|`$. The relationship given above then implies $`g(\theta )=(2\pi \xi _c)^1`$, i.e., a constant as claimed.

We now use this result to investigate the distribution of zeros near the transition point in the complex-$`\alpha `$ plane. In the $`\xi `$ plane, the circle of zeros passes along the line $`\mathrm{Re}\xi =\xi _c`$ which maps onto a curve in the complex $`\alpha =u+iv`$ plane that has $`(u\frac{1}{2})^2v^2=\frac{1}{4}\xi _c`$. The solution of this equation that describes a line passing through the transition point $`\alpha =\alpha _c`$ is $`u=\frac{1}{2}(v^2+\frac{1}{4}\xi _c)^{1/2}`$.

When the phase transition is first order, $`\beta <\frac{1}{2}`$ and $`\xi _c=\beta (1\beta )<\frac{1}{4}`$. Then, the curve of zeros passes smoothly through the transition point $`\alpha =\beta `$. One can show that the line density of zeros at this point in the complex-$`\alpha `$ plane is $`(12\beta )/[2\pi \beta (1\beta )]`$, i.e., nonzero. The fact that the density of zeros is nonzero on the real axis at a first-order phase transition is well-known in the equilibrium Lee-Yang theory Yang and Lee (1952).

At the second order phase transition ($`\beta \frac{1}{2}`$ and $`\xi _c=\frac{1}{4}`$) the line of zeros is given by $`u=\frac{1}{2}|v|`$. This means the zeros approach the transition point $`\alpha =1/2`$ at an angle of $`\frac{\pi }{4}`$ to the real axis, meeting at a right angle. If one defines $`\mathrm{}`$ as the displacement along the straight line of zeros from the transition point, one finds that it varies with $`y=\mathrm{Im}\xi `$ as $`\mathrm{}=\sqrt{y}`$. This implies that the density of zeros behaves as $`g(\mathrm{})=\frac{2}{\pi }\frac{\mathrm{d}y}{\mathrm{d}\mathrm{}}=\frac{4\mathrm{}}{\pi }`$. In the equilibrium theory, it is well known that the density of zeros vanishes as a power-law towards the transition point on the real line (see, e.g., Grossmann and Lehmann (1969)), a result we have recovered in this nonequilibrium case.

To summarize, we have shown that the zeros for the normalization of the ASEP accumulate at the phase transition points in the complex plane of transition rates. As in the case of the Lee-Yang theory, a first-order transition is manifested by a nonzero density of zeros at the accumulation point, whereas the density of zeros decays as a power-law towards a continuous transition point.

A point that we have so far ignored is that the normalization (1) is only defined up to a multiplicative factor that depends on the method used to solve Eqn. (2). This additional factor is analogous to that which would appear in an equilibrium partition function after a uniform shift of the energy scale. Although this factor could itself be a polynomial of the transition rates and so introduce additional zeros, one would not expect these spurious zeros to be physically relevant to the phase behavior. For example one way to define the normalization is to use Cramer’s rule to solve (2) (see Evans and Blythe (2001)). Then for the ASEP the normalization is given by $`(\alpha \beta )^N`$ times the expression (3) and additional zeros are introduced only at the origin.

Finally we remark on the applicability of the Lee-Yang theory to other nonequilibrium systems. Firstly, a generalization of the ASEP that allows particles also to hop to the left at a rate $`q`$ has been solved Sasamoto (1999). The phase behavior when $`q<1`$ is very similar to that of the ASEP, and the rate $`q`$ enters into the expressions for the current in a simple way. The result of this would be for the pattern of zeros in the $`\alpha `$-plane to shrink towards the origin as $`q1`$. At this point, there is a transition to a regime where the current ceases to flow in the thermodynamic limit. It would be of interest to see whether the zeros of the normalization in the complex-$`q`$ plane accumulate at $`q=1`$ since this would provide some evidence for the generality of the Lee-Yang picture when applied to complex transition rates. Further evidence is provided by a study of the percolation probability on finite directed percolation lattices Dammer et al. (2002) in which a similar accumulation phenomenon to that reported here was observed. Also, one might learn something of the approach to the thermodynamic limit by considering how the distribution of zeros varies as the system size is increased.

###### Acknowledgements.

We thank B. Derrida for insightful discussion. R.A.B. thanks the EPSRC for financial support under grant GR/R53197.


text:
# Computing Homotopic Shortest Paths Efficiently

## 0 Introduction

Finding Euclidean shortest paths in simple polygons is a well-studied problem. The funnel algorithm of Chazelle  and Lee and Preparata  finds the shortest path between two points in a simple polygon. Hershberger and Snoeyink  unify earlier results for computing shortest paths in polygons. They optimize a given path among obstacles in the plane under the Euclidean and link metrics and under polygonal convex distance functions. Related work has been done in addressing a classic VLSI problem, the continuous homotopic routing problem . For this problem it is required to route wires with fixed terminals among fixed obstacles when a sketch of the wires is given, i.e., each wire is given a specified homotopy class. If the wiring sketch is not given or the terminals are not fixed, the problem is NP-hard .

Some of the early work on continuous homotopic routing was done by Cole and Siegel  and Leiserson and Maley . They show that in the $`L_{\mathrm{}}`$ norm a solution can be found in $`O(k^3\mathrm{log}n)`$ time and $`O(k^3)`$ space, where $`n`$ is the number of wires and $`k`$ is the maximum of the input and output complexities of the wiring. Maley  shows how to extend the distance metric to arbitrary polygonal distance functions (including Euclidean distance) and presents a $`O(k^4\mathrm{log}n)`$ time and $`O(k^4)`$ space algorithm. The best result so far is due to Gao et al.  who present a $`O(kn^2\mathrm{log}(kn))`$ time and $`O(kn^2)`$ space algorithm. Duncan et al.  and Efrat et al.  present an $`O(kn+n^3)`$ algorithm for the related fat edge drawing problem: given a planar weighted graph $`G`$ with maximum degree 1 and an embedding for $`G`$, find a planar drawing such that all the edges are drawn as thick as possible and proportional to the corresponding edge weights.

The topological notion of homotopy formally captures the notion of deforming paths. Let $`\alpha ,\beta :[0,1]R^2`$ be two continuous curves parameterized by arc-length. Then $`\alpha `$ and $`\beta `$ are homotopic with respect to a set of obstacles $`VR^2`$ if $`\alpha `$ can be continuously deformed into $`\beta `$ while avoiding the obstacles; more formally, if there exists a continuous function $`h:[0,1]\times [0,1]R^2`$ with the following three properties:

1. $`h(0,t)=\alpha (t)`$ and $`h(1,t)=\beta (t)`$, for $`0t1`$
2. $`h(\lambda ,0)=\alpha (0)=\beta (0)`$ and $`h(\lambda ,1)=\alpha (1)=\beta (1)`$ for $`0\lambda 1`$
3. $`h(\lambda ,t)V`$ for $`0\lambda 1,0<t<1`$

Let $`\mathrm{\Pi }=\{\pi _1,\pi _2,\mathrm{},\pi _n\}`$ be a set of disjoint, simple polygonal paths and let the endpoints of the paths in $`\mathrm{\Pi }`$ define the set $`T`$ of at most $`2n`$ fixed points in the plane. Note that we allow a path to degenerate to a single point. We call the fixed points of $`T`$ “terminals,” and call the interior vertices of the paths “bends,” and use “points” in a more generic sense, e.g. “a point in the plane,” or “a point on a path.” We assume that no two terminals/bends lie on the same vertical line.

Our goal is to replace each path $`\pi _i\mathrm{\Pi }`$ by a shortest path $`\sigma _i`$ that is homotopic to $`\pi _i`$ with respect to the set of obstacles $`T`$, see Fig. 1 and Fig. 3. Note that $`\sigma _i`$ is unique. Let $`\mathrm{\Sigma }=\{\sigma _1,\mathrm{}\sigma _n\}`$ be the set of resulting paths. Observe that these output paths may \[self\] intersect by way of segments lying on top of each other, but will be non-crossing in the sense that a slight perturbation of the bends makes the paths simple and disjoint.

Let $`k_{\mathrm{𝑖𝑛}}`$ be the number of edges in all the paths of $`\mathrm{\Pi }`$. Let $`k_{\mathrm{𝑜𝑢𝑡}}`$ be the number of edges in all the paths of $`\mathrm{\Sigma }`$. We will measure the complexity of our algorithms in terms of $`n`$ and $`k=\mathrm{max}\{k_{\mathrm{𝑖𝑛}},k_{\mathrm{𝑜𝑢𝑡}}\}`$. In fact, the only relationship guaranteed among the parameters $`n,k_{\mathrm{𝑖𝑛}},k_{\mathrm{𝑜𝑢𝑡}}`$ is $`k_{\mathrm{𝑜𝑢𝑡}}nk_{\mathrm{𝑖𝑛}}`$. Note that $`k`$ may be arbitrarily large compared to $`n`$. Clearly this is the case since, for example, a path can wind around a set of terminals arbitrarily many times. However, there exist non-trivial cases in which $`k`$ can be much larger than $`n`$. Even after shortest paths have been computed for each wire, $`k`$ can be as large as $`k=\mathrm{\Omega }(2^n)`$, see Fig. 1. The algorithm of Hershberger and Snoeyink  runs in time $`O(nk)`$. The algorithm presented in this paper runs in time $`O(n^{2+\epsilon }+k\mathrm{log}^2n)`$, which is an improvement for $`n<k`$.

Although $`k`$ can be arbitrarily large compared to $`n`$, one easily forms the intuition that, because the paths are simple and disjoint, $`k`$ can be large in a non-trivial way only because path sections are repeated over and over. For example, a path may spiral arbitrarily many times around a set of points, but each wrap around the set is the same.

Our method makes essential use of this observation. We do not begin by explicitly searching for repeated path sections—this seems difficult modulo homotopic equivalence. Instead we begin in section 2 by applying vertical shortcuts to the paths (homotopically) so that each left and each right local extreme point occurs at a terminal. These terminals must then be part of the final shortest paths, and we have decomposed the paths into $`x`$-monotone pieces with endpoints at terminals. In section 3 we argue that the number of homotopically distinct $`x`$-monotone pieces is at most $`O(n)`$. We will bundle together all the homotopically equivalent pieces. Routing one representative from each such bundle using the straightforward “funnel” technique takes $`O(k+n^2)`$ time total. In section 3 we reduce this using a “shielding technique” where we again exploit the fact that the paths are disjoint and use the knowledge gained in routing one shortest path to avoid repeating work when we route subsequent paths. The final step of the algorithm is to unbundle, and recover the final paths by putting together the appropriate pieces. We summarize the main steps of the algorithm in Fig. 2.

Before we turn to the remainder of the paper, which consists of one section for each of steps 1, 2, and 3, it is worth noting that we rely on some powerful techniques. In step 1, we use simplex range search queries to perform homotopic simplifications on the paths. To do this, we use Chazelle’s cutting trees data structure . In step 2 we need to identify homotopically equivalent monotone path pieces. We use the efficient trapezoidization algorithm of Bar-Yehuda and Chazelle  to perform this step. Step 3 uses ideas from “funnel” algorithms for shortest paths . Step 4 is straight-forward.

We will find it convenient to regard each terminal as a small diamond. Otherwise when a vertical shortcut goes through a terminal we will need to specify whether the path actually goes to the left or the right of the terminal, which is awkward to say and to draw.

## 1 Shortcutting to Divide Paths into Monotone Pieces

We begin by applying vertical shortcuts to reduce each path to a sequence of $`x`$-monotone path sections, see Fig. 3. A vertical shortcut is a vertical line segment $`ab`$ joining a point $`a`$ on some path $`\pi `$ with a point $`b`$ also on $`\pi `$, and with the property that the subpath of $`\pi `$ joining $`a`$ and $`b`$, $`\pi _{ab}`$, is homotopic to the line segment $`ab`$.

We will only do elementary vertical shortcuts where the subpath $`\pi _{ab}`$ consists of \[portions of\] 2 line segments or 3 line segments with the middle one vertical, and the other two non-vertical. We distinguish left shortcuts which are elementary vertical shortcuts where $`\pi _{ab}`$ contains a point to the left of the line through $`ab`$; right shortcuts where $`\pi _{ab}`$ contains a point to the right of the line through $`ab`$; and collinear shortcuts where $`\pi _{ab}`$ lies in the line through $`ab`$, see Fig. 4 and Fig. 5. Collinear shortcuts are always applied after left/right shortcuts (and only then), and prevent consecutive vertical segments.

We will in fact only apply maximal elementary vertical shortcuts, where the subpath $`\pi _{ab}`$ cannot be increased. In particular, this means that for left and right shortcuts, either $`a`$ or $`b`$ is a bend, or the line segment $`ab`$ hits a terminal.

A local left \[right\] extreme of a path is a point or, more generally, a vertical segment, where the $`x`$-coordinate of the path reaches a local min \[max\]. Observe that every left \[right\] shortcut “cuts off” a local left \[right\] extreme. Conversely, every left \[right\] local extreme provides a left \[right\] shortcut unless the local extreme is locked at a terminal, meaning that the left \[right\] extreme contains the left \[right\] point of the terminal’s diamond. Figure 4 shows 4 cases of left extremes; in (a–c) they provide left shortcuts, but in case (d) the left extreme is locked at a terminal.

We will use range searching to detect, for a local left extreme, $`l`$, the maximal left shortcut that can be performed there. In particular, given $`l`$, we first identify the maximal potential shortcut in the absence of terminals. To do this, take the segment preceding $`l`$ and the segment following $`l`$ along the path. Take their right endpoints, and let $`t`$ be the leftmost of these two points. Then $`t`$ determines the maximal potential shortcut that can be performed at $`l`$. This potential shortcut forms either a triangle (for example, triangle $`abl`$ in Figure 4(a)) or a trapezoid (for example, the trapezoid determined by $`a`$, $`b`$, and $`l`$, in Figure 4(b)). If a range query tells us that the triangle/trapezoid is free of terminals, then it forms the true maximal shortcut from $`l`$. Otherwise, the range query should tell us the leftmost terminal inside the triangle/trapezoid, and that determines the maximal left shortcut from $`l`$. Note that the ranges we need to query are triangles with one vertical side, and trapezoids with two vertical sides.

Our algorithm is to perform elementary vertical shortcuts until none remain, at which time all local left and right extremes are locked at terminals, so the path divides into monotone pieces. See Section 1.2 for justification. Doing elementary vertical shortcuts in an arbitrary order may result in crossing paths as shown in Fig. 7; to guarantee non-crossing paths we will do the elementary vertical shortcuts in alternating left and right phases, where a left \[right\] phase means we do left \[right\] shortcuts (and the consequent collinear shortcuts) until no more are possible. See Section 1.3 for justification. We summarize these steps in Fig. 6.

The remainder of this section consists of four subsections: Subsection 1.1 justifies organizing the above algorithm around the sets of left and right extremes, rather than explicitly searching each time for a shortcut to perform. Subsection 1.2 establishes the fact that performing elementary vertical shortcuts enables us to divide the paths into $`x`$-monotone pieces with endpoints at terminals. Subsection 1.3 argues that doing elementary vertical shortcuts in left/right phases prevents crossings. Finally, Subsection 1.4 deals with implementation and run-time analysis.

### 1.1 Correctness

The point of this section is to show that the above algorithm correctly maintains the sets $`L`$ and $`R`$ of local left and right extremes not locked at terminals. Consider a left phase. The set $`R`$ is explicitly updated whenever it changes. For the correctness of $`L`$ we use:

###### Claim 1.1

Doing one left shortcut does not alter other local left extremes, nor the shortcuts that will be performed there.

Proof: A left shortcut only removes left portions of non-vertical segments.

We note that a similar claim fails for collinear shortcuts: doing one may prevent others.

### 1.2 Correctness of division into monotone pieces

###### Claim 1.2

At the end of a left \[right\] phase of shortcuts, each local left \[right\] extreme is locked at a terminal.

We note that more than one left and one right phase may be required, since the right phase may add new members to $`L`$, see Fig. 8.

###### Claim 1.3

Let $`\pi `$ be a path, and let $`\mu `$ be a result of performing left and right phases of elementary vertical shortcuts on $`\pi `$ until no more are possible. Suppose that the local left and right extremes of $`\mu `$ are locked at the terminals $`t_{i_1},\mathrm{},t_{i_l}`$ in that order. Let $`\sigma `$ be a shortest path homotopic to $`\pi `$. Then the local left \[right\] extremes of $`\sigma `$ are locked at exactly the same ordered list of terminals, and furthermore, the portion of $`\sigma `$ between $`t_{i_j}`$ and $`t_{i_{j+1}}`$ is a shortest path homotopic to the portion of $`\mu `$ between those same terminals.

Proof Sketch: Because $`\pi `$ and $`\mu `$ are homotopic, $`\sigma `$ is the shortest path homotopic to $`\mu `$. We can thus go from $`\mu `$ to $`\sigma `$ using “rubber band” deformations that only shorten the path, and such deformations cannot loosen a left \[right\] extreme from the terminal it is locked at.

### 1.3 Non-crossing paths

The purpose of this section is to prove the following:

###### Lemma 1.4

Each left \[right\] phase of shortcuts preserves the property that paths are non-crossing.

By symmetry, we can concentrate on a left phase. Note that a left phase, as we have described it, is non-deterministic. At each stage we choose one member $`l`$ from the set $`L`$ of current local left extremes, and use it to perform a left shortcut. To prove the Lemma we will first show that for each left phase there is some sequence of choices that preserves the property that the paths are non-crossing. We will then argue that the end result of a phase does not depend on the choices made.

###### Claim 1.5

Suppose we have a set of non-crossing paths. Let $`L`$ be the current set of local left extremes not locked at terminals, and let $`lL`$ be a rightmost element of $`L`$. Performing the shortcut for $`l`$ (together with any consequent collinear shortcuts) leaves the paths non-crossing.

Observe that we can—at least in theory—complete a left phase of shortcuts using this “rightmost” order. In practice we choose not to do this simply because of the extra time required to maintain a heap.

###### Lemma 1.6

The end result of a left phase does not depend on the sequence of choices made during the phase.

Proof: We begin with the claim that the set $``$ of all local left extremes that appear in $`L`$ over the course of the phase is independent of the choices made during the phase.

This implies that the set of left shortcuts performed during the course of the phase is also independent of the choices made during the phase. However, the set of collinear shortcuts is not independent. In particular, the order in which we perform left shortcuts affects the set of collinear shortcuts, see Fig. 9.

Consider one left phase. Let $`L_0`$ be the initial set of local left extremes not locked at terminals. Let $``$ be the union of $`L`$ over the course of the phase. Suppose two sequences of choices $`C_1`$ and $`C_2`$ during a left phase yield sets $`_1`$ and $`_2`$. We would like to show that $`_1=_2`$. Consider $`l_1`$. We will prove $`l_2`$ by induction on the number of left shortcuts performed in the phase before $`l`$ enters $`_1`$. If this number is 0 then $`lL_0`$ and we are done. Otherwise, $`l`$ is added to $`L`$ as a result of some left shortcut and consequent collinear shortcuts. Any vertical segment used in the collinear shortcuts may, in its turn, have arisen as a result of some left shortcut and consequent collinear shortcuts. Tracing this process, we find that $`l`$ is formed from a set of left shortcuts linked by vertical segments, all in the same vertical line as $`l`$, see Fig. 9. All these left shortcuts arose from local left extremes that entered $`_1`$ before $`l`$ did, and thus, by induction, are in $`_2`$. By Claim 1.1 left shortcuts are performed at each of these during the choice sequence $`C_2`$—though possibly in a different order than in $`C_1`$. The consequent collinear shortcuts will merge all the verticals forming $`l`$, and cannot merge more than that because the segments attached before and after $`l`$ are not vertical ($`l`$ is a local left extreme). Thus $`l`$ is in $`_2`$.

This proves that the set of local left extremes, and thus the set of left shortcuts is independent of choices made during the phase. Any vertical segment that is in the final set of paths output by the phase arises through left shortcuts plus consequent collinear shortcuts. Since any set of choices leads to the same set of left shortcuts, though possibly in different order, the consequent collinear shortcuts will arrive at the same final vertical shortcuts—i.e. the same final paths.

### 1.4 Implementation and run-time analysis

In order to perform range queries we need the cutting trees data structure of Chazelle . The cutting trees can be constructed in $`O(n^{2+\epsilon })`$ time and space, where $`\epsilon `$ is an arbitrarily small constant, and they support simplex range queries in time $`O(\mathrm{log}n)`$ . (Note that a trapezoid is a union of two triangles.) In case a triangle that we query is not empty, we need to find the rightmost/leftmost terminal inside it. To do this, create a segment tree $`𝒯`$ on the $`x`$-projections of the terminals, and maintain a cutting tree on each node of $`𝒯`$. This enables us to find the rightmost/leftmost terminal inside a query triangle in time $`O(\mathrm{log}^2n)`$, without increasing the asymptotic space and preprocessing time. Details are standard and omitted from this abstract.

###### Claim 1.7

The number of elementary vertical shortcuts that can be applied to a set of paths with a total of $`k`$ segments is at most $`2k`$.

Proof: Assume that no two terminals and/or bends line up vertically. Consider the set of vertical lines through bends and through the left and right sides of each terminal’s diamond. An elementary shortcut operates between two of these vertical lines. If a left \[right\] shortcut has its leftmost \[rightmost\] vertical at a bend, then after the shortcut, this vertical disappears forever, see Fig. 5. There are thus at most $`k`$ such elementary shortcuts. Consider, on the other hand, a left shortcut that has its leftmost vertical at a terminal. This only occurs when two previous left shortcuts are stopped at the terminal, and then combined in a collinear shortcut. See the right hand pictures of Fig. 5. Thus an original edge of the path has disappeared. Note that elementary shortcuts never fragment an edge of the path into two edges, but only shorten it from one end or the other. Thus there are at most $`k`$ elementary shortcuts of this type. Altogether, we obtain a bound of $`2k`$ elementary shortcuts.

This implies the following simple corollary:

###### Corollary 1.8

The running time of the shortcutting phase, not counting the preprocessing to construct cutting trees, is $`O(k\mathrm{log}^2n)`$.

Including preprocessing, step 1 takes time $`O(n^{2+\epsilon }+k\mathrm{log}^2n)`$.

## 2 Bundling Homotopically Identical Paths

Let $`M=\{\mu _1\mathrm{}\mu _m\}`$ be the set of $`x`$-monotone paths obtained from step 1 of our algorithm. In the second step of the algorithm we bundle homotopically equivalent paths in $`M`$. More precisely, we take one representative path for each equivalence class of homotopically equivalent paths in $`M`$. This is justified because the paths in each equivalence class have the same homotopic shortest path. Because the paths in $`M`$ are non-crossing and $`x`$-monotone, it is easier to detect homotopic equivalence: two paths are homotopically equivalent if they have the same endpoints, and, between these endpoints no terminal lies vertically above one path and vertically below the other.

In order to perform the bundling we use a trapezoidization of $`M`$, see Fig. 10. We apply the trapezoidization algorithm of Bar-Yehuda and Chazelle  to the paths obtained after the shortcuts of step 1, but before these paths are chopped into monotone pieces.

###### Claim 2.1

We can perturb the paths output from step 1 so that they become disjoint and the Bar-Yehuda and Chazelle algorithm can be applied.

The trapezoidization algorithm takes time $`O(k+n(\mathrm{log}n)^{1+\epsilon })`$, which is bounded by the time taken by step 1. Once we have a trapezoidization of $`M`$, we bundle homotopically equivalent paths as follows. While scanning each $`\mu M`$, we check if it is homotopically equivalent to the path “below” it, by examining all the edges of the trapezoidization that are incident to $`\mu `$ from below. If all these trapezoidization edges reach the same path $`\mu _j`$ and none pass through a terminal on the way to $`\mu _j`$, and $`\mu _i`$ and $`\mu _j`$ have the same terminals as endpoints, then we mark $`\mu _i`$ as a duplicate. Let $`R=\{\rho _1\mathrm{}\rho _r\}`$ be the paths of $`M`$ that are not marked as duplicates.

###### Lemma 2.2

The number of paths in $`R`$ is bounded by $`2n`$.

Proof: Note that every terminal $`t`$ is either a right endpoint of paths in $`R`$ or a left endpoint of paths in $`R`$, but not both. This is because a terminal cannot have local left extremes and local right extremes locked at it without the paths crossing.

For each terminal $`tT`$, associate its bottommost incident path $`\varphi (t)`$ (its “floor”), and the first path hit by a vertical ray going up from $`t`$, $`\gamma (t)`$ (its “roof”). We claim that every path in $`R`$ is $`\varphi (t)`$ or $`\gamma (t)`$ for some $`tT`$. This proves that the number of paths in $`R`$ is at most $`2n`$.

Consider a path $`\rho R`$ with left and right terminals $`s`$ and $`t`$, respectively. Suppose that $`\rho `$ is not a roof. Then every vertical ray extending downward from a point of $`\rho `$ must hit the same path $`\sigma `$. (If two rays hit different paths, then in the middle some ray must hit a terminal.) Furthermore, since $`\rho `$ is not the bottommost path incident to $`s`$, $`\sigma `$ must hit $`s`$. Similarly $`\sigma `$ must hit $`t`$. But then $`\sigma `$ and $`\rho `$ are homotopically equivalent.

## 3 Shortest Paths

In this section we find shortest paths homotopic to the $`O(n)`$ monotone paths $`R=\{\rho _1,\mathrm{},\rho _r\}`$ produced in the previous section. Let $`\rho _i^{}`$ denote the shortest path homotopic to $`\rho _i`$. We route each path using a funnel technique. The funnel algorithm of  and  operates on a triangulation of the $`n`$ points (the terminals in our case), and follows the path through the triangulation maintaining a current “funnel” containing all possible shortest paths to this point. The algorithm takes time proportional to the number of edges in the path plus the number of intersections between the triangulation edges and the path.

Rather than a triangulation, we will use a trapezoidization formed by passing a vertical line through each of the $`n`$ terminals, see Fig. 11(a). Then, since each path $`\rho _i`$ is $`x`$-monotone, it has $`O(n)`$ intersections with trapezoid edges, and the funnel algorithm takes time $`O(n+k_i)`$, where $`k_i`$ is the number of edges in $`\rho _i`$. This gives a total over all paths of $`O(n^2+k)`$. Recall that $`k`$ is $`\mathrm{max}\{k_{\mathrm{𝑖𝑛}},k_{\mathrm{𝑜𝑢𝑡}}\}`$, where $`k_{\mathrm{𝑖𝑛}}`$ is the number of edges in all the input paths, and $`k_{\mathrm{𝑜𝑢𝑡}}`$ is the number of edges in all the output paths. Note that the number of edges in $`R`$ is bounded by $`k_{\mathrm{𝑖𝑛}}`$.

With respect to our whole algorithm, this time of $`O(n^2+k)`$ is dominated by the $`O(n^{2+\epsilon })`$ time required to create the range query data structure in step 1. However, it is interesting to see how much improvement we can make in this step alone. In the remainder of this section we describe a randomized algorithm to route the $`O(n)`$ monotone paths of $`R`$ in time $`O(n\mathrm{log}n+k)`$, and finally mention a deterministic algorithm with $`O(n\sqrt{n}+k)`$ running time.

Both methods use a “shielding” technique. We begin by describing this idea intuitively. First note that the $`\rho _i`$’s can be routed independently, since none affects the others. The initial paths $`\rho _i`$ are non-crossing, and so are the final shortest paths, $`\rho _i^{}`$. If $`\rho _j`$ is below $`\rho _i`$, and we have already computed $`\rho _j^{}`$, then $`\rho _j^{}`$ behaves as a barrier that “shields” $`\rho _i`$ from terminals that are vertically below $`\rho _j^{}`$, see Fig. 11(b).

To utilize shielding we will modify the basic trapezoidization described above as we discover shortest paths $`\rho _i^{}`$. In particular, the upward vertical ray through terminal $`t`$, $`u(t)`$, will be truncated at the lowest shortest path $`\rho _i^{}`$ that is strictly above $`t`$ and does not bend at $`t`$. The downward vertical ray through terminal $`t`$, $`d(t)`$, will be truncated in an analogous way.

The shortest paths found so far, together with the truncated vertical rays $`u(t)`$ and $`d(t)`$ for each terminal $`t`$, partition the plane into trapezoids, each bounded from left and right by the vertical rays, and from above and below by shortest paths. To route a new path $`\rho _i`$ through this modified trapezoidization we use the following algorithm.

Routing $`\rho _i`$ with shielding

1. Identify the first trapezoid that $`\rho _i^{}`$ will traverse. This can be done in $`O(\mathrm{log}n)`$ time because the shortest paths observe the same vertical ordering as the original $`\rho _j`$’s.
2. Traverse from left to right the sequence of trapezoids that $`\rho _i^{}`$ will pass through. (Note that $`\rho _i`$ itself may pass through different trapezoids, see Fig. 11(b).) We construct the funnel for $`\rho _i^{}`$ as we do this traversal. Suppose that our path enters trapezoid $`\tau `$. To leave $`\tau `$ on the right we have two cases. If the right side of $`\tau `$ is a point, then it is a terminal $`t`$, and we are locked between two paths that terminate or bend at $`t`$. Then the funnel collapses to this point, and we proceed to the next trapezoid if $`\rho _i`$ continues. Otherwise the right side of $`\tau `$ is a vertical through some terminal $`t`$, and (unless $`\rho _i`$ ends at $`t`$) we have a choice of two trapezoids to enter, the upper one with left side $`u(t)`$ or the lower one with left side $`d(t)`$. We follow path $`\rho _i`$ until it crosses the infinite vertical line through $`t`$. If it passes above $`t`$ then we enter the upper trapezoid, and otherwise we enter the lower trapezoid. We update the funnel to include $`t`$.
3. When we reach the right endpoint of $`\rho _i`$, the funnel gives us the shortest homotopic path $`\rho _i^{}`$.
4. We update the trapezoidization as follows. For each vertical segment $`u(t)`$ or $`d(t)`$ that is intersected by $`\rho _i^{}`$ we chop the segment at its intersection point with $`\rho _i^{}`$, provided that the intersection point is not $`t`$ (i.e. that $`\rho _i^{}`$ does not bend at or terminate at $`t`$).

Without yet discussing the order in which we route the paths, we can say a bit about the timing of this shielding method. As we traverse a path $`\rho _i`$ we spend time proportional to the size of $`\rho _i`$ plus the number of trapezoids traversed by $`\rho _i^{}`$. When $`\rho _i^{}`$ leaves a trapezoid, say at the line segment $`u(t)`$ above terminal $`t`$, it may happen that $`\rho _i^{}`$ will bend at $`t`$ or terminate at $`t`$. In this case $`t`$ is part of the output path $`\rho _i^{}`$, and we can charge the work for this trapezoid to the output size. If, on the other hand, $`\rho _i^{}`$ does not bend or terminate at $`t`$, then it crosses $`u(t)`$ and we chop $`u(t)`$ there. In this case we charge the work for this trapezoid to the chop. Thus the total time spent by the algorithm is $`O(k+C)`$ where $`C`$ is the total number of chops performed at the $`n`$ verticals.

For shielding to be effective we need to route paths in an order that makes $`C`$ grow more slowly than $`n^2`$. We first analyze the randomized algorithm where we choose the next path to route at random with uniform probability from the remaining paths.

###### Claim 3.1

If paths are routed in random order then the expected number of times we chop a segment $`u(t)`$ or $`d(t)`$ is $`O(\mathrm{log}r)`$, and thus $`C`$ is $`O(n\mathrm{log}n)`$, and the routing takes time $`O(k+n\mathrm{log}n)`$.

Proof: This follows from a standard backward analysis. See for example  for other proofs along these lines. Let $`l`$ be a vertical line through a point $`t`$, and let $`u`$ denote a ray emerging vertically from $`t`$. Assume that $`m`$ paths intersecting $`u`$ have been inserted up to now, and we are about to insert a new one. Then $`u`$ will be chopped if and only if the new path creates an intersection point below all existing intersection points on $`u`$, but above $`t`$ itself. Since the order of the insertions is random, the probability of the new intersection point being below all other intersection points is $`1/m`$. Summing over all paths yields the claimed bound.

Finally, we mention that we can achieve a routing time of $`O(k+n\sqrt{n})`$ deterministically. Suppose that the paths $`\{\rho _1,\mathrm{},\rho _r\}`$ are in order from top to bottom. (We get this for free from step 2 of our algorithm.) Partition the paths into blocks $`B_1,\mathrm{},B_\sqrt{r}`$ each of size $`\sqrt{r}`$, and route them one block at a time from $`B_\sqrt{r}`$ to $`B_1`$. Within each block we process the paths in order. Since the largest increasing \[decreasing\] sequence with this ordering has size $`\sqrt{r}`$, therefore the number of chops at each vertical is $`\sqrt{r}`$. Thus $`C`$ is $`O(n\sqrt{n})`$, and this routing method takes time $`O(k+n\sqrt{n})`$.

## 4 Conclusion and Open Problems

For any set of $`n`$ disjoint paths joining pairs of terminals in the plane we can find shortest paths homotopic with respect to the set of terminals in time $`O(n^{2+\epsilon }+k\mathrm{log}^2n)`$ where $`k`$ is the sum of input and output sizes of the paths. If $`k`$ is larger than $`n`$ this is better than the Hershberger-Snoeyink algorithm which runs in time $`O(nk)`$.

More generally, we can use any other range search method, and obtain a running time of $`O(P+kQ+n(\mathrm{log}n)^{1+\epsilon }+T)`$ where $`P`$ is the preprocessing time for the range search data structure on $`n`$ points, $`Q`$ is the time for a simplex optimization query (find the minimum/maximum $`x`$-coordinate point in a triangle), and $`T`$ is the time for the shortest path method of step 3—$`O(n\mathrm{log}n)`$ randomized, or $`O(n\sqrt{n})`$ deterministic. Whether $`k`$ is large or small compared to $`n`$ determines which trade-off between $`P`$ and $`Q`$ is preferred. For example, the partition tree method of Matoušek  yields an emptiness query time of $`O(\sqrt{n})`$ with near-linear time preprocessing, and would be preferable for $`k=O(n^{3/2})`$. With any approach to range searching, it may be possible to improve the query time by taking into account the fact that one of the edges of each query triangle is vertical.

Another open problem is whether we can do without range queries, and somehow use a trapezoidization of the original paths, since this can be found so efficiently with the Bar-Yehuda and Chazelle algorithm.


text:
# Thresholds to Chaos and Ionization for the Hydrogen Atom in Rotating Fields

## I Introduction

The chaotic ionization of the hydrogen atom in a variety of external fields is a fundamental problem in nonlinear dynamics and atomic physics. Early work, in particular, focused on the ionization mechanism in a linearly polarized (LP) microwave field casa87 ; koch95 . This problem is noteworthy because it showed the general applicability of the “Chirikov Resonance Overlap Criterion” chir79  to real quantum mechanical systems. This empirical criterion predicts the value at which two nearby resonances overlap and large scale stochasticity occurs. There is no doubt that the Chirikov criterion is generally robust and, therefore, it is no surprise that Chirikov’s criteria has been tried as a way to predict transitions to chaos and ionization dynamics in more complicated circumstances, e.g., for the hydrogen atom in circularly polarized (CP) microwave fields. However, the success of the Chirikov criterion in describing the LP problem stands in contrast with what seems to have been a somewhat mixed performance when applied to the hydrogen atom in rotating microwave fields. Partly for this reason, a good deal of controversy has surrounded the ionization mechanism and also the interpretation of the Chirikov criterion when applied to ionization in rotating fields. In this article we propose to examine the ionization of the hydrogen atom in a circularly polarized microwave field (CP). Our analysis will use an advanced method, the renormalization map koch99 ; chan01R , so as to include higher order resonances beyond the Chirikov approach.

Interest in the CP microwave problem began with experiments by the Gallagher group fu90  which revealed a strong dependence of the ionization threshold on the degree of polarization. They explained their CP results proposing that in a rotating frame, ionization proceeds in roughly the same way as for a static field, i.e., a static field has the same effect whether its coordinate system is rotating or not. In a Comment, Nauenberg naue90  argued that the ionization mechanism was substantially more complicated and that the effect of rotation on the ionization threshold had to be taken into account. Fully three-dimensional numerical simulations by Griffiths and Farrelly grif92  were able to provide quantitative agreement with the experimental results. Griffiths and Farrelly grif92  and Wintgen wint91  independently proposed similar models for ionization based on the Runge-Lenz vector.

Subsequent theoretical work can be divided into three main categories: $`(i)`$ classical simulations kapp93 , $`(ii)`$ quantum simulations, and $`(iii)`$ resonance overlap studies. The two main papers on resonance overlap are by Howard howa92  who published the first such study of this system, and by Sacha and Zakrzewski sach97 . In both cases, the Hamiltonian was written in an appropriate rotating frame, a choice of ‘zero-order’ actions was made and the Chirikov machinery invoked. Somewhat surprisingly, the paper by Sacha and Zakrzewski sach97  disagrees in a number of key areas with the results of Howard howa92  as well as with some conclusions drawn from numerical studies farr95 . In a paper by Brunello et al. brun97 , it was shown both numerically and analytically that the actual ionization threshold observed in an experiment must take into account the manner in which the field was turned on; essentially, in some experiments electrons are switched during the field “turn on” directly into unbound regions of phase space. That is, the underlying resonance structure of the Hamiltonian is almost irrelevant because ionization is accomplished by the ramp-up of the field. For this reason, the application of resonance overlap criteria in rotating frames can be quite intricate and this provides some explanation for the apparent inadequacy of the Chirikov criterion - if ionization occurs during the field turn-on time then, of course, resonance overlap is irrelevant.

In this article, we study the hydrogen atom driven by a CP microwave field together with a magnetic field lying perpendicular to the polarization plane (CP$`\times B`$). The magnetic field has been introduced to prevent ionization in the plane. This provides the opportunity to eliminate difficulties associated with the turn on of the field, thus opening up the way to a study of resonance overlap in a more controlled manner in the rotating frame. As noted, without the added magnetic field all the electrons may have gone before the resonances have had a chance to overlap.

The paper is organized as follows: Section II introduces the Hamiltonian and provides a discussion of resonance overlap in the rotating frame in both the Chirikov and renormalization approaches. In order to obtain qualitative and quantitative results about the dynamics of the hydrogen atom in CP$`\times `$B fields, we compare numerically, in Sec. III, Chirikov’s resonance overlap criterion with our renormalization group transformation method. The Chirikov method is found to provide good qualitative results and is useful because of its simplicity. The renormalization transformation is used first to check the qualitative features obtained by the resonance overlap criterion, and to obtain more quantitative results about the dynamics by expanding the Chirikov approach. Conclusions are in Sec. IV.

## II Hydrogen atom in CP$`\times `$B fields

We consider a hydrogen atom perturbed by a microwave field of amplitude $`F`$ and frequency $`\omega _f`$, circularly polarized in the orbital plane, and a magnetic field with cyclotron frequency $`\omega _c`$. The Hamiltonian in atomic units and in the rotating frame of the microwave field is Bfrie91  :

$$H(p_x,p_y,x,y)=\frac{p_x^2+p_y^2}{2}\frac{1}{\sqrt{x^2+y^2}}\mathrm{\Omega }(xp_yyp_x)+Fx+\frac{\omega _c^2}{8}(x^2+y^2),$$
(1)
where $`\mathrm{\Omega }=\omega _f\omega _c/2`$. Following Ref. howa92 , we rewrite Hamiltonian (1) in the action-angle variables $`(J,L,\theta ,\psi )`$ of the problem with $`\omega _c=0`$. The angles $`\theta `$ and $`\psi `$ are conjugate respectively to the principal action $`J`$ and to the angular momentum $`L`$. The Hamiltonian becomes :

$$H(J,L,\theta ,\psi )=H_0(J,L,\theta )+FJ^2\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}V_n(e)\mathrm{cos}(n\theta +\psi ),$$
(2)
where the integrable Hamiltonian $`H_0`$ is

$$H_0(J,L,\theta )=\frac{1}{2J^2}\mathrm{\Omega }L+J^4\frac{\omega _c^2}{8}\underset{n,m=\mathrm{}}{\overset{+\mathrm{}}{}}V_mV_{mn}\mathrm{cos}n\theta .$$
(3)
The coefficients $`V_n`$ of the expansion are the following ones :

$`V_0(e)={\displaystyle \frac{3e}{2}},`$ (4)
$`V_n(e)={\displaystyle \frac{1}{n}}\left[𝒥_n^{}(ne)+{\displaystyle \frac{\sqrt{1e^2}}{e}}𝒥_n(ne)\right],\text{ for }n0,`$ (5)
where $`𝒥_n`$ is the $`n`$th Bessel function of the first kind and $`𝒥_n^{}`$ its derivative. The parameter $`e`$ is given by

$$e=\left(1\frac{L^2}{J^2}\right)^{1/2}.$$

The positions $`x`$ and $`y`$ are given by the following formulas :

$`x={\displaystyle \underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}}V_n(e)\mathrm{cos}(n\theta +\psi ),`$
$`y={\displaystyle \underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}}V_n(e)\mathrm{sin}(n\theta +\psi ).`$
The Hamiltonian (2) can be rescaled in order to withdraw the dependence on $`\mathrm{\Omega }`$. We rescale time by a factor $`\mathrm{\Omega }`$ \[ we divide Hamiltonian (2) by $`\mathrm{\Omega }`$\]. We rescale the actions $`J`$ and $`L`$ by a factor $`\lambda =\mathrm{\Omega }^{1/3}`$, i.e. we replace $`H(J,L,\theta ,\psi )`$ by $`\lambda H(J/\lambda ,L/\lambda ,\theta ,\psi )`$. We notice that this rescaling does not modify the eccentricity $`e`$. The resulting Hamiltonian becomes :

$$H=H_0+F^{}J^2\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}V_n\mathrm{cos}(n\theta +\psi ),$$
(6)
where

$$H_0=\frac{1}{2J^2}L+J^4\frac{\omega _c^2}{8}\underset{n,m=\mathrm{}}{\overset{+\mathrm{}}{}}V_mV_{mn}\mathrm{cos}n\theta .$$
(7)
and $`F^{}=F\mathrm{\Omega }^{4/3}`$ is the rescaled amplitude of the field and $`\omega _c^{}=\omega _c/\mathrm{\Omega }`$. In what follows we assume that $`\mathrm{\Omega }=1`$. Furthermore, for simplicity we assume that $`e`$ is a parameter of the system equal to the initial eccentricity of the orbit in the Keplerian problem ($`\omega _c=0`$ and $`F=0`$).
In this paper, we consider the high-scaled frequency regime for co-rotating orbits : $`\mathrm{\Omega }/\omega _K>1`$ where the Kepler frequency $`\omega _K`$ is equal to $`1/J^3`$, i.e., we consider part of phase space corresponding to $`0<J<1`$.

### II.1 Study of the integrable part of Hamiltonian (2)

Several key features of the dynamics of Hamiltonian (2) can be understood by looking at the integrable part. We consider the mean value with respect to $`\theta `$ of the integrable Hamiltonian (7) :

$$\stackrel{~}{H}_0=\frac{1}{2J^2}L+J^4\frac{\omega _c^2}{8}V^2,$$
(8)
where $`V^2=_{n=\mathrm{}}^+\mathrm{}V_n^2`$, or another way for considering $`\stackrel{~}{H_0}`$ is to assume that $`J^4\frac{\omega _c^2}{8}_{n0,m}V_mV_{mn}\mathrm{cos}n\theta `$ is part of the perturbation of Hamiltonian (6). One of the main features of Hamiltonian $`\stackrel{~}{H_0}`$ is that it does not satisfy the standard twist condition for $`\omega _c0`$. Since the Hessian of this Hamiltonian is

$$\frac{^2\stackrel{~}{H_0}}{J^2}=\frac{3}{J^4}\left(\frac{1}{2}J^6\omega _c^2V^21\right),$$

the phase space is divided into two main parts : a positive twist region where $`^2\stackrel{~}{H_0}/J^2>0`$ for $`J(\sqrt{2}/\omega _cV)^{1/3}`$, and a negative one where $`^2\stackrel{~}{H_0}/J^2<0`$ for $`J(\sqrt{2}/\omega _cV)^{1/3}`$. Both regions are separated by a twistless region ($`^2\stackrel{~}{H_0}/J^2=0`$). We notice that the singularity at $`J=0`$ is located inside the negative twist region, that the energy in the negative twist region is negative, and that the positive twist and twistless regions do not exist in the absence of magnetic field. Furthermore, the size of the negative twist region decreases like $`\omega _c^{1/3}`$ as one increases the magnetic field $`\omega _c`$, i.e., it shrinks to the singularity of the Hamiltonian ($`J=0`$).

The phase space of $`\stackrel{~}{H_0}`$, as well as the one of Hamiltonian (7), is foliated by invariant tori, the main difference being that the invariant tori for $`\stackrel{~}{H_0}`$ are flat in these coordinates. We consider a motion with frequency $`\omega `$, i.e. such that the dynamics is $`\theta (t)=\omega t+\theta (0)mod2\pi `$ and $`J(t)=J(0)`$. The associated invariant torus is located at $`J_\omega `$ such that the frequency $`\frac{\stackrel{~}{H}_0}{J}`$ at $`J=J_\omega `$ is equal to $`\omega `$. The equation determining $`J_\omega `$ is :

$$\frac{\omega _c^2}{2}V^2J_\omega ^6\omega J_\omega ^3+1=0.$$
(9)
There are two real positive solutions of this equation :

$$J_\omega ^\pm =\left(\frac{\omega \pm \sqrt{\omega ^22\omega _c^2V^2}}{\omega _c^2V^2}\right)^{1/3}.$$

The condition of existence of an invariant torus with frequency $`\omega `$ for $`\stackrel{~}{H_0}`$ is $`\omega \sqrt{2}\omega _cV`$. There are two invariant tori with frequency $`\omega `$: one located at $`J_\omega ^{}`$ in the negative twist region is a continuation of the torus with frequency $`\omega `$ in the absence of magnetic field since $`lim_{\omega _c0}J_\omega ^{}=\omega ^{1/3}`$; the other torus located at $`J_\omega ^+`$ in the positive twist region is created far from the singularity $`J=0`$ as soon as the field is non-zero. Figure 1 depicts the position of these tori. We notice that as soon as $`\omega _c0`$ there is creation of a set of invariant tori far from $`J=0`$ in the positive twist region.
If we increase $`\omega _c`$, the position of the positive twist torus decreases whereas the position of the negative twist torus increases. Both tori collide at $`\omega _c=\omega /\sqrt{2}V`$ to a twistless invariant torus (of the same frequency $`\omega `$) located at $`(\omega /2)^{1/3}`$. The approximate location of this twistless torus as a function of $`\omega _c`$ (for fixed parameter $`e`$) is plotted in Fig. 1. As one increases $`\omega _c`$, a large portion of the invariant tori with $`\omega [0,1]`$ disappears.

Remark : First Order Delaunay normalization lanc97
Averaging Hamiltonian (6) over $`\theta `$ gives

$$H_\theta =\frac{1}{2J^2}L+J^4\frac{\omega _c^2}{8}V^2\frac{3e}{2}J^2F\mathrm{cos}\psi .$$

Using the expansion of $`V^2`$ to the second order of $`e`$, $`V^21+\frac{3e^2}{2}`$, and the fact that the previous Hamiltonian does not depend on $`\theta `$ (thus $`J`$ is constant), the Hamiltonian reduces to

$$𝒦=L+\frac{3\omega _c^2}{16}e^2J^4\frac{3}{2}eJ^2F\mathrm{cos}\psi ,$$

which is the Hamiltonian studied in Ref. lanc97  to find bifurcation of equilibrium points.

### II.2 Primary main resonances and Chirikov resonance overlap

The positive and negative twist regions have their own sets of primary resonances given by Hamiltonian (6). The approximate locations of these primary resonances $`m`$:1, denoted $`J_m^\pm `$, are obtained by the condition $`m\dot{\theta }+\dot{\psi }0`$. There are two such resonances located approximately at

$$J_m^\pm =\left(\frac{1\pm \sqrt{12m^2\omega _c^2V^2}}{m\omega _c^2V^2}\right)^{1/3}.$$
(10)
The resonance located at $`J_m^{}`$ is the continuation of the resonance in the case $`\omega _c=0`$. The condition of existence of these resonances is $`m\sqrt{2}\omega _cV<1`$. Similarly to collisions of invariant tori, collisions of periodic orbits occur when increasing $`\omega _c`$. Figure 2 depicts the different domains of existence of real $`J_m`$ for first primary resonances $`(m=1,\mathrm{},5)`$ in the plane of parameters $`e`$-$`\omega _c`$. We notice that the most relevant parameter in this problem is the magnetic field. The variations of the dynamics induced by the eccentricity $`e`$ are smooth.

Between two consecutive resonances $`m`$:1 and $`m`$+1:1 (if they exist), regular and chaotic motions occur for a typical value of the field $`F`$. In order to estimate the chaos threshold between these resonances, i.e. the value for which there is no longer any rotational invariant torus acting like a barrier in phase space, Chirikov overlap criterion provides an upper bound but usually quite far from realistic values (obtained by numerical integration for instance). For quantitatively accurate thresholds, a modified criterion is applied : the 2/3-rule criterion.

In order to compute the resonance overlap value of the field $`F`$ between $`m`$:1 and $`m`$+1:1 primary resonances, we follow the procedure described in Ref. meer82 . First, we change the frame to a rotating one at the phase velocity of resonance $`m`$:1. We apply the following canonical change of variables : $`(𝑨^{},𝝋^{})=(\stackrel{~}{N}^1𝑨,N𝝋)`$ where $`𝑨=(J,L)`$ and $`𝝋=(\theta ,\psi )`$ and

$$N=\left(\begin{array}{cc}m& 1\\ 0& 1\end{array}\right),$$

and $`\stackrel{~}{N}`$ denotes the transposed matrix of $`N`$. Hamiltonian (6) is mapped into

$`H=`$ $`{\displaystyle \frac{1}{2m^2J^2}}J^{}L^{}+{\displaystyle \frac{\omega _c^2}{8}}m^4J^4{\displaystyle \underset{n,n^{}}{}}V_n^{}V_{n^{}n}\mathrm{cos}[n(\theta ^{}\psi ^{})/m]`$
$`+Fm^2J^2{\displaystyle \underset{n}{}}V_n\mathrm{cos}[(n\theta ^{}+(mn)\psi ^{})/m].`$
By averaging over the fast variable $`\psi ^{}`$, the Hamiltonian becomes

$$H=\frac{1}{2m^2J^2}J^{}+\frac{\omega _c^2}{8}m^4J^4V^2+Fm^2V_mJ^2\mathrm{cos}\theta ^{}.$$

The resulting Hamiltonian is integrable and one can compute the width of the resonance $`m`$:1 following Ref. meer82 . We expand the previous Hamiltonian around $`J_m^{}=J_m/m`$ and keep only the quadratic part in $`\mathrm{\Delta }J^{}=J^{}J_m^{}`$ and the constant term in the action $`\mathrm{\Delta }J^{}`$ proportional to $`\mathrm{cos}\theta ^{}`$ :

$$H=\frac{3m^2}{2J_m^4}\left(1\frac{1}{2}\omega _c^2V^2J_m^6\right)\mathrm{\Delta }J^2+FJ_m^2V_m\mathrm{cos}\theta ^{}.$$

We rescale energy by a factor $`3m^2J_m^4(1\omega _c^2V^2J_m^6/2)`$ :

$$H=\frac{1}{2}\mathrm{\Delta }J^2\frac{FJ_m^6V_m}{3m^2\left(1\frac{1}{2}\omega _c^2V^2J_m^6\right)}\mathrm{cos}\theta ^{}.$$

Depending on the positive or negative twist region, this rescaling is positive or negative respectively.

The width of the resonance $`m`$:1 in the variables $`\mathrm{\Delta }J=m\mathrm{\Delta }J^{}`$ is given by :

$$\mathrm{\Delta }_m=4J_m^3\sqrt{\frac{FV_m(e)}{3}}\left|2\frac{J_m^3}{m}\right|^{1/2},$$

since $`\omega _c^2V^2J_m^6/2=1+J_m^3/m`$. The 2/3-rule criterion for the critical threshold for the overlap between resonance $`m`$:1 and $`m`$+1:1 is reached when the distance between two neighboring primary resonances is equal (up to a factor $`2/3`$) to the sum of the half-widths of these resonances :

$$\frac{2}{3}|J_{m+1}J_m|=\frac{1}{2}(\mathrm{\Delta }_m+\mathrm{\Delta }_{m+1}).$$

The critical threshold is given by :

$$F_m(e,\omega _c)=\frac{(J_{m+1}J_m)^2}{3\left(J_m^3\sqrt{V_m}\left|2\frac{J_m^3}{m}\right|^{1/2}+J_{m+1}^3\sqrt{V_{m+1}}\left|2\frac{J_{m+1}^3}{m+1}\right|^{1/2}\right)^2},$$
(11)
where $`J_m`$ stands for either $`J_m^+`$ or $`J_m^{}`$ given by Eq. (10). Therefore we obtain two critical couplings : one in the positive twist region, $`F_m^+(e,\omega _c)`$, and one in the negative twist region, $`F_m^{}(e,\omega _c)`$. We notice that for $`\omega _c=0`$, since $`J_m^{}=m^{1/3}`$, we obtain the formula of Ref. sach97  for the chaos threshold in the CP problem. For $`\omega _c=1/[\sqrt{2}V(m+1)]`$, the threshold $`F_m(e,\omega _c)`$ vanishes. This case corresponds to the collision of the resonances $`m`$+1:1 (the positive and negative twist ones). Therefore, the Chirikov criterion is valid only for $`\omega 1/[\sqrt{2}V(m+1)]`$, i.e. when the two primary resonances $`m`$:1 and $`m`$+1:1 exist in the system.

We use this criterion in order to study the stability in the different regions of phase space as a function of the magnetic field $`\omega _c`$ (for small values of the field $`\omega _c`$) and the eccentricity of the initial orbit $`e`$. However, since this criterion is purely empirical, we use another method to validate or refine the results : we use the renormalization method which has proved to be a very powerful and accurate method for determining chaos thresholds in this type of models chan01R ; chan01a ; chan02a ; chan00d . We compare the results given by both methods in the region where the criterion applies and we use the renormalization map to compute chaos thresholds when Eq. (11) does not apply (when one of the main primary resonance $`m`$:1 has disappeared).

### II.3 Renormalization method

The Chirikov resonance overlap criterion gives us a value for the chaos threshold between two neighboring primary resonances. This value aims at approximating the value of the field $`F`$ for which there is no longer any barrier in phase space, and for which large-scale diffusion of trajectories occur between these resonances. The renormalization method gives a more local and more accurate information. This method computes the threshold of break-up of an invariant torus with a given frequency $`\omega `$. Then by varying $`\omega `$, we obtain the global information on the transition to chaos.

#### II.3.1 Expansion of the Hamiltonian around a torus

In order to apply the renormalization transformation as defined in Refs. chan99c ; chan00a  for a given torus with frequency $`\omega `$, we expand Hamiltonian (2) in Taylor series in the action $`J`$ around $`J=J_\omega `$.

$$\frac{1}{2J^2}=\underset{k=0}{\overset{+\mathrm{}}{}}(1)^{k+1}\frac{k+1}{2J_\omega ^{k+2}}\mathrm{\Delta }J^k,$$

where $`\mathrm{\Delta }J=JJ_\omega `$. The meanvalue of the quadratic term in $`\mathrm{\Delta }J`$ of $`H_0`$ is equal to $`\frac{3(\omega J_\omega ^32)}{2J_\omega ^4}\mathrm{\Delta }J^2`$. We rescale the action variables such that this quadratic term is equal to 1/2, i.e. we replace $`H(\mathrm{\Delta }J,L,\theta ,\psi )`$ by $`\lambda H(\mathrm{\Delta }J/\lambda ,L/\lambda ,\theta ,\psi )`$ with $`\lambda =3(\omega J_\omega ^32)/J_\omega ^4`$. We notice that this rescaling can be done except for two cases $`J_\omega =(\omega /2)^{1/3}`$ which is the twistless case, and for $`J_\omega =0`$ which is the singularity. Therefore, as it is defined in Refs. chan01R ; chan99c ; chan00a , the renormalization cannot be applied in the twistless region.

The Hamiltonian (2) becomes :

$$H_F=\omega \mathrm{\Delta }JL+\underset{k2}{\overset{+\mathrm{}}{}}H_{k,0}\mathrm{\Delta }J^k+\underset{k=0}{\overset{4}{}}\underset{n>0}{}H_{k,n}\mathrm{\Delta }J^k\mathrm{cos}n\theta +F\underset{k=0}{\overset{2}{}}\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}V_{k,n}\mathrm{\Delta }J^k\mathrm{cos}(n\theta +\psi ),$$
(12)
where

$`H_{2,0}={\displaystyle \frac{1}{2}},`$
$`H_{3,0}={\displaystyle \frac{J_\omega ^3(\omega J_\omega ^3+1)}{9(2\omega J_\omega ^3)^2}},`$
$`H_{4,0}={\displaystyle \frac{J_\omega ^6(11\omega J_\omega ^3)}{108(2\omega J_\omega ^3)^3}},`$
$`H_{k,0}={\displaystyle \frac{(k+1)J_\omega ^{3k6}}{23^{k1}(2\omega J_\omega ^3)^{k1}}}\text{ for }k5,`$
$`H_{0,n}={\displaystyle \frac{3}{4}}\omega _c^2(2\omega J_\omega ^3){\displaystyle \underset{m=\mathrm{}}{\overset{+\mathrm{}}{}}}V_mV_{mn},`$
$`H_{1,n}=\omega _c^2J_\omega ^3{\displaystyle \underset{m=\mathrm{}}{\overset{+\mathrm{}}{}}}V_mV_{mn},`$
$`H_{2,n}={\displaystyle \frac{\omega J_\omega ^31}{2\omega J_\omega ^3}}{\displaystyle \frac{_{m=\mathrm{}}^+\mathrm{}V_mV_{mn}}{V^2}},`$
$`H_{3,n}={\displaystyle \frac{2J_\omega ^3(\omega J_\omega ^31)}{9(2\omega J_\omega ^3)^2}}{\displaystyle \frac{_{m=\mathrm{}}^+\mathrm{}V_mV_{mn}}{V^2}},`$
$`H_{4,n}={\displaystyle \frac{J_\omega ^6(\omega J_\omega ^31)}{54(2\omega J_\omega ^3)^3}}{\displaystyle \frac{_{m=\mathrm{}}^+\mathrm{}V_mV_{mn}}{V^2}},`$
$`V_{0,n}=3V_n{\displaystyle \frac{2\omega J_\omega ^3}{J_\omega ^2}},`$
$`V_{1,n}=2J_\omega V_n,`$
$`V_{2,n}=V_n{\displaystyle \frac{J_\omega ^4}{3(2\omega J_\omega ^3)}},`$
for $`n`$. The resulting Hamiltonian is of the form :

$$H_F(𝑨,𝝋)=𝝎𝑨+V(𝛀𝑨,𝝋),$$
(13)
where $`𝝎=(\omega ,1)`$ and $`𝛀=(1,0)`$, and the coordinates are $`𝑨=(\mathrm{\Delta }J,L)`$ and $`𝝋=(\theta ,\psi )`$. For small $`F`$, the Kolmogorov-Arnold-Moser (KAM) theorem states that if $`\omega `$ satisfies a Diophantine condition, the invariant torus with frequency $`\omega `$ will persist. In fact, the picture that emerges from numerical simulations is that if $`F`$ (in absolute value) is smaller than some critical threshold $`F_c(\omega )`$ then Hamiltonian (2) \[or equivalently (12)\] has an invariant torus with frequency $`\omega `$. If $`F`$ is larger than this critical threshold, the torus is broken. In order to compute numerically the critical function $`\omega F_c(\omega )`$, we use the renormalization method described briefly below.

#### II.3.2 Renormalization method

Without restriction (up to a rescaling of time), we assume that $`\omega [0,1]`$. The renormalization relies upon the continued fraction expansion of the frequency $`\omega `$ of the torus:

$$\omega =\frac{1}{a_0+\frac{1}{a_1+\mathrm{}}}[a_0,a_1,\mathrm{}].$$

This transformation will act within a space of Hamiltonians $`H`$ of the form

$$H(𝑨,𝝋)=𝝎𝑨+V(𝛀𝑨,𝝋),$$
(14)
where $`𝛀=(1,\alpha )`$ is a vector not parallel to the frequency vector $`𝝎=(\omega ,1)`$. We assume that Hamiltonian (14) satisfies a non-degeneracy condition : If we expand $`V`$ into

$$V(𝛀𝑨,𝝋)=\underset{\genfrac{}{}{0pt}{}{𝝂^2}{k0}}{}V_𝝂^{(k)}(𝛀𝑨)^ke^{i𝝂𝝋},$$

we assume that $`V_\mathrm{𝟎}^{(2)}`$ is non-zero. This restriction means that we cannot explore the twistless region by the present renormalization transformation.
The transformation contains essentially two parts : a rescaling and an elimination procedure koch99 .

(1) Rescaling : The first part of the transformation is composed by a shift of the resonances, a rescaling of time and a rescaling of the actions. It acts on a Hamiltonian $`H`$ as $`H^{}=H𝒯`$ (see Ref. chan00a  for details) :

$$H^{}(𝑨,𝝋)=\lambda \omega ^1H(\frac{1}{\lambda }N_a𝑨,N_a^1𝝋),$$
(15)
with

$$\lambda =2\omega ^1(a+\alpha )^2V_\mathrm{𝟎}^{(2)},$$
(16)
and

$$N_a=\left(\begin{array}{cc}a& 1\\ 1& 0\end{array}\right),$$
(17)
and $`a`$ is the integer part of $`\omega ^1`$ (the first entry in the continued fraction expansion). This change of coordinates is a generalized (far from identity) canonical transformation and the rescaling $`\lambda `$ is chosen to ensure a normalization condition (the quadratic term in the actions has a mean value equal to 1/2). For $`H`$ given by Eq. (14), this expression becomes

$$H^{}(𝑨,𝝋)=𝝎^{}𝑨+\underset{𝝂,k}{}V_𝝂^{(k)}(𝛀^{}𝑨)^ke^{i𝝂𝝋},$$
(18)
where

$`𝝎^{}=(\omega ^{},1)\text{ with }\omega ^{}=\omega ^1a,`$ (19)
$`𝛀^{}=(1,\alpha )\text{ with }\alpha ^{}=1/(a+\alpha ),`$ (20)
$`V_𝝂^{(k)}=r_kV_{N𝝂}^{(k)}\text{ with }r_k=(1)^k2^{1k}\omega ^{k2}(a+\alpha )^{2k}\left(V_\mathrm{𝟎}^{(2)}\right)^{1k}.`$ (21)
We notice that the frequency of the torus is changed according to the Gauss map :

$$\omega \omega ^{}=\omega ^1\left[\omega ^1\right],$$
(22)
where $`\left[\omega ^1\right]`$ denotes the integer part of $`\omega ^1`$. Expressed in terms of the continued fraction expansion, it corresponds to a shift to the left of the entries

$$\omega =[a_0,a_1,a_2,\mathrm{}]\omega ^{}=[a_1,a_2,a_3,\mathrm{}].$$

(2) Elimination : The second step is a (connected to identity) canonical transformation $`𝒰`$ that eliminates the non-resonant modes of the perturbation in $`H^{}`$. Following Ref. chan00a , we consider the set $`I^{}^2`$ of non-resonant modes as the set of integer vectors $`𝝂=(\nu _1,\nu _2)`$ such that $`|\nu _2|>|\nu _1|`$. The canonical transformation $`𝒰`$ is such that $`H^{\prime \prime }=H^{}𝒰`$ does not have any non-resonant mode, i.e. it is defined by the following equation :

$$𝕀^{}(H^{}𝒰)=0,$$
(23)
where $`𝕀^{}`$ is the projection operator on the non-resonant modes; it acts on a Hamiltonian (14) as :

$$𝕀^{}H=\underset{\genfrac{}{}{0pt}{}{𝝂I^{}}{k0}}{}V_𝝂^{(k)}(𝛀𝑨)^ke^{i𝝂𝝋}.$$

We solve Eq. (23) by a Newton method following Refs. chan98b ; chan00a .

Thus the renormalization acts on $`H`$ for a torus with frequency $`\omega `$ as $`H^{\prime \prime }=(H)=H𝒯𝒰`$ for a torus with frequency $`\omega ^{}`$.
The critical thresholds are obtained by iterating the renormalization transformation $``$. The main conjecture of the renormalization approach is that if the torus exists for a given Hamiltonian $`H`$, the iterates $`^nH`$ of the renormalization map acting on $`H`$ converge to some integrable Hamiltonian $`H_0`$. This conjecture is supported by analytical results in the perturbative regime koch99 ; koch99b , and by numerical results chan00d ; chan01a ; chan01R . For a one-parameter family of Hamiltonians $`\{H_F\}`$, the critical amplitude of the perturbation $`F_c(\omega )`$ is determined by the following conditions :

$`^nH_F\underset{n\mathrm{}}{}H_0(𝑨)=𝝎𝑨+{\displaystyle \frac{1}{2}}(𝛀𝑨)^2\text{ for }F<F_c(\omega ),`$ (24)
$`^nH_F\underset{n\mathrm{}}{}\mathrm{}\text{ for }F>F_c(\omega ).`$ (25)
The critical threshold $`F_c(\omega )`$ is determined by a bisection procedure.
In order to obtain the value $`F_m`$ of chaos threshold between resonances $`m`$:1 ans $`m`$+1:1, we vary $`\omega `$ between $`1/(m+1)`$ and $`1/m`$. The critical threshold is given by $`F_m=\mathrm{max}_{\omega [1/(m+1),1/m]}F_c(\omega )`$.

## III Numerical computation of chaos thresholds

### III.1 Chaos thresholds in the negative twist region

Figure 3 represents a typical plot of the critical function $`\omega F_c^{}(\omega )`$ in the negative twist region obtained by the renormalization method for $`e=1`$ and $`\omega _c=0.3`$. This function vanishes at (at least) all rational values of the frequency $`\omega `$ since all tori with rational frequency are broken as soon as the field is turned on. The condition of existence of a torus with frequency $`\omega `$ obtained from the integrable case (7) is $`\omega >\sqrt{2}\omega _cV`$, which is in that case $`\omega >0.67`$. We expect the tori with frequency belonging to \[0,0.67\] to be broken by collision with invariant tori in the positive twist region before this value of the field $`\omega _c`$, as it is the case in Hamiltonian (7).

From Fig. 3, we deduce that for $`F>F_c=0.023`$, no invariant tori are left in this region of phase space. The frequency of the last invariant torus is equal to $`(\gamma +1)/(2\gamma +1)0.7236`$ where $`\gamma =(\sqrt{5}1)/2`$. This value of the frequency of the last invariant torus varies with the parameters $`e`$ and $`\omega _c`$ (for $`\omega _c=0`$, see Fig. 2 of Ref. chan02a ). The main feature of this frequency is that it remains noble as the parameters $`e`$ and $`\omega _c`$ are varied, in the sense that the tail of the continued fraction expansion of this frequency is a sequence of 1, or equivalently this frequency expresses like $`(a\gamma +b)/(c\gamma +d)`$ where $`a`$, $`b`$, $`c`$, $`d`$ are integers such that $`adbc=\pm 1`$.

We compute chaos thresholds between two successive primary resonances by two methods : the 2/3-rule criterion and the renormalization map. We compute the critical value of the field $`F`$ for which there is no longer any invariant torus between resonances 1:1 and 2:1, located at $`J_1^{}`$ and $`J_2^{}`$ respectively, as a function of the magnetic field $`\omega _c`$ and the parameter $`e`$ in the negative twist region, i.e. $`F_1^{}(e,\omega _c)=\mathrm{max}_{\omega [1/2,1]}F_c^{}(\omega ,e,\omega _c)`$. Figure 4 represents two computations : for $`e=0.5`$ and for $`e=1`$.

For $`\omega _c=0`$, we obtain the values that have been obtained in Ref. chan02a . The 2/3-rule criterion gives a very good approximation for low values of the field $`\omega _c`$ (typically for $`\omega _c[0,0.15]`$). For small $`\omega _c`$, we develop the critical function $`F_1^{}(\omega _c)`$ given by Eq. (11). The corrections to $`F_c^{}`$ due to the magnetic field are of order $`\omega _c^2`$ since $`J_n^{}=n^{1/3}+O(\omega _c^2).`$

We notice that at some value of $`\omega _c`$, the curves $`F_1^{}(e,\omega _c)`$ given by the Chirikov criterion decrease sharply to zero. We have seen that the approximate condition of existence of the resonance $`m`$:1, derived from the integrable case (8), is $`\omega _c<1/(\sqrt{2}mV)`$. For $`e=1`$ and $`m=2`$, this condition is $`\omega _c<0.23`$ (and for $`e=0.5`$ this condition becomes $`\omega _c<0.30`$). This means that the criterion derived in Sec. II.2 is inapplicable for $`\omega _c`$ larger than these values ($`J_2^{}`$ becomes complex). One way of extending this criterion for larger values of $`\omega _c`$ would be to consider the overlap of higher order resonances. From the renormalization map results, we see that after 2:1 resonance disappear, there is an increase of stability which can be explained by the fact that the invariant tori in this region of phase space are no longer deformed by this neighboring resonance. We expect this region to be fully broken before 1:1 disappears (i.e. all the region between 1:1 and 2:1 primary resonances has disappeared), and this happens at $`\omega _c0.45`$ for $`e=1`$, and at $`\omega _c0.6`$ for $`e=0.5`$. These results are consistent with the sharp decreases of $`F_c`$ found by renormalization on Fig. 4.

Furthermore, we investigate the chaos thresholds for the different primary resonances $`m`$:1 by the Chirikov criterion. Using Eq. (11), we compute the resonance overlap value of $`m`$:1 and $`m`$+1:1 for increasing $`m=1,\mathrm{},10`$. For a fixed eccentricity $`e=1`$, Fig. 5 shows the resonance overlap values as a function of the field $`\omega _c`$

From this figure it emerges that the dynamics in that region is mainly insensitive to the magnetic field for low values of this field. The sharp decreases of the curves $`F_m^{}`$ are due to the disappearance of one of the primary resonances from which the resonance overlap is computed. At these values and for larger $`\omega _c`$, we expect stability enhancement due to the disappearance of this primary resonance. For $`\omega _c1/(\sqrt{2}mV)`$, the critical curve $`F_m^{}`$ sharply decreases to zero due to the disappearance of all the region between primary resonances $`m`$:1 and $`m`$+1:1 into the twistless region. For instance, for $`m=3`$, the critical curve $`F_3^{}(\omega _c)`$ increases slowly with $`\omega _c`$ for $`\omega _c[0,0.11]`$, then for $`\omega _c[0.11,0.15]`$, we expect stabilization enhancement by the magnetic field, and for $`\omega _c`$ larger than 0.15, the region of phase space between 3:1 and 4:1 resonances disappears into the twistless region and $`F_3^{}(\omega _c)`$ vanishes.

The region between resonances of low order ($`m`$ small) seem to be more stable for $`e=1`$. However, this feature varies with the parameter $`e`$ as it has been observed in Ref. chan02a . For $`e[0,0.8]`$ and for $`\omega _c=0`$, the regions near $`m`$:1 with large $`m`$ become very stable (and this stabilization is increased by the field for low values of $`\omega _c`$) and the diffusion of the trajectories is very limited. Therefore, the orbits with eccentricity close to one are the easiest to diffuse in a broad range of phase space. In particular, these orbits can ionize more easily than medium eccentricity ones. This feature is expected to hold for small $`\omega _c`$ (typically for $`\omega _c0.05`$) even if the diffusion is now limited in the negative twist region since we expect the twistless region to be very stable and act as a barrier in phase space cast96 . Figure 6 displays the values of $`F_m^{}(e)`$ for $`e[0,1]`$, and for two values of $`\omega _c`$ : $`\omega _c=0.01`$ and $`\omega _c=0.1`$.

For $`\omega =0.1`$, since a large portion of phase space between primary resonances with $`m`$ large has disappeared into the twistless region, the orbits with medium eccentricity become as easy to ionize as the ones with eccentricity close to one in the remaining part of the negative twist region of phase space. For eccentricities close to one, the critical threshold is dominated by the chaos threshold between resonances 1:1 and 2:1.

In summary, the effect of the magnetic field $`\omega _c`$ is to stabilize the dynamics in the negative twist region by successively breaking up primary resonances. Therefore, between two succesive resonances with $`m`$ small, the effect of the magnetic field is expected to be smooth in the sense that the critical thresholds are only slightly changed (increased) from the chaos thresholds in the absence of magnetic field up to some critical value of $`\omega _c`$ where collisions of periodic orbits occur. Investigating this region by Chirikov resonance overlap yields values which are qualitatively and quantitatively accurate with respect to the renormalization results for $`\omega _c[0,0.1]`$.
For small values of $`\omega _c`$, the qualitative behavior concerning diffusion of trajectories is expected to be the same as the one in the absence of magnetic field even if the diffusion coefficients may be smaller due to the limited region of phase space and the stability enhancement due to the magnetic field. Increasing $`\omega _c`$ makes the orbits with medium eccentricity easier to diffuse in a more and more reduced region of phase space.

### III.2 Chaos thresholds in the positive twist region

The classical dynamics in the positive twist region of phase space is essential for the stochastic ionization process since it contains the region where the action $`J`$ is large, anf for large $`\omega _c`$ this region is the predominant region since the negative twist region shrinks to the singularity $`J=0`$.

Figure 7 show a typical critical function $`\omega F_c^+(\omega )`$ obtained by the renormalization method for $`\omega _c=0.3`$ and for $`e=1`$. This figure is analogous to Fig. 3.

From Fig. 7, we deduce that for $`F>F_c=0.021`$, no invariant tori are left in this region of phase space. For these values of the parameters $`e`$ and $`\omega _c`$, the condition of existence of invariant tori $`\omega \sqrt{2}\omega _cV`$ obtained from the integrable case (7), is $`\omega 0.67`$. We notice that there is a large chaotic zone for $`\omega [0.67,0.8]`$ (the set of $`\omega `$ where $`F_c^+(\omega )`$ is very small). The diffusion of trajectories throughout the $`J`$ large region is prevented by the small set of invariant tori with frequency $`\omega [0.85,0.90]`$ for intermediate values of the amplitude of the microwaves $`F_c[0.015,0.02]`$. The frequency of the last invariant torus to survive in this region of phase space is approximately equal to 0.87. The frequency of the last invariant torus in the positive twist region fluctuates in a very erratic way as one varies $`\omega _c`$. However we observed that this frequency is always between 0.85 and 0.90 (close to the 1:1 resonance).

The changes induced by the magnetic field are stronger in the positive twist region. We compute the critical thresholds between 1:1 and 2:1 resonances in the positive twist region. Figure 8 represents the results obtained by renormalization and by the 2/3-rule criterion. In this region, the critical thresholds determined by 2/3-rule criterion are well below the ones given by the renormalization for significant values of the field.

Figure 9 represents the critical threshold for the break-up of invariant tori with frequency $`\gamma =(\sqrt{5}1)/2`$. In the negative twist region, the golden mean torus is slightly stabilized by the magnetic field. In the positive twist region, the influence of the field is very strong : first, the torus is created then stabilized to high values of the field, then it disappears. This figure shows that in contrast to the negative twist region, the influence of the magnetic field in the positive twist region is very strong. The curve for the positive twist torus appears to have non-smooth variations conversely to the negative twist torus like for instance for $`\omega _c0.065`$. We notice that for the golden mean tori, there is no break-up by collision since the one located in the positive twist region is broken before the expected collision.

For small values of the field $`\omega _c`$, we expand the threshold given by Chirikov criterion (11). Since the resonances located at $`J_n^+`$ do not exist in the absence of the field $`\omega _c`$, we expect $`F_c^+(m,e,\omega _c)`$ to vanish at $`\omega _c=0`$. The corrections to $`F_c^+`$ due to the magnetic field are obtained using the following expansion of $`J_n^+=2^{1/3}n^{1/3}V^{2/3}\omega _c^{2/3}+O(1)`$ Thus we have : $`F_m^+(e,\omega _c)=\alpha _m(e)\omega _c^{2/3}[1+O(\omega ^{2/3})]`$, where $`\alpha _m`$ are some constants depending on $`e`$. This means that the increase of stabilization is very sharp (with infinite slope) for low values of $`\omega _c`$.

Furthermore, using Eq. (11), we compute the 2/3-rule resonance overlap value of $`m`$:1 and $`m`$+1:1 resonances for increasing $`m=1,\mathrm{},10`$. For a fixed eccentricity $`e=1`$ , Fig. 10 shows the resonance overlap values as a function of the field $`\omega _c`$ in the positive twist region.

What emerges from this figure is that this region of phase space is strongly stabilized by the magnetic field, and that the curve $`m=1`$ dominates. However, this feature may vary with the parameter $`e`$ as it has been observed in Ref. chan02a  and in Sec. III.1 for the negative twist region.

These curves are similar to Fig. 6. As a result, increasing the magnetic field makes intermediate eccentricity orbits as easy to ionize (diffusion through a part of phase space where $`J`$ is large) as the ones with eccentricity close to one. This situation is different from the situation without magnetic field.

The overall effect of the magnetic field in the positive twist region is basically the same as the one in the negative twist region, i.e., it breaks up the primary resonances and fills up the remaining region by very stable invariant tori, preventing the diffusion of trajectories throughout phase space.

## IV Conclusion

We have analyzed the classical phase space of the hydrogen atom in crossed magnetic and microwave fields in the high frequency regime. Useful information about the dynamics is provided by the analysis of an integrable part of the Hamiltonian. Accurate information about chaos threshold is obtained by using the renormalization map and the 2/3-rule Chirikov criterion. The global effect of the magnetic field is to stabilize the dynamics and consequently reducing the diffusion of trajectories and the stochastic ionization process.