# Improved upper bounds for the Heilbronn’s Problem for $k$-gons

Rishikesh Gajjala Indian Institute of Science, Bengaluru,
<EMAIL_ADDRESS>Jayanth Ravi

###### Abstract

The Heilbronn triangle problem asks for the placement of $n$ points in a unit
square that maximizes the smallest area of a triangle formed by any three of
those points. In $1972$, Schmidt considered a natural generalization of this
problem. He asked for the placement of $n$ points in a unit square that
maximizes the smallest area of the convex hull formed by any four of those
points. He showed a lower bound of $\Omega(n^{-3/2})$, which was improved to
$\Omega(n^{-3/2}\log{n})$ by Leffman.

A trivial upper bound of $3/n$ could be obtained and Schmidt asked if this can
be improved asymptotically. However, despite several efforts, no asymptotic
improvement over the trivial upper bound was known for the last $50$ years,
and the problem started to get the tag of being notoriously hard. Szemerédi
posed the question of whether one can, at least, improve the constant in this
trivial upper bound. In this work, we answer this question by proving an upper
bound of $2/n+o(1/n)$. We also extend our results to any convex hulls formed
by $k\geq 4$ points.

##  1 Introduction

Given a constant $k\geq 3$ and a set
$\mathcal{P}=\\{P_{1},P_{2},\ldots,P_{n}\\}$ of $n\geq k$ points on the unit
square $[0,1]^{2}$, let $\mathcal{A}_{k}(P)$ be the area of the smallest
convex hull among all convex hulls determined by subsets of $k$ points in
$\mathcal{P}$. The supremum value of $\mathcal{A}_{k}(P)$ over all choices of
$\mathcal{P}$ is denoted by $\Delta_{k}(n)$. The Heilbronn triangle problem
asks for the value of $\Delta_{3}(n)$.

The Heilbronn triangle problem is one of the fundamental problems in discrete
geometry and discrepancy theory and has a rich history. Paul Erdős proved that
$\Delta_{3}(n)=\Omega\left(\frac{1}{n^{2}}\right)$. This was believed to be
the upper bound for some time until Komlós, Pintz and Szemerédi [KPS82] proved
that

$\Delta_{3}(n)=\Omega\left(\frac{\log{n}}{n^{2}}\right)$

Over a series of works, the upper bounds were improved by Roth [Rot72, Rot73,
Rot76, Rot51] and Schmidt [Sch72]. The current best-known upper bound is due
to Komlós, Pintz and Szemerédi [KPS81]

$\Delta_{3}(n)=\mathcal{O}\left(\frac{2^{c\sqrt{\log{n}}}}{n^{8/7}}\right)$

This has been recently claimed to be improved by Cohen, Pohoata and Zakharov
[CPZ23] to $\mathcal{O}\left(n^{-8/7-1/2000}\right)$.

There has also been work on several variants of this problem. Jiang, Li and
Vitany [JLV02] and Benevides, Hoppen, Lefmann and Odermann [BHLO23] studied
the case in which the points were randomly distributed. The problem was also
explored in higher dimensions by placing $n$ points in $d$-dimensional unit
cubes $[0,1]^{d}$ instead of a unit square [Bar01, Bar04, BS07, Bra05, Lef03,
LS02].

Schmidt asked about the value of $\Delta_{k}(n)$ and proved that
$\Delta_{4}(n)=\Omega\left(\frac{1}{n^{1.5}}\right)$ [Sch72]. Bertraln-
Kretzberg, Hofmeister and Lefmann generalized this result to $k$-gons by
proving that $\Delta_{k}(n)=\Omega\left(\frac{1}{n^{\frac{k-1}{k-2}}}\right)$
[BHL00]. This was improved by Lefmann[Lef08] to

$\Delta_{k}(n)=\Omega\left(\frac{{(\log{n})}^{1/k-2}}{n^{1+\frac{1}{k-2}}}\right)$

##  2 Our results

A trivial upper bound of $\Delta_{4}(n)\leq\dfrac{3}{n}$ can be obtained by
subdividing the unit square into squares of side length $\sqrt{\dfrac{3}{n}}$
using the pigeonhole principle. However, despite several efforts to improve
this upper bound (asymptotically) since it was posed in $1972$, there has been
no progress! Szemerédi asked if at least the constants in this upper bound can
be improved [Sze22]. In this work, we answer this question by proving the
following theorem.

######  Theorem 2.1.

$\Delta_{4}(n)\leq\dfrac{2}{n}+o\left(\dfrac{1}{n}\right)$

We also generalize our result to general $k$-gons for any constant $k\geq 4$.

######  Theorem 2.2.

$\Delta_{k}(n)\leq\dfrac{k-2}{n}+o\left(\dfrac{1}{n}\right)$

##  3 Convex quadrilaterals: Proof of Theorem 2.1

We solve a more general problem by having the points on a unit rectangle
(instead of a unit square). Given a set
$\mathcal{P}=\\{P_{1},P_{2},\ldots,P_{n}\\}$ of $n\geq 3$ points on the unit
rectangle $[0,d]\times[0,d^{-1}]$ and $k\leq n$, let
$\mathcal{A^{\prime}}_{k}(P)$ be the minimum area of the convex hull
determined by a set of $k$ points in $\mathcal{P}$ for any $0<d\leq 1$. The
supremum value of $\mathcal{A^{\prime}}_{k}(P)$ over all choices of
$\mathcal{P}$ is denoted by $\Delta_{k}^{\prime}(n)$. It is easy to see that
by definition

$\Delta_{k}(n)\leq\Delta_{k}^{\prime}(n)$

For $k=4$, when there are $n$ points, in the argument to obtain a trivial
bound of $\Delta_{4}^{\prime}(n)\leq\dfrac{3}{n-1}$, we partition the unit
rectangle into at most ${(n-1)}/{3}$ smaller rectangles. This would guarantee
that there exists a small rectangle containing at least $4$ points. Naturally,
one can also extend this idea to make sure there are at most
$\left(\dfrac{n-1}{n^{\prime}-1}\right)$ smaller rectangles and force
$n^{\prime}\geq 4$ points into one rectangle. This gives us a relation between
$\Delta_{4}^{\prime}(n)$ and $\Delta_{4}^{\prime}(n^{\prime})$, which is
formalized in Observation 3.1.

######  Observation 3.1.

For $4\leq n^{\prime}\leq n$,

$\Delta_{4}^{\prime}(n)\leq{\left\lfloor{\dfrac{n-1}{n^{\prime}-1}}\right\rfloor}^{-1}\Delta_{4}^{\prime}(n^{\prime})$

###### Proof:

Partition the unit area into a grid with
${\left\lfloor{\dfrac{n-1}{n^{\prime}-1}}\right\rfloor}$ rectangles of area
${\left\lfloor{\dfrac{n-1}{n^{\prime}-1}}\right\rfloor}^{-1}$. Since there are
$n$ points and
${\left\lfloor{\dfrac{n-1}{n^{\prime}-1}}\right\rfloor}\leq\dfrac{n-1}{n^{\prime}-1}$
rectangles, by the pigeonhole principle, one of the smaller rectangles (with
their boundary included) has at least $n^{\prime}$ points. Therefore, there
always exists $n^{\prime}$ points within an area at most
${\left\lfloor{\dfrac{n-1}{n^{\prime}-1}}\right\rfloor}^{-1}$. It now follows
by a scaling argument that there exist four points within an area at most
${\left\lfloor{\dfrac{n-1}{n^{\prime}-1}}\right\rfloor}^{-1}\Delta_{4}^{\prime}(n^{\prime})$.

When $n^{\prime}=4$, this gives the trivial bound of
$\Delta_{4}^{\prime}(n)\leq\dfrac{3}{n-3}\approx 3/n$ as expected. We can do
slightly better by tuning the value of $n^{\prime}$ to be $6$. We start with
finding the exact value of $\Delta_{4}^{\prime}(6)$.

######  Observation 3.2.

$\Delta_{4}^{\prime}(6)=1/2$

###### Proof:

Let $C$ be the centre of the rectangle $[0,d]\times[0,d^{-1}]$ and
$\mathcal{P}$ be a set of any six points in $[0,d]\times[0,d^{-1}]$. Pick an
arbitrary point $P\in\mathcal{P}$ and extend the line segment $PC$ into a
line. The extended line $PC$ cuts the rectangle $[0,d]\times[0,d^{-1}]$ into
two convex parts, and by symmetry, both these parts have the same area, i.e.,
$1/2$. From the pigeonhole principle, at least $3$ of the remaining $5$ points
lie on one side of the extended line $PC$. Therefore, $4$ points (including
$P$) exist, which are contained in a convex shape whose area is $1/2$.
Therefore, $\Delta_{4}^{\prime}(6)\leq 1/2$.

(0,0)(0,$0.5d^{-1}$)(0,$d^{-1}$)($d$,0)($d,0.5d^{-1}$)($d,d^{-1}$) Figure 3.1:
$\Delta_{4}^{\prime}(6)\geq 1/2$

It may be noted the bound of $1/2$ is achieved in Figure 3.1. Therefore
$\Delta_{4}^{\prime}(6)\geq 1/2$

######  Corollary 3.3.

$\Delta_{4}(n)\leq\dfrac{1}{2}\left\lfloor{\dfrac{n-1}{5}}\right\rfloor^{-1}\leq\dfrac{2.5}{n-5}$

###### Proof:

By substituting $n^{\prime}=6$ in Observation 3.1 and using Observation 3.2,
we get

$\Delta_{4}^{\prime}(n)\leq{\left\lfloor{\dfrac{n-1}{5}}\right\rfloor}^{-1}\Delta_{4}^{\prime}(6)=\dfrac{1}{2}{\left\lfloor{\dfrac{n-1}{5}}\right\rfloor}^{-1}\leq\dfrac{2.5}{n-5}$

We will now extend the idea of Observation 3.2 to all $n$ of the form
$2^{s}+2$ for any $s\geq 2$, i.e., for all $n$ of such form, we prove
$\Delta_{4}(n)\leq\dfrac{2}{n-2}$ in Theorem 3.6 using Observation 3.4 and
Lemma 3.5.

######  Observation 3.4.

For any point $P$ in a convex polygon $C$, there exists a line through $P$
which partitions $C$ into two halves of equal area.

###### Proof:

Let $A$ be the area of $C$. Pick any arbitrary line $L$ through $P$ and let
the areas of the convex polygons on both of its sides of $L$ be $L_{1}$ and
$L_{2}$ such that $L_{1}\leq A/2\leq L_{2}$. By rotating the line by
${180}^{\circ}$ around $P$, we get $L_{2}\leq 1/2\leq L_{1}$. Since $L_{1}$
changes continuously as a function of the angle of rotation $\theta$, by the
intermediate value theorem, it must have achieved $1/2$ in between for some
$\theta$.

######  Lemma 3.5.

If a convex polygon of area $\Delta$ has $2^{i}\cdot\beta+2$ points, then
there is a convex polygon of area $\Delta/2$ which contains at least
$2^{i-1}\cdot\beta+2$ points.

###### Proof:

Pick one of the $2^{i}\cdot\beta+2$ points arbitrarily, say $P$. From
Observation 3.4, there is a line $L$ through $P$ cutting the polygon into two
halves. Note that by the pigeonhole principle, one of the halves would have at
least $2^{i-1}\cdot\beta+1$ points. By including $P$ in this region, we get a
convex polygon of area $\Delta/2$, containing at least $2^{i-1}\cdot\beta+2$
points.

We first introduce some new notation. Given a set
$\mathcal{Q}=\\{Q_{1},Q_{2},\ldots,Q_{n}\\}$ of $n\geq 3$ points on an
arbitrary convex object $\mathcal{C}$ of unit area, let
$\mathcal{A}_{k}^{\mathcal{C}}(Q)$ be the minimum area of the convex hull
determined by some $k$ points in $\mathcal{Q}$. The supremum value of
$\mathcal{A}_{k}^{\mathcal{C}}(Q)$ over all choices of $\mathcal{Q}$ is
denoted by $\Delta^{\mathcal{C}}_{k}(n)$. Let $\Delta_{k}^{\prime\prime}(n)$
denote the supremum value of $\Delta^{\mathcal{C}}_{k}(n)$ over all convex
objects of unit area. It is easy to see that by definition

$\Delta_{k}(n)\leq\Delta_{k}^{\prime}(n)\leq\Delta_{k}^{\prime\prime}(n)$

######  Theorem 3.6.

If $n=2^{s}+2$ for some integer $s>0$, then

$\Delta_{4}(n)\leq\Delta_{4}^{\prime}(n)\leq\dfrac{1}{2^{s-1}}=\dfrac{2}{n-2}$

###### Proof:

Let $\mathcal{Q}=\\{Q_{1},Q_{2},\ldots,Q_{n}\\}$ be a set of $n$ points. We
will prove a stronger statement of

$\Delta^{\prime\prime}_{4}(n)\leq\dfrac{1}{2^{s-1}}$

######  Observation 3.7.

For every $i\in[0,s-1]$, there exists a convex polygon of area at most
$\dfrac{1}{2^{i}}$ which has at least $2^{s-i}+2$ points from $\mathcal{Q}$.

###### Proof:

We prove this by induction on $i$. When $i=0$, the claim is true by
definition. Suppose there exists a convex polygon of area at most
$\dfrac{1}{2^{i}}$ with at least $2^{s-i}+2$ points from $\mathcal{Q}$, then
from Lemma 3.5, for $i<s-1$, there exists a convex polygon of area at most
$\dfrac{1}{2^{i+1}}$ with at least $2^{s-(i+1)}+2$ points from $\mathcal{Q}$

By substituting $i=s-1$ in Observation 3.7, Theorem 3.6 follows.

One may note that this would give an upper bound of $\approx 2/n$ for many
arbitrarily large $n$ of the form $2^{s}+2$. However, there are also several
arbitrarily large $n$ of the form, say, $2^{s}+1$, for which this bound is not
useful. We fix this using Observation 3.1.

######  Corollary 3.8.

$\Delta_{4}(n)\leq\dfrac{2}{n}+o\left(\dfrac{1}{n}\right)$

###### Proof:

For every $n$, there exists some $i$ such that

$2^{i-1}+2\leq n<2^{i}+2$

Let $n^{\prime}=2^{\left\lceil{*}\right\rceil{0.5i}}+2$. From Theorem 3.6,

$\Delta_{4}^{\prime}(n^{\prime})\leq\dfrac{2}{n^{\prime}-2}<\dfrac{2}{n^{\prime}-1}$

${\left\lfloor{\dfrac{n-1}{n^{\prime}-1}}\right\rfloor}>\frac{n-1}{n^{\prime}-1}-1\geq\frac{n-n^{\prime}}{n^{\prime}-1}$

From Observation 3.1,

$\Delta_{4}^{\prime}(n)\leq\Delta_{4}^{\prime}(n^{\prime}){\left\lfloor{\dfrac{n-1}{n^{\prime}-1}}\right\rfloor}^{-1}\leq\dfrac{2}{n^{\prime}-1}\cdot\frac{n^{\prime}-1}{n-n^{\prime}}=\dfrac{2}{n-n^{\prime}}$
$=\dfrac{2}{n}+\dfrac{2n^{\prime}}{n(n-n^{\prime})}=\dfrac{2}{n}+\mathcal{O}\left(\dfrac{1}{n^{1.5}}\right)$

(0,0)(0,$0.5d^{-1}$)(0,$d^{-1}$)($0.5d$,0)($0.5d$,$0.5d^{-1}$)($0.5d$,$d^{-1}$)($d$,0)($d,0.5d^{-1}$)($d,d^{-1}$)
Figure 3.2: $\Delta_{4}^{\prime}(9)\geq 1/4$

From Figure 3.2, it is easy to see that $\Delta_{4}^{\prime}(9)\geq 1/4$. We
conjecture that the other direction is also true, i.e.,
$\Delta_{4}^{\prime}(9)\leq 1/4$

######  Conjecture 3.9.

$\Delta_{4}^{\prime}(9)=1/4$

If true, Conjecture 3.9 would directly imply that
$\Delta_{4}(n)\leq\dfrac{2}{n-8}$ from Observation 3.1 by picking $n^{\prime}$
to be $9$. We also note that finding the exact values of $\Delta_{k}{(n)}$ is
of independent interest and has been well studied for $k=3$ [CY02, DZ22,
ZC08].

One may note that our analysis will extend to general convex figures with unit
area (instead of just unit squares), i.e.,

$\Delta^{\prime\prime}_{4}(n)\leq\dfrac{2}{n}+o\left(\dfrac{1}{n}\right)$

##  4 Convex $k$-gons: Proof of Theorem 2.2

In this section, we extend our proofs for $k=4$ to give upper bounds on
$\Delta_{k}(n)$ for any constant $k$.

######  Proposition 4.1 (Analogue of Observation 3.1).

For $n^{\prime}\leq n$,

$\Delta_{k}^{\prime}(n)\leq{\left\lfloor{\dfrac{n-1}{n^{\prime}-1}}\right\rfloor}^{-1}\Delta_{k}^{\prime}(n^{\prime})$

###### Proof:

The proof of Proposition 4.1 directly extends from Observation 3.1.

######  Theorem 4.2 (Analogue of Theorem 3.6).

If $n=2^{s}(k)-2^{s+1}+2$ for some integer $s>0$, then

$\Delta_{k}(n)\leq\Delta_{k}^{\prime}(n)\leq\dfrac{1}{2^{s}}=\dfrac{k-2}{n-2}$

###### Proof:

Let $\mathcal{Q}=\\{Q_{1},Q_{2},\ldots,Q_{n}\\}$ be a set of $n$ points. We
will prove a stronger statement of

$\Delta^{\prime\prime}_{k}(n)\leq\dfrac{1}{2^{s}}$

######  Observation 4.3.

For every $i\in[0,s]$, there exists a convex polygon of area at most
$\dfrac{1}{2^{i}}$ with at least $2^{s-i}(k)-2^{s+1-i}+2$ points from
$\mathcal{Q}$.

###### Proof:

We prove this by induction on $i$. When $i=0$, the claim is true by
definition. Suppose there exists a convex polygon of area at most
$\dfrac{1}{2^{i}}$ with at least $2^{s-i}(k)-2^{s+1-i}+2$ points from
$\mathcal{Q}$, then from Lemma 3.5, for $i<s$, there exists a convex polygon
of area at most $\dfrac{1}{2^{i+1}}$ with at least
$2^{s-(i+1)}(k)-2^{s+1-(i+1)}+2$ points from $\mathcal{Q}$.

By substituting $i=s$ in Observation 4.3, Theorem 4.2 follows.

######  Corollary 4.4.

$\Delta_{k}(n)\leq\dfrac{k-2}{n}+o\left(\dfrac{1}{n}\right)$

###### Proof:

For every $n$, there exists some $i$ such that

$2^{i-1}(k)-2^{i}+2\leq n<2^{i}(k)-2^{i+1}+2$

Let
$n^{\prime}=2^{\left\lceil{*}\right\rceil{0.5i}}(k)-2^{\left\lceil{*}\right\rceil{0.5i}+1}+2$.
From Theorem 4.2,

$\Delta_{k}^{\prime}(n^{\prime})\leq\dfrac{k-2}{n^{\prime}-2}<\dfrac{k-2}{n^{\prime}-1}$

${\left\lfloor{\dfrac{n-1}{n^{\prime}-1}}\right\rfloor}>\frac{n-1}{n^{\prime}-1}-1\geq\frac{n-n^{\prime}}{n^{\prime}-1}$

From Proposition 4.1,

$\Delta_{k}^{\prime}(n)\leq\Delta_{k}^{\prime}(n^{\prime}){\left\lfloor{\dfrac{n-1}{n^{\prime}-1}}\right\rfloor}^{-1}\leq\dfrac{k-2}{n^{\prime}-1}\cdot\frac{n^{\prime}-1}{n-n^{\prime}}=\dfrac{k-2}{n-n^{\prime}}$
$=(\dfrac{1}{n}+\dfrac{n^{\prime}}{n(n-n^{\prime})})\cdot(k-2)=\dfrac{k-2}{n}+\mathcal{O}\left(\dfrac{1}{n^{1.5}}\right)$

######  Conjecture 4.5.

$\Delta_{k}^{\prime}(\alpha(k-1))=\dfrac{1}{2(\alpha-1)}$ for $\alpha\leq k-1$

One can prove that
$\Delta_{k}^{\prime}(\alpha(k-1))\geq\dfrac{1}{2(\alpha-1)}$ by placing the
points in the corners of a $\alpha\times(k-1)$ grid. We conjecture that this
is, in fact, the optimal placement. When $\alpha=2$, this is indeed true by
Lemma 3.5. Conjecture 3.9 is a special case of Conjecture 4.5 when $k=4$. If
true, Conjecture 4.5 would imply an upper bound of $\approx{k}/{(2n)}$

One may notice that Conjecture 4.5 can not be extended to $\alpha\leq k-1$, as
at least $k$ points become collinear in that case.

## Acknowledgements

RG thanks Endre Szemerédi for introducing this problem to him during the $9$th
Heidelberg Laureate Forum. RG thanks Saladi Rahul for his comments on the
manuscript.

## References

  * [Bar01] Gill Barequet  “A Lower Bound for Heilbronn’s Triangle Problem in _d_ Dimensions”  In _SIAM J. Discret. Math._ 14.2, 2001, pp. 230–236  DOI: 10.1137/S0895480100365859
  * [Bar04] Gill Barequet  “The on-line Heilbronn’s triangle problem”  In _Discret. Math._ 283.1-3, 2004, pp. 7–14  DOI: 10.1016/J.DISC.2004.03.001
  * [BHL00] Claudia Bertram-Kretzberg, Thomas Hofmeister and Hanno Lefmann  “An Algorithm for Heilbronn’s Problem”  In _SIAM J. Comput._ 30.2, 2000, pp. 383–390  DOI: 10.1137/S0097539798348870
  * [BHLO23] Fabricio S Benevides, Carlos Hoppen, Hanno Lefmann and Knut Odermann  “Heilbronn triangle-type problems in the unit square [0, 1] 2”  In _Random Structures & Algorithms_ 62.3 Wiley Online Library, 2023, pp. 585–599
  * [Bra05] Peter Braß  “An Upper Bound for the d-Dimensional Analogue of Heilbronn’s Triangle Problem”  In _SIAM J. Discret. Math._ 19.1, 2005, pp. 192–195  DOI: 10.1137/S0895480103435810
  * [BS07] Gill Barequet and Alina Shaikhet  “The On-Line Heilbronn’s Triangle Problem in d Dimensions”  In _Discret. Comput. Geom._ 38.1, 2007, pp. 51–60  DOI: 10.1007/S00454-007-1323-X
  * [CPZ23] Alex Cohen, Cosmin Pohoata and Dmitrii Zakharov  “A new upper bound for the Heilbronn triangle problem”  In _arXiv preprint arXiv:2305.18253_ , 2023
  * [CY02] Francesc Comellas and J Luis A Yebra  “New lower bounds for Heilbronn numbers”  In _the electronic journal of combinatorics_ , 2002, pp. R6–R6
  * [DZ22] Lydia Dehbi and Zhenbing Zeng  “Heilbronn’s Problem of Eight Points in the Square”  In _Journal of Systems Science and Complexity_ 35.6 Springer, 2022, pp. 2452–2480
  * [JLV02] Tao Jiang, Ming Li and Paul M.. Vitányi  “The average-case area of Heilbronn-type triangles”  In _Random Struct. Algorithms_ 20.2, 2002, pp. 206–219  DOI: 10.1002/RSA.10024
  * [KPS81] János Komlós, János Pintz and Endre Szemerédi  “On Heilbronn’s triangle problem”  In _Journal of the London Mathematical Society_ 2.3 Citeseer, 1981, pp. 385–396
  * [KPS82] János Komlós, János Pintz and Endre Szemerédi  “A lower bound for Heilbronn’s problem”  In _Journal of the London Mathematical Society_ 2.1 Oxford University Press, 1982, pp. 13–24
  * [Lef03] Hanno Lefmann  “On Heilbronn’s Problem in Higher Dimension”  In _Comb._ 23.4, 2003, pp. 669–680  DOI: 10.1007/S00493-003-0040-1
  * [Lef08] Hanno Lefmann  “Distributions of points in the unit square and large k-gons”  In _Eur. J. Comb._ 29.4, 2008, pp. 946–965  DOI: 10.1016/J.EJC.2007.11.021
  * [LS02] Hanno Lefmann and Niels Schmitt  “A Deterministic Polynomial-Time Algorithm for Heilbronn’s Problem in Three Dimensions”  In _SIAM J. Comput._ 31.6, 2002, pp. 1926–1947  DOI: 10.1137/S0097539701395115
  * [Rot51] Klaus F Roth  “On a problem of Heilbronn”  In _Journal of the London Mathematical Society_ 1.3 Wiley Online Library, 1951, pp. 198–204
  * [Rot72] KF Roth  “On a problem of Heilbronn, II”  In _Proceedings of the London Mathematical Society_ 3.2 Oxford University Press, 1972, pp. 193–212
  * [Rot73] KF Roth  “Estimation of the area of the smallest triangle obtained by selecting three out of n points in a disc of unit area”  In _Proc. of Symposia in Pure Mathematics_ 24, 1973, pp. 251–262  AMS Providence
  * [Rot76] KF Roth  “Developments in Heilbronn’s triangle problem”  In _Advances in Mathematics_ 22.3 Elsevier, 1976, pp. 364–385
  * [Sch72] Wolfgang M. Schmidt  “On a Problem of Heilbronn†”  In _Journal of the London Mathematical Society_ s2-4.3, 1972, pp. 545–550  DOI: 10.1112/jlms/s2-4.3.545
  * [Sze22] Endre Szemerédi  “Personal Communication”, 2022
  * [ZC08] Zhenbing Zeng and Liangyu Chen  “On the Heilbronn optimal configuration of seven points in the square”  In _International Workshop on Automated Deduction in Geometry_ , 2008, pp. 196–224  Springer


Axial-Anomaly in Noncommutative QED and Pauli-Villars Regularization

Mohammad Walid AlMasri

Physics Department, Koç University, Sarıyer 34450,Istanbul,Turkey
Emails<EMAIL_ADDRESS><EMAIL_ADDRESS>

Abstract

We calculate the $U(1)$ axial-anomaly in two and four dimensions using a
modified path integral method coupled to a Pauli-Villars regulator field in
the noncommutative QED. Pauli-Villars regularization method provides us with
unambiguous way to connect the modified path integral formalism with
perturbative axial Ward identities at each step of calculations.

###### Contents

  1. 1 Introduction
  2. 2 Noncommutative $U(1)$ Gauge Theory
  3. 3 Pauli-Villars Regularization Method
    1. 3.1 Two-Dimensional Anomaly
    2. 3.2 Four-Dimensional Anomaly
  4. 4 Conclusion
  5. A Fujikawa Method In Four-Dimensions

##  1 Introduction

Quantum anomaly occurs when the conservation laws are valid in the classical
level but it is violated at the quantized level. One basic feature of quantum
field theory is the principle of gauge symmetry. Its breakdown violates the
corresponding conserved current. For example, in some theories of fermions
with chiral symmetry, the quantization may lead to the breaking of chiral
symmetry which in turn implies a non-conserved chiral charges. The axial
anomaly was first studied by Jack Steinberger in 1949 when he calculated the
Feynman Triangle diagram in “Pion-Nucleon” $\pi$-N Model which contains a
$\gamma_{5}$ vertex. in 1963, Johnson pointed out that in a massless
2-dimensional QED one cannot have both the conservation of gauge current and
axial current. It was proved by performing the Feynman triangle diagrams which
made up of one-axial current and two-vector currents containing an UV-
divergence. They found that conservation of vector current implies the
breaking of conversation law for axial current. In 1969 Adler-Bardeen showed
that anomaly is already given by 1-Loop calculations. Thus, axial anomaly
receives no contributions from radiative corrections. This means that due to
loop corrections, the symmetry might be spoiled and this symmetry will not be
the symmetry of effective action.
Non-commutative quantum field theory has been studied extensively in the last
two decades [1, 2, 3]. It has been shown that an UV/IR mixing arises in
noncommutative theories. It also emerges in the context of perturbative string
theory in the presence of D-branes with non-vanishing $B$-Field background [4,
5].
The $U(1)$ axial anomaly in the non-commutative QED has been calculated and it
is well understood, [6, 7]. In this work, We will calculate the $U(1)$ axial
anomaly in two and four dimensions using Pauli-Villars regularization method.

The layout of this paper as follows, In first section of this work, we shall
give a brief review on Non-commutative $U(1)$ gauge theory. In the second
section, we will calculate the axial anomaly in two and four dimensions using
our modified path integral coupled to Pauli-Villars heavy fields. At the end,
we conclude and give some advantages of our method. The paper is followed by a
an appendix Fujikawa’s method in calculating axial anomaly in NC-QED.

##  2 Noncommutative $U(1)$ Gauge Theory

Non-commutative field theories are defined on space-time which the following
commutation relation between its coordinates is satisfied
$[\widehat{x}_{i},\widehat{x}_{j}]=i\;\theta^{ij}$, where $\theta^{ij}$ is a
$\mathbb{R}^{4}\bigotimes\mathbb{R}^{4}$ real antisymmetric constant. The
usual product is promoted to $\star$-Moyal product in the non-commutative
theories. Let $C^{\infty}(\mathbb{R}^{2n})$ be the space of smooth functions,
$f:\mathbb{R}^{2n}\rightarrow\mathbb{C}^{n}$, for $f,g$ $\epsilon$
$C^{\infty}(\mathbb{R}^{2n})$ the Moyal $\star$ product is defined as

$f\;\star\;g(x)=e^{i\frac{\theta_{\mu\nu}}{2}\;\frac{\partial}{\partial{\zeta_{\mu}}}\;\frac{\partial}{\partial{\eta_{\nu}}}}\;f(x+\zeta)\;g(x+\eta)\mid_{\zeta=\eta=0}$
(2.1)

Where $\theta_{\mu\nu}$ is a real antisymmetric constant.The previous Moyal
Star $\star$ product satisfies the $C^{\star}$-algebra. by expanding the
previous formula up to first order, we find that

$f\;\star
g=f\;q+\frac{i}{2}\;\theta^{\mu\nu}\;\partial_{\mu}f\;\partial_{\nu}g+\mathcal{O}(\theta^{2})$
(2.2)

The star product has many properties like

$f\;\star g\mid_{\theta}=g\;\star f\mid_{-\theta}$ (2.3)

Suppose $h$ is another smooth function , we have the following cyclic property

$\int_{-\infty}^{\infty}\;d^{D}x\;(f\star g\star
h)=\int_{-\infty}^{\infty}\;d^{D}x\;(h\star f\star
g)=\int_{-\infty}^{\infty}\;d^{D}x\;(g\star h\star f)$ (2.4)

The regular commutator is generalized in the non-commutative case to the Moyal
bracket

$[f,g]_{\star}=f\star\;g-g\star\;f$ (2.5)

The deformed Jacobi identity is

$[f,[g,h]_{\star}]_{\star}+[h,[f,g]_{\star}]_{\star}+[g,[h,f]_{\star}]_{\star}=0$
(2.6)

Remarkably, the commutator in Moyal product involves both ordinary commutators
and anti-commutators[8].

$f\star\;g-g\star\;f=[f,g]_{(\star,even)}+\\{f,g\\}_{(\star,odd)}$ (2.7)
$[f,g]_{(\star,even)}=[f,g]+(\frac{i}{2})^{2}\;\theta^{\mu\nu}\theta^{\rho\sigma}[\partial_{\mu}\partial_{\rho}f,\partial_{\nu}\partial_{\sigma}g]+\mathcal{O}(\theta^{4})$
(2.8)
$\\{f,g\\}_{(\star,odd)}=\frac{i}{2}\theta^{\mu\nu}\\{\partial_{\mu}f,\partial_{\nu}g\\}+(\frac{i}{2})^{3}\theta^{\mu\nu}\theta^{\rho\sigma}\theta^{\kappa\lambda}\\{\partial_{\mu}\partial_{\rho}\partial_{\kappa}f,\partial_{\nu}\partial_{\sigma}\partial_{\lambda}g\\}+\mathcal{O}(\theta^{5})$
(2.9)

Let us now define symmetric and antisymmetric star product as

$(f\star g)_{s}=\frac{1}{2}(f\star g+g\star
f)=fg+(\frac{i}{2})^{2}\;\theta^{\mu\nu}\theta^{\rho\sigma}\;\partial_{\mu}\partial_{\rho}f\;\partial_{\nu}\partial_{\sigma}g+\mathcal{O}(\theta^{4})$
(2.10)

$(f\star g)_{a}=\frac{1}{2}(f\star g-g\star
f)=(\frac{i}{2})\theta^{\mu\nu}\;\partial_{\mu}f\partial_{\nu}g+(\frac{i}{2})^{3}\theta^{\mu\nu}\theta^{\rho\sigma}\theta^{\kappa\lambda}\;\partial_{\mu}\partial_{\rho}\partial_{\kappa}f\;\partial_{\nu}\partial_{\sigma}\partial_{\lambda}g+\mathcal{O}(\theta^{5})$
(2.11)

The covariant derivative can be written as

$D_{\mu}=\partial_{\mu}+i\;g\;[A_{\mu},\;\;]_{\star}$ (2.12)

Where $A_{\mu}$ is the $U(1)$ gauge connection. The two-form field strength
tensor is given by the commutator

$[D_{\mu},D_{\nu}]_{\star}=i\;g\;F_{\mu\nu}(x)$ (2.13)
$F_{\mu\nu}(x)=\partial_{\mu}A_{\nu}(x)-\partial_{\nu}A_{\mu}(x)+ig[A_{\mu}(x),A_{\nu}(x)]_{\star}$
(2.14)

The field strength $F_{\mu\nu}$ transforms like

$F_{\mu\nu}(x)\rightarrow
F_{\mu\nu}^{{}^{\prime}}(x)=U(x)\star\;F_{\mu\nu}(x)\;\star U^{-1}(x)$ (2.15)

From the equation (2.14) we see that the two-form field strength tensor has a
similar form to field strength in Yang-Mills gauge theory, so noncommutative
QED exhibits properties like the usual Yang-Mills theories which has non-
Abelian gauge group at the commutative level.
The pure ( not-coupled to matter) non-commutative $U(1)$ gauge action in any
arbitrary dimension $D$ is

$S_{\emph{gauge}}[A_{\mu}]=-\frac{1}{4}\;\int d^{D}x\;F_{\mu\nu}(x)\star
F^{\mu\nu}(x)$ (2.16)

The non-commutative gauge action transforms as

$S_{\emph{gauge}}\rightarrow S^{{}^{\prime}}_{\emph{gauge}}=-\frac{1}{4}\;\int
d^{D}x\;U(x)\star\;F_{\mu\nu}\star\;U(x)^{-1}\star\;U(x)\star\;F^{\mu\nu}(x)\star\;U(x)^{-1}$
(2.17)

Note that the $U(1)$ Lie group is deformed in the non-commutative case. It can
be expanded as [5]

$\displaystyle
U(x)_{\star}=(e^{i\;a})_{\star}=\sum_{n=o}^{\infty}\;\frac{(ia)^{n}_{\star}}{n!}$
(2.18)

On other hand ,the unitarity condition still valid in the non-commutative as
it should be since we are seeking a unitary class of transformations

$U(x)\star\;U^{-1}(x)=U^{-1}(x)\star\;U(x)=\mathbb{I}$ (2.19)

Where $\mathbb{I}$ is the $N\;\times\;N$ identity matrix , where $N$ is the
dimension of the Lie group representation. In the fundamental representation
of $U(1)$, the identity matrix $\mathbb{I}$ is simply the number one.
The local gauge transformations of the gauge field are

$A_{\mu}\rightarrow A^{{}^{\prime}}_{\mu}(x)=U(x)\star A_{\mu}\star
U^{-1}(x)+\frac{i}{g}(\partial_{\mu}U)\star U^{-1}(x)$ (2.20)

where $U$ is the unitary transformation matrix.it has the form
$U(x)=(e^{ig\;a(x)})_{\star}$ and it can be expanded as follows

$U(x)=(e^{ig\;a(x)})_{\star}=1+ig\;a(x)-\frac{g^{2}}{2}a(x)\star
a(x)+\mathcal{O}(a^{3})$ (2.21)

The Fermionic (matter) action can be expressed as

$S_{F}[\overline{\psi},\psi]=\int
d^{D}x\;[\;i\overline{\psi}(x)\;\gamma^{\mu}\star\;D_{\mu}\;\psi(x)-m\;\overline{\psi}(x)\;\star\psi(x)\;]$
(2.22)

The covariant derivative that acts on $\psi$ is

$D_{\mu}\psi=\partial_{\mu}\psi(x)+ig\;A_{\mu}\star\psi(x)$ (2.23)

The transformation of matter fields are

$\displaystyle\psi^{{}^{\prime}}(x)=U(x)\star\psi(x),$ (2.24)
$\displaystyle\overline{\psi^{{}^{\prime}}}(x)=\overline{\psi}(x)\star
U^{-1}(x).$

Up to first order in parameter $a(x)$, the previous matter fields transform
like

$\displaystyle\psi^{{}^{\prime}}(x)\rightarrow\psi(x)+ig\;a(x)\star\psi(x),$
(2.25)
$\displaystyle\overline{\psi^{{}^{\prime}}}(x)\rightarrow\overline{\psi}(x)-ig\;\overline{\psi}(x)\star
a(x).$

In the noncommutative case, the conserved current is

$j^{\mu}=\overline{\psi}(x)\star\gamma^{\mu}\psi(x)=\overline{\psi}\gamma^{\mu}\psi+\frac{i}{2}\;\theta^{\mu\nu}\;\partial_{\mu}\overline{\psi}\gamma^{\mu}\partial_{\nu}\psi+\dots$
(2.26)

Finally, we can write the total NC-QED action after considering the
matter+gauge+ghost contributions altogether as

$\displaystyle S_{\emph{NC-QED}}[A_{\mu},\overline{\psi},\psi]=\int
d^{D}x\;-\frac{1}{4}\;F_{\mu\nu}(x)\star
F^{\mu\nu}(x)+\overline{\psi}(x)\star(i\not{D}-m)\;\psi(x)$ (2.27)
$\displaystyle+\mathcal{L}_{\emph{gauge-fixing}}+\mathcal{L}_{\emph{ghost}}$

Where the gauge-fixing action $S_{\emph{gauge-fixing}}$ and the ghost action
$S_{\emph{ghost}}$ are given respectively as

$S_{\emph{gauge-fixing}}[A_{\mu}]=-\frac{1}{2\xi}\int
d^{D}x\;(\partial_{\mu}A_{\mu})^{2},$ (2.28)

$S_{\emph{ghost}}[A_{\mu},C,\overline{C}]=\int
d^{D}x\;\partial_{\mu}\overline{C}(\partial_{\mu}C+ig[A_{\mu},C]_{\star}).$
(2.29)

Where $C$ and $\overline{C}$ are the Faddeev-Popov ghosts.

The corresponding Euler-Lagrange equations for the total NC-QED 2.27 are

$i\gamma^{\mu}\partial_{\mu}\psi-\;m\psi-g\gamma^{\mu}A_{\mu}\star\psi=0,$
(2.30)
$i\partial_{\mu}\overline{\psi}\gamma^{\mu}+m\overline{\psi}+g\overline{\psi}\gamma^{\mu}\star
A_{\mu}=0,$ (2.31)
$D_{\mu}F^{\mu\nu}+\frac{1}{\xi}\partial^{\nu}\partial_{\mu}A^{\mu}=gj^{\nu}.$
(2.32)

##  3 Pauli-Villars Regularization Method

The main idea of Pauli-Villars regularization scheme is to introduce a heavy
fields (regulators) and sending their masses to infinity at the end of
calculations. Thus integrals of the form

$\mathcal{I}=\int\frac{d^{4}k}{(2\pi)^{4}}\;\frac{1}{(k^{2}-m^{2}+i\epsilon)}$
(3.1)

will become in the Pauli-Villars regularization scheme like

$\mathcal{I}_{\emph{Pauli-Villars scheme
}}=\int\frac{d^{4}k}{(2\pi)^{4}}\big{[}\;\frac{1}{(k^{2}-m^{2}+i\epsilon)}-\frac{1}{(k^{2}-M^{2}+i\epsilon)}\big{]}$
(3.2)

Where $k$ is the momentum of the particle, $m$ its mass and $M$ is the
fictitious large mass which corresponds to Pauli-Villars field which are
usually referred to as ghost.

The calculation of chiral anomaly in the commutative gauge theories using
Pauli-Villars method was elaborated in [9, 10].

###  3.1 Two-Dimensional Anomaly

The modified NC - massless Dirac action in two-dimensions is

$S_{F}=\int
d^{2}x\;[\overline{\psi}(x)\star\;i\gamma^{\mu}D_{\mu}\;\star\psi(x)+\overline{\eta}(x)\star\;(i\gamma^{\mu}D_{\mu}-M)\;\star\eta(x)]$
(3.3)

Where $\psi(x)$, $\overline{\psi}(x)$ are the two-dimensional spinor fields
and $\eta(x)$, $\overline{\eta}(x)$ are Bose-like (unphysical ) spinor fields.
These unphysical fields are represented as complex numbers rather than
Grassmann numbers like the usual physical spinors. We adopt Pauli matrices as
basis throughout our calculations of the axial anomaly in two dimensions

$\displaystyle\sigma_{1}=\gamma_{1}=\left(\begin{array}[]{ c c }0&1\\\
1&0\end{array}\right),\;\;\;\sigma_{2}=\gamma_{2}=\left(\begin{array}[]{ c c
}0&-i\\\
i&0\end{array}\right),\;\;\;\sigma_{3}=\gamma_{3}=\left(\begin{array}[]{ c c
}1&0\\\ 0&-1\end{array}\right)$ (3.10)

$\gamma_{3}$ is the chirality operator and it is defined to be
$\gamma_{3}=-i\;\gamma_{1}\gamma_{2}$.
The trace formulas in two-dimensions are

${\rm Tr}(\gamma_{\mu}\gamma_{\nu})=2g_{\mu\nu}$ (3.11) ${\rm
Tr}(\gamma_{3}\gamma_{\mu}\gamma_{\nu})=2i\;\epsilon_{\mu\nu}$ (3.12)

The Euclidean effective action is

$e^{-W}=\int\mathcal{D}\psi\;\mathcal{D}\overline{\psi}\;\mathcal{D}\eta\;\mathcal{D}\overline{\eta}\;e^{-S_{F}},$
(3.13)

and the corresponding conserved vector current is

$J_{\mu}=\overline{\psi}(x)\star\;\gamma_{\mu}\;\psi(x)+\overline{\eta}(x)\star\;\gamma_{\mu}\;\eta(x).$
(3.14)

The associated axial-vector current is

$J^{3}_{\mu}=\overline{\psi}(x)\star\;\gamma_{\mu}\gamma^{3}\;\psi(x)+\overline{\eta}(x)\star\;\gamma_{\mu}\gamma^{3}\;\eta(x).$
(3.15)

Under chiral transformation and up to first order in $a(x)$, the matter fields
transform like

$\displaystyle\psi^{{}^{\prime}}(x)\rightarrow\psi(x)+i\;\gamma_{3}\;a(x)\;\star\psi(x),$
(3.16)
$\displaystyle\overline{\psi^{{}^{\prime}}}(x)\rightarrow\overline{\psi}(x)+i\;\overline{\psi}(x)\star\;\gamma_{3}\;a(x),$
$\displaystyle\eta^{{}^{\prime}}(x)\rightarrow\eta(x)+i\;\gamma_{3}\;a(x)\star\;\eta(x),$
$\displaystyle\overline{\eta^{{}^{\prime}}}(x)\rightarrow\overline{\eta}(x)+i\;\overline{\eta}(x)\star\;\gamma_{3}\;a(x).$

The variation of modified NC- action is

$\displaystyle\delta{S_{F}}$ $\displaystyle=-\int
d^{2}x\;\big{[}\overline{\psi}(x)\star\;\gamma^{\mu}\gamma_{3}\;\partial_{\mu}a(x)\star\;\psi(x)+$
(3.17)
$\displaystyle\overline{\eta}(x)\star\;\gamma^{\mu}\gamma_{3}\;\partial_{\mu}a(x)\star\;\eta(x)+2iM\;\overline{\eta}(x)\;\gamma_{3}\star\;a(x)\star\;\eta(x)\big{]}.$

By virtue of cyclic property associated with $\star$-product (2.4), we can
write

$\int
d^{2}x\big{[}\overline{\psi}(x)\star\;\gamma^{\mu}\gamma_{3}\;\partial_{\mu}a(x)\star\;\psi(x)+\overline{\eta}(x)\star\;\gamma^{\mu}\gamma_{3}\;\partial_{\mu}a(x)\star\;\eta(x)\big{]}=\int
d^{2}x\;a(x)\star\;\partial_{\mu}J^{\mu}_{3}(x).$ (3.18)

The axial Ward identity can be computed by imposing the following condition

$\frac{\delta{W}}{\delta{a(x)}}\big{|}_{a(x)=0}=0,$ (3.19)

Which is equivalent to

$\partial_{\mu}<J^{\mu}_{3}(x)>=\lim_{M\rightarrow\infty}2iM\;<\overline{\eta}(x)\;\gamma_{3}\star\;\eta(x)>.$
(3.20)

The averaged regularized chiral current in two dimensions is

$\displaystyle<J^{\mu}_{3}(x)>=$ (3.21)
$\displaystyle\frac{\int\mathcal{D}\psi\mathcal{D}\overline{\psi}\mathcal{D}\eta\mathcal{D}\overline{\eta}(\overline{\psi}\gamma^{\mu}\gamma_{3}\star\psi+\overline{\eta}\gamma^{\mu}\gamma_{3}\star\eta)e^{-\int
d^{2}x\overline{\psi}\star i\not{D}\star\psi-\int
d^{2}x\overline{\eta}\star(i\not{D}-M)\star\eta}}{\int\mathcal{D}\psi\mathcal{D}\overline{\psi}\mathcal{D}\eta\mathcal{D}\overline{\eta}e^{-\int
d^{2}x\overline{\psi}\star i\not{D}\star\psi-\int
d^{2}x\overline{\eta}\star(i\not{D}-M)\star\eta}},$

Where $\eta$, $\overline{\eta}$ are the regulator and antiregulator Bose-like
fields respectively.
The anomaly in two dimensions is

$<\overline{\eta}(x)\;\gamma_{3}\;\star\eta(x)>=\frac{\int\mathcal{D}\eta\;\mathcal{D}\overline{\eta}\;(\overline{\eta}(x)\;\gamma_{3}\;\star\eta(x))\;\star
e^{-\int
d^{2}x\;\overline{\eta}(i\not{D}-M)\star\;\eta}}{\int\mathcal{D}\eta\;\mathcal{D}\overline{\eta}\;e^{-\int
d^{2}x\;\overline{\eta}\star(i\not{D}-M)\star\;\eta}}$ (3.22)

The eigenvalues of Dirac operator are defined as

$\not{D}\star\phi_{n}(x)=\lambda_{n}\phi_{n}(x)$ (3.23)

Where $\phi_{n}$ are the eigenfunctions that correspond to the eigenvalue
$\lambda_{n}$. They satisfy the following orthonormal conditions

$\displaystyle\sum\phi^{\dagger}_{n}(x)\star\;\phi(y)=\delta^{4}(x-y),$ (3.24)
$\displaystyle\int
d^{4}x\;\phi^{\dagger}_{n}(x)\star\;\phi_{m}(x)=\delta_{nm}.$

We expand the Pauli-Villars fields using the eigenfunctions of Dirac operator
as a basis

$\displaystyle\eta(x)=\sum\;c_{n}\;\phi_{n}(x),$ (3.25)
$\displaystyle\overline{\eta}(x)=\sum\;\overline{c}_{n}\;\phi^{\dagger}_{n}(x).$

Where $c$ and $\overline{c}$ are ordinary complex numbers not a Grassmann
numbers since the Pauli-Villars field $\eta$ and $\overline{\eta}$ obey Bose
statistics. We insert the previous field expansion in the equation (3.22), we
get

$<\overline{\eta}(x)\;\gamma_{3}\star\;\eta(x)>=\frac{\prod_{n}\;\int
dc_{n}\;\int
d\overline{c}_{n}\;\sum_{\ell,m}\;\overline{c}_{\ell}\;\phi^{\dagger}_{\ell}(x)\;\gamma_{3}\;c_{m}\star\;\phi(x)_{m}\;e^{-\sum_{k}\;\overline{c}_{k}\;c_{k}(i\lambda_{k}-M)}}{\prod_{n}\int
dc_{n}\int
d\overline{c}_{n}\;e^{-\sum_{k}\;\overline{c}_{k}\;c_{k}\;(i\lambda_{k}-M)}}$
(3.26)

After performing the integrals, we find the simple expression for the anomaly

$<\overline{\eta}(x)\;\gamma_{3}\star\;\eta(x)>=\sum_{\ell}\;\frac{\phi^{\dagger}_{\ell}(x)\;\gamma_{3}\;\star\phi_{\ell}(x)}{i\;\lambda_{\ell}-M}.$
(3.27)

The sum can be rewritten as
$\sum_{\ell}\;\phi^{\dagger}_{\ell}(x)\;\gamma_{3}\;\star\phi_{\ell}(x)=\mathrm{Tr}{\gamma_{3}\;\delta{(0)}}$
since $\gamma_{3}$ is a tracless matrix and $\delta(0)$ is infinite. Note that
the previous sum is not well defined, so we need to regularize it in order to
get a finite expression at the end of calculations. One way to regularize it
is by using Dirac operators,

$\lim_{M\rightarrow\infty}\;2iM\;<\overline{\eta}(x)\;\gamma_{3}\star\;\eta(x)>=\lim_{M\rightarrow\infty}\;2iM\;{\rm
Tr}\big{[}\gamma_{3}\langle x|\frac{1}{i\not{D}-M}|x\rangle\big{]}$ (3.28)

After we multiply the numerator and denominator of the previous equation by
$-i\not{D}-M$, we get

$\displaystyle\lim_{M\rightarrow\infty}\;2iM\;<\overline{\eta}(x)\;\gamma_{3}\star\;\eta(x)>=$
(3.29) $\displaystyle\lim_{M\rightarrow\infty}2M\;{\rm
Tr}\big{[}\gamma_{3}\langle
x|\frac{\not{D}}{\not{D}^{2}+M^{2}}|x\rangle\big{]}-\lim_{M\rightarrow\infty}2iM^{2}\;{\rm
Tr}\big{[}\gamma_{3}\;\langle x|\frac{1}{\not{D}^{2}+M^{2}}|x\rangle\big{]}$

We perform a Dyson expansion of the term $\frac{1}{\not{D}^{2}+M^{2}}$

$\displaystyle\lim_{M\rightarrow\infty}\;2iM\;<\overline{\eta}(x)\;\gamma_{3}\star\;\eta(x)>$
(3.30) $\displaystyle=-\lim_{M\rightarrow\infty}2iM^{2}\;{\rm
Tr}\int\frac{d^{2}k}{(2\pi)^{2}}\;\frac{\big{[}\gamma_{3}\;\;g\;\sigma_{\mu\nu}\;\;F^{\mu\nu}(x)_{\star}\big{]}}{2(k^{2}+M^{2})^{2}}$
$\displaystyle=\frac{-g\epsilon^{\mu\nu}}{2\pi}\;F_{\mu\nu}(x)_{\star}$

Where we have used the trace properties of $\gamma$ matrices in two dimensions
and taking the limit when $M\rightarrow\infty$.

We have used the following identity during the previous calculations

$\displaystyle\not{D}^{2}=D_{\mu}\gamma^{\mu}\;\star\;D_{\nu}\gamma^{\nu}=\frac{1}{4}\big{(}[\gamma^{\mu},\gamma^{\nu}]+\\{\gamma^{\mu},\gamma^{\nu}\\}\big{)}\big{(}\\{D_{\mu},D_{\nu}\\}_{\star}+[D_{\mu},D_{\nu}]_{\star}\big{)}$
(3.31) $\displaystyle=D_{\mu}\star
D^{\mu}+\frac{g}{2}\;\sigma_{\mu\nu}\;F^{\mu\nu}_{\star}$

Where $\sigma_{\mu\nu}=\frac{i}{2}[\gamma_{\mu},\gamma_{\nu}]$ and
$F_{\mu\nu}(x)_{\star}$ is the two-form non-commutative field strength.

Finally, the divergence of $J^{\mu}_{3}$ is given by

$\partial_{\mu}J^{\mu}_{3}=-\frac{g\epsilon_{\mu\nu}}{2\pi}F^{\mu\nu}(x)_{\star}$
(3.32)

###  3.2 Four-Dimensional Anomaly

The modified NC -massive Dirac action in four dimensions is

$S_{F}=\int
d^{4}x\;[\overline{\psi}(x)\star\;(i\gamma^{\mu}D_{\mu}-m)\;\psi(x)+\overline{\eta}(x)\star\;(i\gamma^{\mu}D_{\mu}-M)\;\eta(x)].$
(3.33)

Spinors act in a Complex vector space $\mathbb{V}$ endowed with a certain
algebraic structure. The bases of spinor vector space satisfy the anti-
commuting relation $\\{\gamma^{\mu},\gamma^{\nu}\\}=2\;g^{\mu\nu}$, which is a
Clifford algebra $\mathcal{C}\ell(1,3)$. $\gamma$ matrices are defined as
Endomorphisms $End(\mathbb{V})$ of the complex vector space $\mathbb{V}$.
We have considered the Dirac representation for $\gamma$ matrices

$\displaystyle\gamma^{0}=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&-1&0\\\
0&0&0&-1\end{pmatrix}\quad\gamma^{1}=\begin{pmatrix}0&0&0&1\\\ 0&0&1&0\\\
0&-1&0&0\\\ -1&0&0&0\end{pmatrix}$ (3.34)
$\displaystyle\gamma^{2}=\begin{pmatrix}0&0&0&-i\\\ 0&0&i&0\\\ 0&i&0&0\\\
-i&0&0&0\end{pmatrix}\quad\gamma^{3}=\begin{pmatrix}0&0&1&0\\\ 0&0&0&-1\\\
-1&0&0&0\\\ 0&1&0&0\end{pmatrix}$

and the chirality is defined as

$\gamma_{5}=i\;\frac{\epsilon^{\mu\nu\rho\sigma}}{4!}\;\gamma_{\mu}\;\gamma_{\nu}\;\gamma_{\rho}\;\gamma_{\sigma}=\;i\gamma^{0}\;\gamma^{1}\;\gamma^{2}\;\gamma^{3},$
(3.35)

Where $\epsilon^{\mu\nu\rho\sigma}$ is totally antisymmetric tensor.

The trace formulas in four dimensions are

${\rm Tr}(\gamma_{5}\;\gamma_{\mu}\;\gamma_{\nu})=o$ (3.36) ${\rm
Tr}(\gamma_{5}\;\gamma_{\mu}\;\gamma_{\nu}\;\gamma_{\rho}\;\gamma_{\sigma})=-4i\;\epsilon_{\mu\nu\rho\sigma}$
(3.37) ${\rm Tr}(\gamma_{\mu}\gamma_{\nu})=4\;g_{\mu\nu}$ (3.38)

alongside with the fact that trace of any odd number of $\gamma$ matrices is
zero.
The Euclidean effective action is

$e^{-W}=\int\mathcal{D}\psi\;\mathcal{D}\overline{\psi}\;\mathcal{D}\eta\;\mathcal{D}\overline{\eta}\;e^{-S_{F}}$
(3.39)

and the corresponding conserved vector current and axial-vector current are
respectively

$J_{\mu}=\overline{\psi}(x)\star\;\gamma_{\mu}\;\psi(x)+\overline{\eta}(x)\star\;\gamma_{\mu}\;\eta(x)$
(3.40)

$J^{5}_{\mu}=\overline{\psi}(x)\star\;\gamma_{\mu}\gamma^{5}\;\psi(x)+\overline{\eta}(x)\star\;\gamma_{\mu}\gamma^{5}\;\eta(x)$
(3.41)

Under chiral rotation and up to first order in $a(x)$, the matter fields
transform like

$\displaystyle\psi^{{}^{\prime}}(x)\rightarrow\psi(x)+i\;\gamma_{5}\;a(x)\;\star\psi(x),$
(3.42)
$\displaystyle\overline{\psi^{{}^{\prime}}}(x)\rightarrow\overline{\psi}(x)+i\;\overline{\psi}(x)\star\;\gamma_{5}\;a(x),$
$\displaystyle\eta^{{}^{\prime}}(x)\rightarrow\eta(x)+i\;\gamma_{5}\;a(x)\star\;\eta(x),$
$\displaystyle\overline{\eta^{{}^{\prime}}}(x)\rightarrow\overline{\eta}(x)+i\;\overline{\eta}(x)\star\;\gamma_{5}\;a(x).$

The variation of modified Pauli-Villars action is

$\displaystyle\delta{S_{F}}$ $\displaystyle=-\int
d^{4}x\;\big{[}\overline{\psi}(x)\star\;\gamma^{\mu}\gamma_{5}\;\partial_{\mu}a(x)\star\;\psi(x)+2im\;\overline{\psi}(x)\;\gamma_{5}\star\;a(x)\star\;\psi(x)+$
(3.43)
$\displaystyle\overline{\eta}(x)\star\;\gamma^{\mu}\gamma_{5}\;\partial_{\mu}a(x)\star\;\eta(x)+2iM\;\overline{\eta}(x)\;\gamma_{5}\star\;a(x)\star\;\eta(x)\big{]}$

Using the cyclic property of $\star$-product (2.4)

$\int
d^{4}x\big{[}\overline{\psi}(x)\star\;\gamma^{\mu}\gamma_{5}\;\partial_{\mu}a(x)\star\;\psi(x)+\overline{\eta}(x)\star\;\gamma^{\mu}\gamma_{5}\;\partial_{\mu}a(x)\star\;\eta(x)\big{]}=\int
d^{4}x\;a(x)\star\;\partial_{\mu}J^{\mu}_{5}(x)$ (3.44)

The following condition gives the so-called axial-Ward identity

$\frac{\delta{W}}{\delta{a(x)}}\big{|}_{a(x)=0}=0$ (3.45)

The axial-Ward (AW) identity is

$\partial_{\mu}<J^{\mu}_{5}(x)>=2im\;<\overline{\psi}(x)\;\gamma_{5}\star\;\psi(x)>+\lim_{M\rightarrow\infty}2iM\;<\overline{\eta}(x)\;\gamma_{5}\star\;\eta(x)>$
(3.46)

The real field $\psi$ obeys Fermi statistics and the ghost field( regulator)
$\eta$ obeys a Bose statistics. Thus,the regulator is used only for fermionic
loops and that is what we need since our aim is to calculate the axial
anomaly.
The second term in the right-hand side of previous equation is

$<\overline{\eta}(x)\;\gamma_{5}\;\star\eta(x)>=\frac{\int\mathcal{D}\eta\;\mathcal{D}\overline{\eta}\;(\overline{\eta}(x)\;\gamma_{5}\;\star\eta(x))\;\star
e^{-\int
d^{4}x\;\overline{\eta}(i\not{D}-M)\star\;\eta}}{\int\mathcal{D}\eta\;\mathcal{D}\overline{\eta}\;e^{-\int
d^{4}x\;\overline{\eta}(i\not{D}-M)\star\;\eta}}$ (3.47)

We expand the fields in orthonormal basis like the case of two-dimensional
anomaly

$\displaystyle\eta(x)=\sum\;c_{n}\;\phi_{n}(x),$ (3.48)
$\displaystyle\overline{\eta}(x)=\sum\;\overline{c}_{n}\;\phi^{\dagger}_{n}(x).$

We insert the expansion in equation (3.47)

$\displaystyle<\overline{\eta}(x)\;\gamma_{5}\star\;\eta(x)>=$ (3.49)
$\displaystyle\frac{\prod_{n}\;\int dc_{n}\;\int
d\overline{c}_{n}\;\sum_{\ell,m}\;\overline{c}_{\ell}\;\phi^{\dagger}_{\ell}(x)\;\gamma_{5}\;c_{m}\star\;\phi(x)_{m}\;e^{-\sum_{k}\;\overline{c}_{k}\;c_{k}(i\lambda_{k}-M)}}{\prod_{n}\int
dc_{n}\int
d\overline{c}_{n}\;e^{-\sum_{k}\;\overline{c}_{k}\;c_{k}\;(i\lambda_{k}-M)}}$

After performing the integrals in previous equation, we obtain

$<\overline{\eta}(x)\;\gamma_{5}\star\;\eta(x)>=\sum_{\ell}\;\frac{\phi^{\dagger}_{\ell}(x)\;\gamma_{5}\;\star\phi_{\ell}(x)}{i\;\lambda_{\ell}-M}$
(3.50)

As we have seen in the previous section, the
$\sum_{\ell}\;\phi^{\dagger}_{\ell}(x)\;\gamma_{5}\;\star\phi_{\ell}(x)$ is
not well defined. We solve this problem in the virtue of Dirac operators

$\lim_{M\rightarrow\infty}\;2iM\;<\overline{\eta}(x)\;\gamma_{5}\star\;\eta(x)>=\lim_{M\rightarrow\infty}\;2iM\;{\rm
Tr}\big{[}\gamma_{5}\langle x|\frac{1}{i\not{D}-M}|x\rangle\big{]}$ (3.51)

We multiply the numerator and denominator of the previous equation by
$-i\not{D}-M$. This gives

$\displaystyle\lim_{M\rightarrow\infty}\;2iM\;<\overline{\eta}(x)\;\gamma_{5}\star\;\eta(x)>$
(3.52) $\displaystyle=\lim_{M\rightarrow\infty}2M\;{\rm
Tr}\big{[}\gamma_{5}\langle
x|\frac{\not{D}}{\not{D}^{2}+M^{2}}|x\rangle\big{]}-\lim_{M\rightarrow\infty}2iM^{2}\;{\rm
Tr}\big{[}\gamma_{5}\;\langle x|\frac{1}{\not{D}^{2}+M^{2}}|x\rangle\big{]}$

By performing Dyson expansion of the term $\frac{1}{\not{D}^{2}+M^{2}}$ and
due to the trace properties of $\gamma$ matrices in four dimensions and at the
limit when $M\rightarrow\infty$. The axial anomaly is

$\displaystyle\lim_{M\rightarrow\infty}\;2iM\;<\overline{\eta}(x)\;\gamma_{5}\star\;\eta(x)>$
(3.53) $\displaystyle=-\lim_{M\rightarrow\infty}2iM^{2}\;{\rm
Tr}\;\int\frac{d^{4}k}{(2\pi)^{4}}\frac{\big{[}\gamma_{5}\;\;g^{2}\sigma_{\mu\nu}\;\sigma_{\rho\sigma}\;F^{\mu\nu}(x)\star\;F^{\rho\sigma}(x)\big{]}}{4(k^{2}+M^{2})^{3}}$
$\displaystyle=\frac{-g^{2}\epsilon^{\mu\nu\rho\sigma}}{16\pi^{2}}\;{\rm
Tr}\;(F_{\mu\nu}(x)\star\;F_{\rho\sigma}(x))$

The previous result matches the result in [6].
Finally, the divergence of $J^{\mu}_{5}$ is given by

$\partial_{\mu}J^{\mu}_{5}=2im\;{\rm
Tr}[\gamma_{5}<x|\frac{1}{i\not{D}-m}|x>]-\frac{g^{2}\epsilon^{\mu\nu\rho\sigma}}{16\pi^{2}}\;{\rm
Tr}\;(F_{\mu\nu}(x)\star\;F_{\rho\sigma}(x))$ (3.54)

##  4 Conclusion

We have calculated the $U(1)$ axial anomaly in two and four dimensions using
the modified path integral formalism coupled to Pauli-Villars regulators in
the presence of noncommutative background. One advantage of this method is
that we can easily convert from the non-perturbative path integral formalism
to the perturbative axial Ward identities at each step in the derivation. This
resilience is maintained in this method while in other methods such transition
can not be established. We also expect that our developed method can be used
to calculate anomalies for any gauge group not only for $U(1)$. Also, it can
be used in calculating chiral anomalies in curved spacetimes. We hope that our
study will give evidence about the power of path integral formalism in
comparison with other perturbative methods.

##  Appendix A Fujikawa Method In Four-Dimensions

We will follow the same general procedures that is performed in [6] with
Gaussian function regularization.
The partition function in NC-QED is

$\mathcal{Z}=\int\mathcal{D}\psi\;\mathcal{D}\overline{\psi}\;e^{-iS_{F}}.$
(A.1)

The chiral transformations for the matter fields $\psi$ and $\overline{\psi}$
are

$\displaystyle\psi\rightarrow\psi^{{}^{\prime}}(x)=U^{5}(x)\star\;\psi(x),$
(A.2)
$\displaystyle\overline{\psi}(x)\rightarrow\overline{\psi^{{}^{\prime}}}(x)=\overline{\psi}(x)\star\;U^{5}(x).$

Where $U^{5}(x)=(e^{i\gamma_{5}\;a(x)})_{\star}$.
The NC-massless Dirac action is

$S_{F}[\psi,\overline{\psi}]=i\;\int
d^{4}x\;\overline{\psi}(x)\;\gamma^{\mu}\;\star\partial_{\mu}\psi(x).$ (A.3)

under chiral rotation the previous action transforms like

$S_{F}\rightarrow S_{F}^{{}^{\prime}}=S_{F}-\int
d^{4}x\;\overline{\psi}(x)\star\;\partial_{\mu}a(x)\;\gamma^{\mu}\gamma^{5}\star\psi(x)$
(A.4)

From the cyclic property of $\star$-Product (2.4), we can write

$\int
d^{4}x\;\overline{\psi}(x)\star\partial_{\mu}a(x)\;\gamma^{\mu}\gamma^{5}\;\star\psi(x)=i\;\int
d^{4}x\;a(x)\;\star\partial_{\mu}J^{\mu}_{5}(x)$ (A.5)

Where $J^{\mu}_{5}(x)$ is the axial-current. From previous property, the
massless Dirac action becomes

$S^{{}^{\prime}}_{F}=S_{F}-i\int
d^{4}x\;a(x)\star\;\partial_{\mu}J^{\mu}_{5}(x).$ (A.6)

The Partition function under Chiral rotaion reads as

$\displaystyle\mathcal{Z}^{{}^{\prime}}$
$\displaystyle=\int\mathcal{D}\psi^{{}^{\prime}}\;\mathcal{D}\overline{\psi}(x)^{{}^{\prime}}\;e^{-iS^{{}^{\prime}}_{F}}$
(A.7)
$\displaystyle=\int\mathcal{D}\psi^{{}^{\prime}}\;\mathcal{D}\overline{\psi}^{{}^{\prime}}\;e^{-iS_{F}-\int
d^{4}x\;a(x)\star\;\partial_{\mu}J^{\mu}_{5}(x)}$

Under chiral transformations, Dirac fields transforms like

$\displaystyle\overline{\psi}^{{}^{\prime}}(x)\approx\overline{\psi}(x)+i\;\gamma_{5}\;\overline{\psi}(x)\star
a(x),$ (A.8)
$\displaystyle\psi^{{}^{\prime}}(x)\approx\psi(x)+i\;\gamma_{5}\;a(x)\;\star\psi(x).$

The fermionic measure changes like

$\displaystyle\mathcal{D}\psi^{{}^{\prime}}\;\mathcal{D}\overline{\psi}^{{}^{\prime}}$
$\displaystyle=(\mathcal{J})^{-2}\;\mathcal{D}\psi\;\mathcal{D}\overline{\psi}$
(A.9) $\displaystyle=e^{i\int
d^{4}x\;a(x)\star\;\mathcal{A}(x)}\;\mathcal{D}\psi\;\mathcal{D}\overline{\psi}$

Where $\mathcal{J}$ is the Jacobian.
The sum of exponential is not well defined
$\sum_{n}^{\infty}\;\phi^{\dagger}_{n}(x)\;\gamma_{5}\star\phi(x)_{n}=\mathrm{Tr}\gamma_{5}\;\delta(0)$.
This is because $\gamma_{5}$ is a traceless matrix and $\delta(0)$ is
infinite. Thus, we need to regularize it by introducing a Gaussian cutoff

$\sum_{n}\;\phi^{\dagger}_{n}(x)\;\gamma_{5}\star\phi_{n}(x)=\lim_{M\rightarrow\infty}\sum_{n}\;\phi^{\dagger}_{n}(x)\;\star
e^{-\frac{D_{\mu}\;\gamma^{\mu}\;\star
D_{\nu}\;\gamma^{\nu}}{M^{2}}}\;\star\phi_{n}(x)$ (A.10)

Where $M$ is a large mass and Dirac delta function is defined in the same way
as in the commutative case.
It is convenient to perform Fourier transformations
$\phi_{n}(x)=\int\frac{d^{4}k}{(2\pi)^{2}}\;e^{ik.x}\;\widetilde{\phi}_{n}(k)$
of fields. The Gaussian sum can be expressed as

$\sum_{n}\;\phi^{\dagger}_{n}(x)\star\;\gamma_{5}\phi_{n}(x)=\lim_{M\rightarrow\infty}\int\frac{d^{4}k}{(2\pi)^{4}}\;\mathrm{Tr}\;e^{-ik.x}\;\gamma_{5}\;e^{-\frac{D_{\mu}\;\gamma^{\mu}\;\star
D_{\nu}\;\gamma^{\nu}}{M^{2}}}\;e^{ik.x},$ (A.11)

where the trace in previous relation is taken over the $\gamma$ matrices and
over the Lie algebra of the action.
We can do the following algebraic manipulation of Dirac operators

$\displaystyle-
D_{\mu}\gamma^{\mu}\;\star\;D_{\nu}\gamma^{\nu}=-\frac{1}{4}\big{(}[\gamma^{\mu},\gamma^{\nu}]+\\{\gamma^{\mu},\gamma^{\nu}\\}\big{)}\big{(}\\{D_{\mu},D_{\nu}\\}_{\star}+[D_{\mu},D_{\nu}]_{\star}\big{)}$
(A.12) $\displaystyle=-D_{\mu}\star
D^{\mu}-\frac{g}{2}\;\sigma_{\mu\nu}\;F^{\mu\nu}$

where the two-form field strength is defined as
$ig\;F_{\mu\nu}=[D_{\mu},D_{\nu}]_{\star}$ and
$\sigma_{\mu\nu}=\frac{i}{2}\;[\gamma_{\mu},\gamma_{\nu}]$. We insert the
previous identity inside the Gaussian sum A.11. Also, we perform a scaling of
momenta $k_{\mu}\rightarrow Mk^{{}^{\prime}}_{\mu}$, note that the
differential Dirac operators will make a plane wave shift. Finally we arrive
to

$\displaystyle\sum_{n}\;\phi^{\dagger}_{n}(x)\star\;\gamma_{5}\;\phi_{n}(x)$
(A.13)
$\displaystyle=\lim_{M\rightarrow\infty}\;\int\frac{M^{4}d^{4}k^{{}^{\prime}}}{(2\pi)^{4}}\;e^{k^{2}}\;\mathrm{Tr}\big{[}\gamma_{5}\;e^{-\frac{(D_{mu}+iM\;k^{{}^{\prime}}_{\mu})\star(D^{\mu}+iM\;k^{{}^{\prime}\mu})}{M^{2}}}\;\star
e^{-\frac{i\gamma_{\mu}\;\gamma_{\nu}F^{\mu\nu}}{2M^{2}}}\big{]}$

We expand the exponents in previous relation. By using the $\gamma$ trace
properties with $M\rightarrow\infty$ gives the singlet anomaly

$\mathcal{A}(x)=\frac{i\;g^{2}\epsilon^{\mu\nu\rho\sigma}}{16\pi^{2}}\;\mathrm{Tr}(F_{\mu\nu}(x)\star\;F_{\rho\sigma}(x)),$
(A.14)

and the divergence of chiral current is

$\partial_{\mu}J^{\mu}_{5}(x)=\frac{-g^{2}\epsilon^{\mu\nu\rho\sigma}}{16\pi^{2}}\;\mathrm{Tr}\;(F_{\mu\nu}(x)\star\;F_{\rho\sigma}(x)).$
(A.15)

## References

  * [1] A. Connes, M. R. Douglas, and A. Schwarz _Noncommutative Geometry and Matrix Theory: Compactification on Tori_ , _JHEP_ 9802:003 (1998).
  * [2] N. Seiberg and E. Witten, _Noncommutative Geometry And String Theory_ , _JHEP_ 9909:032 (1999).
  * [3] S.Minwalla, M. Van Raamsdonk, and N. Seiberg, _Noncommutative Perturbative Dynamics_ , _JHEP_ 0001:028 (2000).
  * [4] F.Ardalan, H.Arfaei and M.M.Sheikh-Jabbari, _Noncommutative geometry from strings and branes_ , _JHEP_ 9902:016 (1999).
  * [5] S. Ferrara and M.A. Lledo, _Some Aspects of Deformations of Supersymmetric Field Theories_ , _JHEP_ 0005 (2000).
  * [6] F.Ardalan and N.Sadooghi, _Axial Anomaly In Noncommutative QED on $R^{4}$_, _Int.J.Mod.Phys._ , A16 (2001).
  * [7] J. M. Gracia-Bondia and C. P. Martin , _Chiral Gauge Anomalies on Noncommutative $\mathbb{R}^{4}$_, _Phys.Lett.B_ , 479 (2000). .
  * [8] A.H.Chamseddine, _$SL(2,\mathbb{C})$ Gravity with Complex Vierbein and Its Noncommutative Extension_, _Phys.Rev. D_ 69  (2004).
  * [9] K.Fujikawa, _Generalized Pauli-Villars Regularization and the Covariant Form of Anomalies_ , _Nucl.Phys. B_ 428 (1994).
  * [10] R.Banerjee, _Anomalie and the Index Theorem In The Path Integral approach_ , _HD-THEP_ 25 Heidelberg University (1990).


# Orbits of charged particles with an azimuthal initial velocity in a dipole
magnetic field

Hanrui Pang  School of Physical Science and Technology, Southwest Jiaotong
University, Chengdu 610031, People’s Republic of China  Siming Liu
<EMAIL_ADDRESS>School of Physical Science and Technology, Southwest
Jiaotong University, Chengdu 610031, People’s Republic of China  Rong Liu
School of Physical Science and Technology, Southwest Jiaotong University,
Chengdu 610031, People’s Republic of China

###### Abstract

Nonintegrable dynamical systems have complex structures in their phase space.
Motion of a test charged particle in a dipole magnetic field can be reduced to
a 2 degree-of-freedom (2 d.o.f.) nonintegrable Hamiltonian system. We carried
out a systematic study of orbits of charged particles with an azimuthal
initial velocity in a dipole field via calculation of their Lyapunov
characteristic exponents (LCEs) and escape times for a dimensionless energy
less and greater than 1/32, respectively. Meridian plane periodic orbits
symmetric with respect to the equatorial plane are then identified. We found
that 1) symmetric periodic orbits can be classified into several classes based
on their number of crossing points on the equatorial plane; 2) the initial
conditions of these classes locate on closed loops or closed curves going
through the origin; 3) most isolated regions of stable quasi-periodic orbits
are associated asymmetric stable periodic orbits; 4) classes of asymmetric
periodic orbits either go through the origin or terminate at flat equatorial
plane orbits with the other end approaching centers of spiral structures; 5)
there are apparent self-similarities in the above features with the decrease
of energy.

Lyapunov characteristic exponents ; 4D phase space; Fractal structure; Dipole
magnetic field; Periodic orbits

††preprint: APS/123-QED

##  I  Introduction

The behavior of nonintegrable dynamical systems is complicated due to presence
of complex structures in their phase space [1]. In particular, stable
trajectories wherever exist can affect properties of the whole system
significantly [2]. Although chaotic trajectories are ergodic, they may also
show nontrivial behaviors [3, 4]. As one of the simplest nonintegrable
Hamiltonian systems, motion of charged particles in the dipole magnetic field
has been studied extensively during the past century [5, 6, 7]. In particular,
significant efforts have been put to study periodic orbits in the Meridian
plane [8, 9, 10, 11, 12], and it has been known that stable periodic orbits
can affect statistical properties of the whole system significantly [4, 2].
Stable periodic orbits also bear significance in physical systems since their
associated stable quasi-periodic orbits can be realized and have distinct
characteristics [13, 14].

With the aid of Lyapunov characteristic exponents (LCE), stable orbits can be
readily identified since the value of their maximum LCE (mLCE) is zero.
However, the distribution of stable orbits in the phase space is complex,
showing fractal structures [15, 16, 7]. To understand this distribution, one
may identify classes of periodic orbits [12, 17]. In this paper, we focus on
studying orbits with an azimuthal initial velocity, i.e., initially static in
the Meridian plane. We first show that periodic orbits in the Meridian plane
have two static points half-a-period apart and those symmetric with respect to
the equatorial plane cross it perpendicularly twice in each period at the same
distance from the origin. Classes of such symmetric open periodic orbits are
then identified and used to explore distribution of stable orbits in the phase
space.

In §II we present the Hamiltonian equations for the Meridian plane motion of a
charged particle moving in a dipole magnetic field and prove general
properties of periodic orbits with open ends in §III. Our main results are
shown in §IV. Conclusions are drawn in §V.

##  II Motion of charged particles in the Meridian plane of a dipole magnetic
field

Since the dipole magnetic field is axisymmetric, the azimuthal canonical
momentum $p_{\phi}$ and the energy $H$ of a charged particle moving in it are
conserved if one ignores the radiation reaction force, which is a quantum
phenomenon and cannot be treated self-consistently in classical mechanism [18,
19]. Following the standard procedure [6, 7], one can introduce a
characteristic length with $p_{\phi}$, the charge of the particle $e$, and
magnetic moment of the dipole $M$, and obtain the dimensionless Hamiltonian
for motion in the Meridian plane $[\rho,z]$:

$H=\frac{1}{2}\left({p_{z}}^{2}+{p_{\rho}}^{2}\right)+V,$ (1)

where $z$ is the symmetric axis, $\rho$ is the radius in the equatorial plane,
and the effective potential

$V=\frac{1}{2}\left[\frac{1}{\rho}-\frac{\rho}{(z^{2}+\rho^{2})\,^{\frac{3}{2}}}\right]^{2}$
(2)

corresponds to the kinetic energy associated with the azimuthal motion,
$p_{\rho}$, and $p_{z}$ are the momentum in the $\rho$, and $z$ directions,
respectively. We then have the Hamiltonian differential equations of motion:

$\displaystyle\dot{z}$ $\displaystyle=\frac{\partial H}{\partial
p_{z}}=p_{z}\,,$ (3) $\displaystyle\dot{\rho}$ $\displaystyle=\frac{\partial
H}{\partial p_{\rho}}=p_{\rho}\,,$ (4) $\displaystyle\dot{p_{z}}$
$\displaystyle=-\frac{\partial H}{\partial
z}=-\frac{3z\rho\left[\frac{1}{\rho}-\frac{\rho}{(z^{2}+\rho^{2})^{\frac{3}{2}}}\right]}{(z^{2}+\rho^{2})^{\frac{5}{2}}}\,,$
(5) $\displaystyle\dot{p_{\rho}}$ $\displaystyle=-\frac{\partial
H}{\partial\rho}\ ,$ (6)
$\displaystyle=\left[\frac{1}{\rho^{2}}-\frac{3\rho^{2}}{(z^{2}+\rho^{2})^{\frac{5}{2}}}+\frac{1}{(z^{2}+\rho^{2})^{\frac{3}{2}}}\right]\left[\frac{1}{\rho}-\frac{\rho}{(z^{2}+\rho^{2})^{\frac{3}{2}}}\right],$

where an upper dot $\dot{}$ indicates a derivative with respect to the time
$t$.

##  III General Characteristics of Open Periodic Orbits

Periodic orbits are the simplest orbits and are ubiquitous in the phase space
[20]. Stable quasi-periodic orbits are always associated with stable periodic
orbits [17, 3]. It is therefore essential to understand properties of periodic
orbits to uncover the phase space structure of a dynamical system.

Particles with an azimuthal initial velocity have a zero initial velocity with
$p_{z,0}=0,p_{\rho,0}=0$ in the $[z,\rho]$ plane. Periodic orbits in this
Meridian plane have open paths and some general characteristics that are
essential to search for such orbits numerically. To prove these
characteristics, we define $\ddot{z}\equiv f(z,\rho)$ and $\ddot{\rho}\equiv
g(z,\rho)$, and the minimum positive period $T$. Static orbits with $T=0$ are
trivial and will not be considered here. The Hamiltonian equations then lead
to

$\begin{split}f(-z,\rho)&=-f(z,\rho)\,,\\\ g(-z,\rho)&=g(z,\rho)\,,\\\
\end{split}$ (7)

and

$\begin{split}z(t+T)=z(t),\ \ \dot{z}(t+T)=\dot{z}(t)\,,\\\
\rho(t+T)=\rho(t),\ \ \dot{\rho}(t+T)=\dot{\rho}(t)\,.\\\ \end{split}$ (8)

Since open orbits start with zero initial velocities, it is evident that
$z(t)=z(-t),\ \rho(t)=\rho(-t)$. Let $t=-t^{*}$, we have

$\begin{cases}\frac{d^{2}z}{{dt^{*}}^{2}}=f(z,\rho)\\\
\frac{d^{2}\rho}{{dt^{*}}^{2}}=g(z,\rho)\end{cases}{\rm with}\ \
\begin{cases}z|_{t^{*}=0}=z_{0},\ \ \frac{dz}{dt^{*}}|_{t^{*}=0}=0\,,\\\
\rho|_{t^{*}=0}=\rho_{0},\ \ \frac{d\rho}{dt^{*}}|_{t^{*}=0}=0\,,\end{cases}$
(9)

which are the same equations followed by $z(t)$ and $\rho(t)$. Then we have
$z(t)=z(t^{*}),\rho(t)=\rho(t^{*})$, i.e.,

$\begin{split}z(t)=z(-t)\,,\\\ \rho(t)=\rho(-t)\,.\end{split}$ (10)

###  III.1 Open Periodic Orbits

For open periodic orbits, it can be shown that there are two static points in
each period at times $t=0$ and $t=\frac{T}{2}$, respectively. From equation
(10), we have

$\begin{split}\dot{z}(t)=-\dot{z}(-t)\,,\\\
\dot{\rho}(t)=-\dot{\rho}(-t)\,.\end{split}$ (11)

Let $t=\frac{T}{2}$, equation (8) then gives

$\begin{split}\dot{z}(\frac{T}{2})=-\dot{z}(-\frac{T}{2})=-\dot{z}(\frac{T}{2})\,,\\\
\dot{\rho}(\frac{T}{2})=\dot{\rho}(-\frac{T}{2})=-\dot{\rho}(\frac{T}{2})\,.\end{split}$
(12)

Therefore

$\begin{split}\dot{z}(\frac{T}{2})=0\,,\\\
\dot{\rho}(\frac{T}{2})=0\,,\end{split}$ (13)

and we define $z_{1}\equiv z(\frac{T}{2}),\ \rho_{1}\equiv\rho(\frac{T}{2})$.
Clearly, $(z_{1},\rho_{1})\neq(z_{0},\rho_{0})$. Otherwise, $T$ will no longer
be the minimum positive period since $T/2<T$ is also a positive period.

Furthermore, we can prove that $t=\frac{T}{2}$ is the only $t$ that satisfies
$\dot{z}(t)=0,\ \dot{\rho}(t)=0$ in $(0,T)$. Suppose that there exists a
$t_{1}\in(0,T),t_{1}\neq\frac{T}{2}$ such that $\dot{z}(t_{1})=0,\
\dot{\rho}(t_{1})=0$, we have

$\begin{cases}\frac{d^{2}z}{{dt}^{2}}=f(z,\rho)\\\
\frac{d^{2}\rho}{{dt}^{2}}=g(z,\rho)\end{cases}{\rm with}\ \
\begin{cases}z|_{t=t_{1}}=z(t_{1}),\ \ \frac{dz}{dt}|_{t=t_{1}}=0\,,\\\
\rho|_{t=t_{1}}=\rho(t_{1}),\ \ \frac{d\rho}{dt}|_{t=t_{1}}=0\,.\end{cases}$
(14)

Let $t=-t^{*}+2t_{1}$, we have

$\begin{cases}\frac{d^{2}z}{{dt^{*}}^{2}}=f(z,\rho)\\\
\frac{d^{2}\rho}{{dt^{*}}^{2}}=g(z,\rho)\end{cases}{\rm with}\ \
\begin{cases}z|_{t^{*}=t_{1}}=z(t_{1}),\ \
\frac{dz}{dt^{*}}|_{t^{*}=t_{1}}=0\,,\\\ \rho|_{t^{*}=t_{1}}=\rho(t_{1}),\ \
\frac{d\rho}{dt^{*}}|_{t^{*}=t_{1}}=0\,.\end{cases}$ (15)

Equations (14) and (15) have the same solutions: $z(t)=z(t^{*}),\
\rho(t)=\rho(t^{*})$. Therefore

$\begin{split}z(t)=z(-t+2t_{1})=z(t-2t_{1})=z(t+T-2t_{1})\,,\\\
\rho(t)=\rho(-t+2t_{1})=\rho(t-2t_{1})=\rho(t+T-2t_{1})\,.\end{split}$ (16)

Then $0<|T-2t_{1}|\ {\rm or}\ 2t_{1}<T$ is a positive period.

###  III.2 Symmetric Open Periodic Orbits

For open periodic orbits symmetric with respect to the equatorial plane,
$(z_{0},\rho_{0})=(-z_{1},\rho_{1})$. When $t=\frac{T}{4}$, it can be shown
that they pass through the equatorial plane perpendicularly with $p_{\rho}=0$
and $z=0$ and there is only one such perpendicular crossing within
$\frac{T}{2}$. Let $z=-z^{*},\ t=t^{*}+\frac{T}{2}$ and considering equation
(7), we have

$\begin{cases}\frac{d^{2}z^{*}}{{dt^{*}}^{2}}=f(z^{*},\rho)\\\
\frac{d^{2}\rho}{{dt^{*}}^{2}}=g(z^{*},\rho)\end{cases}{\rm with}\ \
\begin{cases}z^{*}|_{t^{*}=0}=-z_{1}=z_{0},\ \
\frac{dz^{*}}{dt^{*}}|_{t^{*}=0}=0\,,\\\ \rho|_{t^{*}=0}=\rho_{1}=\rho_{0},\ \
\frac{d\rho}{dt^{*}}|_{t^{*}=0}=0\,,\end{cases}$ (17)

which are the same equations followed by $(z(t),\ \rho(t))$. Therefore
$z(t)=z^{*}(t^{*}),\ \rho(t)=\rho(t^{*})$, and we have

$\begin{split}z(t)=-z(t-\frac{T}{2})\,,\\\
\rho(t)=\rho(t-\frac{T}{2})\,.\end{split}$ (18)

Taking the time derivative of $\rho(t)$, one has

$\begin{split}\dot{\rho}(t)=\dot{\rho}(t-\frac{T}{2})\,.\end{split}$ (19)

Let $t=\frac{T}{4}$, equations (10), (11), (18), and (19) give

$\begin{split}z(\frac{T}{4})=z(-\frac{T}{4})=-z(-\frac{3T}{4})=-z(\frac{T}{4})\,,\\\
\dot{\rho}(\frac{T}{4})=\dot{\rho}(-\frac{T}{4})=-\dot{\rho}(\frac{T}{4})\,.\end{split}$
(20)

Therefore

$\begin{split}z(\frac{T}{4})=0\,,\\\ \dot{\rho}(\frac{T}{4})=0\,.\end{split}$
(21)

One can prove that $t=\frac{T}{4}$ is the only $t\in[0,\frac{T}{2})$ that
satisfies $z(t)=0,\ \dot{\rho}(t)=0$. Suppose there exists a
$t_{1}\in[0,\frac{T}{2})$, and $t_{1}\neq\frac{T}{4}$ such that $z(t_{1})=0,\
\dot{\rho}(t_{1})=0$, we have

$\begin{cases}\frac{d^{2}z}{{dt}^{2}}=f(z,\rho)\\\
\frac{d^{2}\rho}{{dt}^{2}}=g(z,\rho)\end{cases}{\rm with}\ \
\begin{cases}z|_{t=t_{1}}=0,\ \ \frac{dz}{dt}|_{t=t_{1}}=v_{z}(t_{1})\,,\\\
\rho|_{t=t_{1}}=\rho_{t_{1}},\ \ \frac{d\rho}{dt}|_{t=t_{1}}=0\,.\end{cases}$
(22)

Let $t=-t^{*}+2t_{1},\ z=-z^{*}$, equation (7) implies:

$\begin{cases}\frac{d^{2}z^{*}}{{dt^{*}}^{2}}=f(z^{*},\rho)\\\
\frac{d^{2}\rho}{{dt^{*}}^{2}}=g(z^{*},\rho)\end{cases}{\rm with}\ \
\begin{cases}z^{*}|_{t^{*}=t_{1}}=0,\ \
\frac{dz^{*}}{dt^{*}}|_{t^{*}=t_{1}}=v_{z}(t_{1})\,,\\\
\rho|_{t^{*}=t_{1}}=\rho_{t_{1}},\ \
\frac{d\rho}{dt^{*}}|_{t^{*}=t_{1}}=0\,.\end{cases}$ (23)

Equations (22) and (23) have the same solutions: $z(t)=z^{*}(t^{*}),\
\rho(t)=\rho(t^{*})$. Equations (10), (18), and (19) then lead to

$\begin{split}z(t)=-z(-t+2t_{1})=-z(t-2t_{1})=z(t+\frac{T}{2}-2t_{1})\,,\\\
\rho(t)=\rho(-t+2t_{1})=\rho(t-2t_{1})=\rho(t+\frac{T}{2}-2t_{1})\,.\end{split}$
(24)

Then $0<|\frac{T}{2}-2t_{1}|\leq\frac{T}{2}<T$ is a positive period.

One can actually show that any open orbit that crosses the equatorial plane
perpendicularly must be symmetric open periodic orbits. Equations (10), (22),
and (23) imply

$\begin{split}z(t)=-z(-t+2t_{1})=-z(t-2t_{1})\,,\\\
\rho(t)=\rho(-t+2t_{1})=\rho(t-2t_{1})\,.\end{split}$ (25)

Then we have

$\begin{split}z(t)=z(t-4t_{1})\,,\\\ \rho(t)=\rho(t-4t_{1})\,.\end{split}$
(26)

Therefore $T=4t_{1}$ is a positive period and $z(0)=-z(T/2)$,
$\rho(0)=\rho(T/2)$.

##  IV Results

Figure 1:  The distribution of the mLCE and the escape time of particles with
an azimuthal initial velocity. The cyan lines indicate the maximum energy of
periodic orbits [11]. The radial acceleration is positive beyond the outer
boundary. Besides several blocks of stable quasi-periodic orbits associated
with stable periodic orbits in the equatorial plane, there are many small
isolated regions of stable quasi-periodic orbits. See text for details.

Using the same numerical codes used by Liu  _et al._ [7], i.e. RK45 formula
being used to solve the relevant ordinary differential equations with default
absolute and relative error tolerance of $10^{-9}$, we have the main result
shown in Figure 1. For $H<1/32$, the $(\rho,\ z)$ plane is scanned to obtain
the mLCE of each orbit starting with $\rho=\rho_{0},\ z=z_{0},\
p_{\rho}=p_{z}=0$. It should be pointed out that at a given energy, most
chaotic orbits should have the same value of the mLCE [16]. Finite integration
time and stickiness of orbits near stable orbits make the numerical results
not as uniform as expected [21, 22, 23].

Since $\dot{r}>0$ with $r=(\rho^{2}+z^{2})^{1/2}$ for $r>(4H)^{-1/3}$ ,[11] we
only consider the domain with $r\leq(4H)^{-1/3}$, beyond which all particles
escape to infinity. In the Meridian plane, this upper limit is equivalent to
the magnetic field line passing through the point ($z=0,\ \rho=2$), then we
have $r\leq 2\cos^{2}\lambda$, where $\cos\lambda=\rho/r$. For $H\geq 1/32$,
we calculate the escape time $t_{\rm esc}$ that it takes for the particle to
reach $\rho=2$. Particles with $t_{\rm esc}>10^{4}$ are considered as trapped
particles, whose orbits are likely stable. It is interesting to note that
there are spiral structures in the escape time and most isolated regions of
stable orbits locate at some turning points of spiral structures. Such
structures also appear to emerge repeatedly approaching the origin. We also
indicate the contours for the maximum energy $H\simeq 0.80819$ that periodic
orbits can have [11] and $H=1/32$ in Figure 1.

The distribution of stable orbits (dark red) shows complex structures [24].
Besides stable orbits at low energies and those associated stable orbits in
the equatorial plane [8, 12], there are many small islands of stable orbits
spread over the whole phase space. The boundaries of regions of stable orbits
appear to be fractal, and there are apparent self-similarities with the
decrease of energy. We will try to understand these structures by identifying
classes of symmetric open periodic orbits and analyzing their relations to
regions of stable orbits.

Figure 2:  Same as Fig. 1 with the open-path periodic orbits s1 indicated by
gray lines. From left and right to middle, we have the families $f_{2}$,
$f_{6}$,… and $f_{0}$, $f_{4}$, … defined by Markellos  _et al._ [11],
respectively. Regions of stable orbits at the low energy end of each family
have similar structures.

###  IV.1 Classes of Symmetric Open Periodic Orbits

Figure 3:  Same as Fig. 1 with the open-path periodic orbits s2 and s1
indicated by red and grey lines, respectively.

Figure 4:  Upper: Same as Fig. 1 with the open-path periodic orbits s3
indicated by red lines. Lower: Same as the upper panel except that the color
for $H<1/32$ indicates the crossing time $t_{\rm cross}$ defined in the text.

s1: Given the general characteristics of symmetric open periodic orbits
discussed in the previous section, one can readily identify classes of such
orbits numerically. We first search for the principle families introduced by
Markellos  _et al._ [11], members of which cross the equatorial plane once in
each half a period. We search for such orbits by following the radial velocity
$p_{\rho}$ at the first crossing of the equatorial plane of each orbit.
Principle symmetric open-path periodic orbits have a null radial velocity at
the first crossing of the equatorial plane.

Figure 2 shows families of such orbits (grey lines), which we labelled as
"s1". They corresponds to families identified by Markellos  _et al._ [11] with
even indexes. However, in contrast to these families in the $[\rho,\ p_{z}]$
plane, each family locates on one side of the "thalweg", i.e., the bottom of
the potential $V$. From high to low energies, we have $f_{0}$, $f_{4}$, … and
$f_{2}$, $f_{6}$, .. on the right and the left side of the "thalweg",
respectively. The bulk of stable orbits associated with the low energy part of
each family have similar features at the boundary, and $f_{0}$ goes through a
prominent region of stable orbits with $H>1/32$, which agrees to the results
of Markellos  _et al._ [12]. Most stable orbits associated with the high-
energy stable periodic orbits in the equatorial plane [8] are not associated
with families in s1, except for $f_{2}$, both of whose ends on the equatorial
plane are stable.

s2: With the principle families identified above, one can readily identify the
second class of symmetric open period orbits, which cross the equatorial plane
perpendicularly at the second crossing from the initial location. Numerically,
when scanning the $(\rho_{0},z_{0})$ plane, we trace the radial velocity
$p_{\rho}$ at the second crossing of the equatorial plane of each orbit.
Periodic orbits in s2 have $p_{\rho}=0$ at the second crossing of the
equatorial plane. The orbits obtained this way include those in s1, which have
been identified above and can be readily separated from those in s2.

Figure 3 shows s2 (red) together with s1 (grey). For each family in s1, there
are 2 pairs of families in s2, with two families on each side of the thalweg,
that go through the origin and terminate at 2 pairs of stable orbits in the
equatorial plane at the other ends. Each pair of these stable orbits actually
correspond to the two ends of one orbit in the equatorial plane. One of these
families associated with $f_{0}$ in s2 crosses $f_{0}$ in the highest energy
region of stable orbits. Another crosses $f_{0}$ at a lower energy. The orbit
in s2 reduces to s1 at the crossing point of these two families. There are
fine structures at the crossings between families in s1 and s2. The other two
families on the other side of the thalweg do not go through families in s1.

The other families appear as loops and locate between these family pairs, and
the number of loops increases quickly with the decrease of energy. These loops
appear to trace spiral structures in the escape time and do not cross families
in s1. Moreover, families in one class do not go through each other.

s3: Similarly one can identify the third class of symmetric open periodic
orbits that cross the equatorial plane perpendicularly at the third crossing
from each of their initial location. Families in this class are shown in
Figure 4. Similar to the class s2, for each family in s1, there are 4 pairs of
families in s3 that go through the origin and terminate at 4 pairs of stable
periodic orbits in the equatorial plane at the other ends. Two of these
families go through $f_{0}$ in the highest energy region of stable orbits,
while the other two cross $f_{0}$ at lower energies. Fractal structures are
seen near these crossing points. The other four families at the inner side of
the potential do not cross families in s1.

Compared to s2, there are many more loops between these families, that clearly
trace spiral structures in the escape time and approach to their centers. To
reveal similar spiral structures for $H<1/32$, one can introduce $\rho_{\rm
max}=[1-(1-4(2H)^{1/2})^{1/2}]/[2(2H)^{1/2}]$ and $\rho_{\rm
min}=[(1+4(2H)^{1/2})^{1/2}-1]/[2(2H)^{1/2}]$ that correspond to the maximum
and the minimum distance to the origin for motion in the equatorial plane with
an energy $H$, respectively. We define the crossing time $t_{\rm cross}$ as
the first crossing of $\rho=\rho_{\rm max}-(\rho_{\rm max}-\rho_{\rm min})/50$
with a positive $p_{\rho}$. This crossing time is indicated in the lower panel
of Figure 4. Similar to the escape time, initial conditions with the crossing
time greater than $10^{4}$ are indicated with dark red color and are well
correlated with stable orbits in the upper panel of Figure 4. There are spiral
structures in the crossing time that are associated with family loops in s3.
Similar to s2, the family loops do not appear to cross any other families, and
families in s2 do not appear to cross families in s3 either, which is
different from those in the [$\rho,\ p_{z}$] plane, where sub-branches of
periodic orbits may cross each other [12].

###  IV.2 Classes of Asymmetric Open Periodic Orbits

Figure 5:  Examples of stable periodic orbits (thin solid lines) with the
corresponding energy contours (thick solid lines) and the bottom of the
potential $V$ (dotted) also shown.

Figure 6:  More examples of stable periodic orbits.

Although the s3 class of symmetric open periodic orbits go through the most
prominent regions of stable orbits, the chain of isolated regions of stable
orbits at the turns of the main spiral structure are clearly not associated
with s3 (See Fig. 1). We therefore randomly plot a periodic orbit from each of
these regions to reveal their properties. Figure 5 shows some of these orbits.
The upper panel gives examples of symmetric periodic orbits in the highest
energy region of stable orbits. Since all classes of symmetric periodic orbits
go through this region, one can readily find orbits for each of these three
classes presented above. Besides these trajectories, we also plot the
corresponding energy contours to better demonstrate the orbital behavior.

The middle and lower panels show the upper and lower chains of islands of
stable orbits. Based on the number of crossing through the equatorial plane
and through the bottom of the potential $V$: $r=\cos^{2}\lambda$ in each half
a period, we label these orbits with the letter "a" (for asymmetric) with the
lower and upper indexes indicating these numbers, respectively. It is evident
that orbits with an odd number of the upper index have a static point on each
side of the thalweg while those with an even number of the upper index have
their static points on one side of the thalweg. These two kinds of asymmetric
orbits belong to different families. The two orbits $a_{1}^{3}$ in the lower
panel show that for stable orbits near $H=1/32$, orbits with $H<1/32$ are very
similar to those with $H>1/32$. The critical energy $H=1/32$ is not an
important parameter determining properties of these orbits.

On the other hand, along families of symmetric open periodic orbits, the
number of crossing of the bottom of the potential may change [11]. Moreover
Markellos and Klimopoulos [10] discovered a principle asymmetric class
(surface) in the 3 dimensional sub-phase-space of $z=0$ that contains both
open and close types of periodic orbits, and they both can oscillate around
the bottom of the potential $V$ many times. They also found more than 20
regimes of stable orbits in this asymmetric class. In the 4 dimensional phase
space, the trajectories of this asymmetric class occupy a 3 dimensional volume
and its intersections with the plane of static points for the Meridian motion
should be 1 dimensional. The stable asymmetric orbits in Figure 5 all belong
to the principle asymmetric class identified by Markellos and Klimopoulos
[10]!

Nevertheless, in the plane of static points, they are divided into two
families, one with their two ends on one side of the thalweg and the other
with their two static ends on each side of the thalweg. Figure 5 shows that
these two families follow the main spiral structure and the number of
oscillations around the bottom of the potential increases as these families
approach the center of this spiral structure. Since these two families cross
the equatorial plane once in each half a period, their two ends locate on two
sides of the equatorial plane.

The principle asymmetric class terminates at the equatorial plane, the family
with an odd upper index should terminate at the equatorial plane as well. The
other ends of this family locate on the inner side of the thawhal (below the
equatorial plane). As this family approaches the center of the spiral
structure, its the other ends approach the origin. This family therefore has
another branch on the inner side of the potential. The family with an even
upper indexes crosses $f_{0}$ in the low energy end and extends towards the
origin. Since both ends of this family are on the outer side of the thalweg,
there is a one to one map for points on both sides of the crossing point with
$f_{0}$. As this family approaches the center of the main spiral structure, it
approaches the origin at the other end (below the equatorial plane).

Figure 7:  Upper: The maximum (left: $\lambda_{\rm max}$) and minimum (right:
$\lambda_{\rm min}$) angles to the equatorial plane for $H<1/32$. Lower: the
sum of the maximum and the minimum of $\lambda$.

Besides the main spiral structure, there are many other smaller spiral
structures. The upper left panel of Figure 6 shows some stable asymmetric
periodic orbits around a spiral structure below the main spiral structure and
closer to the equatorial plane (See Fig. 4). Compared with the principle
asymmetric class discussed above, it crosses the equatorial plane twice in
half-a-period. Otherwise their behaviors are very similar to those around the
main spiral structure. There therefore should be two families of asymmetric
orbits around each spiral structure, and one of them crosses $f_{0}$
approaching the origin and the other terminates at the equatorial plane. Since
both of their open ends are on one side of the equatorial plane, as the family
with an even upper index crosses $f_{0}$, its two ends merge into one and
another static point appears near one quarter of its period (i.e. half a
period of the corresponding symmetric orbit). In contrast to the principle
asymmetric class that bifurcates from $f_{0}$ with the same period, the period
doubles as this class bifurcates from $f_{0}$.

###  IV.3 Other Regions of Stable Orbits

We also notice that there are other tiny islands of stable orbits. In the
upper right panel of Figure 6, the upper orbit is similar to the families in
the upper left panel but belongs to the main spiral structure since it locates
above a family in s3 that surrounds the main spiral structure. The lower orbit
belongs to s3. It likely belongs to a branch bifurcating from the family of
the principle asymmetric class that approaches the origin. The period also
doubles at the point of bifurcation.

The lower left panel of Figure 6 shows two stable orbits that belong to
families in s2 that approach the origin. The lower one actually belongs to a
family that is associated with $f_{2}$. The top one in the lower right panel
of Figure 6 shows a stable asymmetric open orbit that appears to be
bifurcating from the family in the principle asymmetric class that terminates
at the equatorial plane. However, the top, middle, and bottom panel are
associated with $f_{0}$, $f_{2}$, and $f_{4}$, respectively, confirming self-
similarity among different families in s1.

###  IV.4 Asymmetry of Open Orbits

Calculation of the mLCE is time consuming. There is a more efficient way to
obtain region of stable asymmetric orbits. Given the asymmetry of the
corresponding periodic orbits, one may just obtain the maximum and minimum
values of $z/r$ for a given long enough computation time. Figure 7 shows the
maximum and minimum of $\lambda=\arcsin{z/r}$ and the sum of these two for a
dimensionless computation time of $10^{4}$. The most prominent isolated region
of stable asymmetric orbits near $(z=0.205\,\rho=1.62)$ can be readily seen.
Moreover, we find some asymmetric stable orbits near the boundary of the
region for the bulk of low energy stable orbits. It is interesting to note
that the center of the spiral structure below the main spiral and closer to
the equator is also highly asymmetric.

##  V Conclusions and Discussions

In this paper, we systematically study orbits of charged particles in a dipole
magnetic field with an azimuthal initial velocity that correspond to the open-
path type orbits in the Meridian plane starting with a zero velocity. To
understand the distribution of stable orbits, we identified three classes of
open periodic orbits symmetric with respect to the equatorial plane. We found
that orbits in the class s1 that cross the equatorial plane once in half a
period locate on several separated families and there are compelling evidence
of similarities among these families. We therefore focus on exploring stable
orbits associated with the highest energy one $f_{0}$.

We found that families in s2 and s3 either approach the origin with the other
end terminating at the equatorial plane or have a loop shape. The former
appears in pairs and one branch of these pairs cross $f_{0}$ at a stable
orbit. These families separate the phase space into several regions with
spiral structures in the escape and crossing times. The loops trace these
spiral structures. Families in s2 and s3 do not seem to cross each other.

Although we did not identify classes of asymmetric open periodic orbits
explicitly, by plotting orbits in isolated regions of stable orbits, we found
that except for the highest energy one and those associated with stable orbits
in the equatorial plane, the most prominent isolated regions of stable orbits
are associated with a class of asymmetric open periodic orbits, that is the
principle asymmetric family identified by [10]. In the plane with zero
meridian velocities, this class can be divided into two families: one crossing
$f_{0}$ and approaching the origin with the two open ends of each orbit
locating on either side of $f_{0}$ and the other one terminating at the
equatorial plane with the other open end of each orbit in a branch on the
other side of the bottom of the potential $V$. The other ends of these two
families trace the main spiral structure and approach to its center. The
numbers of oscillation around the bottom of the potential $V$ increase as one
moves along these families toward the center of the main spiral structure with
the other open end of the corresponding orbit approaching the origin.

Loops in the symmetric open periodic orbit classes also approach the center of
the main spiral structure with the increase in number of crossing of the
equatorial plane in half a period. The center of the main spiral structure
appears to singular in the sense that both symmetric periodic orbits
oscillating around the equatorial plane and asymmetric periodic orbits
oscillating around the bottom of the potential with their two open ends
locating in regions near the origin and this central point approach to this
same point with the increase in the number of oscillations. The nature of the
orbit at the center of spiral structure is ambiguous. One can imagine that
similar families of asymmetric orbits exist around each spiral structures in
the escape and crossing times. There are therefore many singular centers.
Further explorations are warranted.

We found evidence that, besides the principle asymmetric family bifurcating
from $f_{0}$ [12] and tracing the main spiral structure, there are many other
asymmetric families bifurcating from $f_{0}$ and tracing other spiral
structures. There are also branches bifurcating from the asymmetric families
tracing these spiral structures. Such kinds of bifurcation exist in all
families in s1.

###### Acknowledgements.

We appreciate helpful discussions with Ms. Yilin Yang. This work is partially
supported by the National Key R&D program of China under the grant Nos.
G2021166002L, and 2018YFA0404203, NSFC grants U1931204, 12147208, 11947404 and
11761131007, DFG Sino-German Collaboration Project Nos. BU 777/15-1 and MU
4255/1-1, Department of Science and Technology of Sichuan Province
No.2020YFSY0016, and by the SRTP program of the school of physical science and
technology, Southwest Jiaotong University No. 202110613088.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Chirikov [1979] B. V. Chirikov, A universal instability of many-dimensional oscillator systems, Phy. Report 52, 263 (1979).
  * Harsoula and Contopoulos [2018] M. Harsoula and G. Contopoulos, Global and local diffusion in the standard map, Phys. Rev. E 97, 022215 (2018).
  * Lichtenberg and Lieberman [1992] A. Lichtenberg and M. Lieberman,  _Regular and Chaotic Dynamics_ (Springer-Verlag, 1992).
  * Zaslavsky [2002] G. M. Zaslavsky, Chaos, fractional kinetics, and anomalous transport, Phy. Report 371, 461 (2002).
  * Störmer [1930] C. Störmer, Periodische Elektronenbahnen im Felde eines Elementarmagneten und ihre Anwendung auf Brüches Modellversuche und auf Eschenhagens Elementarwellen des Erdmagnetismus. Mit 32 Abbildungen., Zeitschrift für Astrophysik 1, 237 (1930).
  * Dragt [1965] A. Dragt, Trapped orbits in a magnetic dipole field, Reviews of Geophysics 3 (1965).
  * Liu  _et al._ [2022] R. Liu, S. Liu, F. Zhu, Q. Chen, Y. He, and C. Cai, Orbits of charged particles trapped in a dipole magnetic field, Chaos 32, 043104 (2022).
  * De Vogelaere [1950] R. De Vogelaere, Equation de hill et probleme de störmer, Canadian Journal of Mathematics 2, 440–456 (1950).
  * Markellos and Halioulias [1977] V. V. Markellos and A. A. Halioulias, On the totality of periodic motions in the meridian plane of a magnetic dipole, Astrophysics and Space Science 51, 177 (1977).
  * Markellos and Klimopoulos [1977] V. V. Markellos and S. Klimopoulos, A new family of periodic oscillations in the Störmer problem: The principal asymmetric, Astrophysics and Space Science 48, 471 (1977).
  * Markellos  _et al._ [1978a] V. V. Markellos, S. Klimopoulos, and A. A. Halioulias, Periodic Motions in the Meridian Plane of a Magnetic Dipole, I, Celestial Mechanics 17, 215 (1978a).
  * Markellos  _et al._ [1978b] V. V. Markellos, S. Klimopoulos, and C. Goudas, Periodic Motions in the Meridian Plane of a Magnetic Dipole, II, Celestial Mechanics 17, 233 (1978b).
  * Wang  _et al._ [2021] R. Wang, S. Liu, A. Xiong, Q.-H. Chen, and F. Zhu, The Origin of Strong Linear Polarization from Fast Radio Bursts, Astrophys. J. 909, 59 (2021), arXiv:2012.09495 [hep-ph] .
  * Horn-Stanja  _et al._ [2018] J. Horn-Stanja, S. Nißl, U. Hergenhahn, T. Sunn Pedersen, H. Saitoh, E. V. Stenson, M. Dickmann, C. Hugenschmidt, M. Singer, M. R. Stoneking, and J. R. Danielson, Confinement of Positrons Exceeding 1 s in a Supported Magnetic Dipole Trap, Phys. Rev. Lett. 121, 235003 (2018).
  * Umberger and Farmer [1985] D. K. Umberger and J. D. Farmer, Fat fractals on the energy surface, Phys. Rev. Lett. 55, 661 (1985).
  * Lee and Flach [2020] W. S. Lee and S. Flach, Deep learning of chaos classification, Machine Learning: Science and Technology 1, 045019 (2020).
  * Piña and Jiménez Lara [1987] E. Piña and L. Jiménez Lara, On the symmetry lines of the standard mapping, Physica D Nonlinear Phenomena 26, 369 (1987).
  * Price  _et al._ [2022] W. Price, M. Formanek, and J. Rafelski, Radiation reaction and limiting acceleration, Phys. Rev. D 105, 016024 (2022), arXiv:2112.04444 [hep-ph] .
  * Gralla  _et al._ [2009] S. E. Gralla, A. I. Harte, and R. M. Wald, Rigorous derivation of electromagnetic self-force, Phys. Rev. D 80, 024031 (2009), arXiv:0905.2391 [gr-qc] .
  * Zehnder [1973] E. Zehnder, Homoclinic points near elliptic fixed points, Communications on Pure and Applied Mathematics 26, 131 (1973), https://onlinelibrary.wiley.com/doi/pdf/10.1002/cpa.3160260204 .
  * Meiss [2017] J. D. Meiss,  _Differential Dynamical Systems_ (SIAM, 2017).
  * Harsoula and Contopoulos [2018] M. Harsoula and G. Contopoulos, Global and local diffusion in the standard map, Phys. Rev. E 97, 022215 (2018), arXiv:1807.06320 [nlin.CD] .
  * Harsoula  _et al._ [2019] M. Harsoula, K. Karamanos, and G. Contopoulos, Characteristic times in the standard map, Phys. Rev. E 99, 032203 (2019), arXiv:1810.11294 [nlin.CD] .
  * Xie and Liu [2020] Y. Xie and S. Liu, From period to quasiperiod to chaos: A continuous spectrum of orbits of charged particles trapped in a dipole magnetic field., Chaos 30 12, 123108 (2020).