Witnessing a Poincar\'e recurrence with Mathematica

TL;DR

This paper demonstrates how to observe Poincaré recurrence in a separable system using Mathematica, leveraging Diophantine approximation and the LLL algorithm to compute large recurrence times accurately.

Contribution

It introduces a method to witness Poincaré recurrence in separable systems through Diophantine approximation and Mathematica's LatticeReduce, providing exact large recurrence times.

Findings

Recurrence times follow the expected scaling law

The method accurately computes very large recurrence times

Demonstration with a harmonic chain system

Abstract

The often elusive Poincar\'e recurrence can be witnessed in a completely separable system. For such systems, the problem of recurrence reduces to the classic mathematical problem of simultaneous Diophantine approximation of multiple numbers. The latter problem then can be somewhat satisfactorily solved by using the famous Lenstra-Lenstra-Lov\'{a}sz (LLL) algorithm, which is implemented in the Mathematica built-in function \verb"LatticeReduce". The procedure is illustrated with a harmonic chain. The incredibly large recurrence times are obtained exactly. They follow the expected scaling law very well.

Equations

$$h(t)=\cos (t)+\cos (\sqrt{2}t)+\cos (\sqrt{3}t)+\cos (\sqrt{5}t)$$

$$t$$

$$H=H(J_1,\dots ,J_N)$$

$${\dot{\theta }}_i=\frac{\partial H}{\partial J_i}\equiv {\omega }_i.$$

$${\mathbb{T}}^N=\{({\theta }_1,\dots ,{\theta }_N)|{\theta }_i\in [0,2\pi ),1\le i\le N\}.$$

$${\mathbb{T}}^{N-1}=\{({\theta }_1,\dots ,{\theta }_{N-1})|{\theta }_i\in [0,2\pi ),1\le i\le N-1\},$$

$$f({\theta }_1,\dots ,{\theta }_{N-1})=({\theta }_1+2\pi {\alpha }_1,\dots ,{\theta }_{N-1}+2\pi {\alpha }_{N-1}),$$

$${\alpha }_i\equiv \frac{{\omega }_i}{{\omega }_N}.$$

$$\Omega =\{({\theta }_1,\dots ,{\theta }_{N-1})||{\theta }_i|\le \pi ϵ\}$$

$$\text{Vol}(\bigcup _{i=0}^K{\Omega }_i)$$

$${\Omega }_m\bigcap {\Omega }_n\ne 0.$$

$${\Omega }_0\bigcap {\Omega }_q\ne 0.$$

$$\underset{1\le i\le N-1}{\max }|{\theta }_i-{\tilde{\theta }}_i|\le 2\pi ϵ,$$

$$⟨q{\alpha }_i⟩\equiv |q{\alpha }_i-\lfloor q{\alpha }_i\rceil |\le ϵ$$

$$\max \{⟨q{\alpha }_1⟩,\dots ,⟨q{\alpha }_n⟩\}\le ϵ.$$

$$L$$

$${\vec{v}}_i$$

$$d(L)$$

$$|\det (\mathcal{V})|=|\det (M)\det (\mathcal{U})|=|\det (\mathcal{U})|.$$

$$d(L)\le \prod _{i=1}^n|{\vec{u}}_i|.$$

$$B_1$$

$$B_2$$

$$|{\vec{b}}_1|\le 2^{(n-1)/4}[d(L){]}^{1/n}.$$

$$A$$

$${\vec{b}}_1$$

$$q\le |{\vec{b}}_1|\le 2^{n/4}{Q}^{n/(n+1)},$$

$$Q|q{\alpha }_i-p_i|\le |{\vec{b}}_1|\le 2^{n/4}{Q}^{n/(n+1)}.$$

$$q$$

$$|q{\alpha }_i-p_i|$$

$$B$$

$${\vec{b}}_1$$

Full text

I Introduction

The theorem established by Poincaré in 1890 is both stunning and intriguing.poincare It states that for any Hamiltonian system with a finite accessible phase space, the system will return to the proximity of a generic initial state infinitely many times. For instance, it would assert that if a cloud of gas initially confined in the left compartment of a vessel is released into the right empty compartment, then after a sufficiently long time, the gas molecules will reassemble in the left compartment. Moreover, they will do so again and again indefinitely.

The power of this paradoxical theorem lies in its generality. It results from the volume-preserving property (Liouville’s theoremlandau ) of the Hamiltonian flow and the finitenessnon of the phase space. The former is a common to any Hamiltonian system and the latter can be easily satisfied in many cases, say, the cloud of molecules in the scenario above. A key insight of Poincaré is that, because of the volume-preserving property, for any initial region $V_{0}$, its time-evolved copies at times [math], $\tau$, $2\tau$, $\ldots$, denoted respectively as $V_{0}$, $V_{1}$, $V_{2}$, $\ldots$, are all of the same volume as $V_{0}$, and since the total accessible phase space has a finite volume, two of these regions must overlap. Let $V_{m}\bigcap V_{n}\neq 0$ ($m<n$), then we have $V_{0}\bigcap V_{n-m}\neq 0$. That is, some points in $V_{0}$ come back into $V_{0}$ at $(n-m)\tau$. Complete proofs and the rigorous statements can be found in the original paper by Poincaré and many books on dynamical systems and/or ergodic theory.walters

While the proof is not complicated and the conclusion inescapable, the predicted recurrence still seems elusive. The problem is that the proof is not constructive—It merely guarantees the existence of the recurrence, but does not tell us when it will occur. In a reply to Zermelo, who used the Poincaré recurrence to question the validity of the $H$-theorem,zermelo ; steckline Boltzmann argued that,boltzmann as is widely accepted today, the recurrence time is extremely large. But how long exactly is it? Can we predict a Poincaré recurrence as we can predict the recurrence of Halley’s comet? Although this is a question out of curiosity, it is not without pedagogical value. Usually, the recurrence time is estimated rudely, with the calculation based more on probability than on mechanical considerations. It is definitely more persuasive to do a rigorous calculation respecting all the mechanical laws and get the number exactly.

The task is apparently a difficult one for a generic system. However, as shown in previous works, for a special class of systems, completely integrable systems to be concrete, the problem is tractable. For such systems, the formidable problem of integrating the equation of motion can be bypassed, and the recurrence problem can be reduced to the classic problem of simultaneous Diophantine approximation of real numbers.frisch ; hemmer ; kim For this classic problem in number theory, there are tons of deep and beautiful theorems,cassels and on the computational side, there is the celebrated Lenstra-Lenstra-Lovász (LLL) algorithm.LLL Fortunately for our purpose, this algorithm is implemented in Mathematica in the function LatticeReduce. By invoking this function, we can obtain effortlessly the astronomically large recurrence times exactly.

Before proceeding, let us recall how the paper was motivated. In 2014, when the first author was a postdoc in the Max Planck Institute for the Physics of Complex Systems, out of curiosity he offered 100 Euro for a number $t$, such that $t>10$ and the function (see Fig. 1)

[TABLE]

is close to $4$ within $10^{-6}$, i.e., $|h(t)-4|\leq 10^{-6}$. The question was posed because of its apparent relevance to the Poincaré recurrence. For $h$ to be close to $4$, each cosine term should be close to unity, or, the four numbers $\{t,\sqrt{2}t,\sqrt{3}t,\sqrt{5}t\}$ should be close to integer multiples of $2\pi$ simultaneously, a condition essentially the same as that of Poincaré recurrence of a completely integrable system. By brute-force, a colleague found the huge number

[TABLE]

and won the money. But the brute-force method apparently is out of the question if there are more cosine terms or if a much higher precision is required, and a more general and sophisticated approach was really what was wanted. Later, we learnt of the LLL algorithm.

II Poincaré recurrence of completely integrable systems

For a generic multi-degree-of-freedom system, its dynamics is chaotic, hence it is hard to follow its motion for a long time, let alone to predict the Poincaré recurrences in the extremely distant future.

More trackable is the motion of completely integrable systems.landau For such a system executing a motion finite in all coordinates, the picture is very simple in terms of the action variables $J_{i}$ and angle variables $\theta_{i}$ ($1\leq i\leq N$). The Hamiltonian

[TABLE]

depends only on the action variables, hence all the action variables are conserved, while all the angle variables increase linearly with time,

[TABLE]

As the angle variables are defined modulo $2\pi$, the system evolves on a torus, the radiuses of which are determined by the actions and the surface of which is parametrized by the angle variables. Explicitly, the torus is

[TABLE]

The flow defined in (4) is volume-preserving and the total volume of the torus $\mathbb{T}^{N}$ is finite, hence the two prerequisites for Poincaré recurrence are satisfied. For a recurrence, all the angle variables should return close to their initial values. In particular, the variable $\theta_{N}$ should return close to its initial value. Since it has a period of $T=2\pi/\omega_{N}$, we take the snapshots of the continuous evolution defined by (4) with a period of $T$, and consider the resulting discrete dynamical system defined on the torus

[TABLE]

namely,

[TABLE]

where the variables $\alpha_{i}$ are defined as ($1\leq i\leq N-1$)

[TABLE]

Apparently, $f$ is also volume-preserving.

Now take a small neighborhood of the origin

[TABLE]

with $\epsilon\ll 1$. Consider the images of $\Omega$ under the iterated action of $f$, i.e., $\Omega_{i}=f^{i}(\Omega)$, $0\leq i\leq\infty$. They are all of the same volume. Hence, for $K=\lfloor 1/\epsilon^{N-1}\rfloor$, where $\lfloor x\rfloor$ denotes the largest integer no more than $x$, the $K+1$ regions $\{\Omega_{i}|0\leq i\leq K\}$ cannot be mutually disjoint. Otherwise,

[TABLE]

Hence there must be $0\leq m<n\leq K$ such that

[TABLE]

Now by applying $f^{-m}$ to both sides, we get ($q=n-m\leq K$)

[TABLE]

This means recurrence. Specifically, let $(\tilde{\theta}_{1},\ldots,\tilde{\theta}_{N-1})\in\Omega_{0}\bigcap\Omega_{q}$. Since this point belongs to $\Omega_{q}$, there exists a point $(\theta_{1},\ldots,\theta_{N-1})\in\Omega_{0}$ which, when translated $q$ times by $f$, arrives at $(\tilde{\theta}_{1},\ldots,\tilde{\theta}_{N-1})\in\Omega_{0}$, i.e., it comes back into $\Omega_{0}$. The relative error of recurrence is bounded by

[TABLE]

as $\{|\theta_{i}|,|\tilde{\theta_{i}}|\}\leq\pi\epsilon\ll\pi$. As $f$ translates the torus as a whole rigidly, when one point returns to the proximity of its initial position, so do all other points. That is, the system not only recurs, but even recurs uniformly, or in other words, independent of the initial condition. We have thus not only proved the existence of recurrence, but also obtained an upper bound of the recurrence time in terms of the precision $\epsilon$ and the number of degrees of freedom. The scaling law anticipated is that the recurrence time is on the order of $1/\epsilon^{N-1}$ for a completely integrable system with $N$ degrees of freedom.

It is easily seen that the condition (12) is equivalent to the condition that

[TABLE]

for each $1\leq i\leq N-1$. Here, $\lfloor x\rceil$ denotes the integer nearest to the real number $x$ (for instance, $\lfloor 2.4\rceil=2$ and $\lfloor 2.7\rceil=3$), and $\langle x\rangle$ denotes the distance of $x$ to $\lfloor x\rceil$. By definition, $\langle x\rangle\leq 1/2$ for an arbitrary $x$.

Equation (14) is the sufficient condition for a recurrence with a prescribed precision of $\epsilon$. The problem of looking for a Poincaré recurrence is then reduced to the problem of looking for an integer $q$ such that (14) is satisfied simultaneously for all the numbers $\{\alpha_{1},\ldots,\alpha_{N-1}\}$. It turns out that this is a classic problem termed simultaneous Diophantine approximation in number theory.cassels For this problem, there is the celebrated Dirichlet theorem, which states that for $n$ arbitrary real numbers $\{\alpha_{1},\ldots,\alpha_{n}\}$, and any $\epsilon>0$, there exists a positive integer $q\leq 1/\epsilon^{n}$, such that

[TABLE]

This is essentially the same as what is proven above.

III LLL algorithm

The problem is now to find a $q$ such that (14) is satisfied for given $\epsilon$, which will be arbitrarily small. For this purpose, a useful tool is the famous LLL algorithm, which is implemented in Mathematica by the function LatticeReduce. As the name suggests, the algorithm is about the lattice, a most basic concept in solid state physics. Suppose $\{\vec{u}_{1},\ldots,\vec{u}_{n}\}$ are $n$ linearly independent vectors in $\mathbb{R}^{n}$. The lattice determined by them is the subset of $\mathbb{R}^{n}$ defined as

[TABLE]

The vectors $\{\vec{u}_{1},\ldots,\vec{u}_{n}\}$ are called a basis of the lattice. We say they span the lattice. The point is that a basis determines a lattice, but not vice versa. In other words, for a given lattice, the basis is not unique. Actually, for any $n\times n$ matrix $M$, whose entries are all integral and whose determinant is $\pm 1$, the vectors

[TABLE]

constitute another basis of $L$. The proof is simple. By (17), the lattice spanned by the $\vec{v}$’s, denoted as $L^{\prime}$, is contained in $L$. However, by construction, the inverse of $M$ has the same property as $M$, and hence the $\vec{u}$’s are also integral linear combinations of the $\vec{v}$’s, which means that $L$ is contained in $L^{\prime}$. Therefore, they must be equal.

For a lattice, a characteristic quantity is its determinant, which is defined as

[TABLE]

Here $\mathcal{U}\equiv(\vec{u}_{1};\ldots;\vec{u}_{n})$ is the $n\times n$ matrix with $\vec{u}_{i}$ being the $i$th row. Despite of the explicit reference to the basis in its definition, the quantity $d(L)$ does not depend on the choice of the basis actually. The reason is simple—For another basis $\{\vec{v}_{i}\}$ of $L$, the matrix $\mathcal{V}\equiv(\vec{v}_{1};\ldots;\vec{v}_{n})$ relates to the matrix $\mathcal{U}$ as $\mathcal{V}=M\mathcal{U}$ by (17), and therefore,

[TABLE]

The physical or geometric meaning of $d(L)$ is the volume of the parallelepiped subtended by the $n$ basis vectors. Hence, we have the Hadamard’s inequalitymatrix

[TABLE]

The equality is achieved if and only if $\vec{u}_{i}$ are orthogonal to each other.

Although all basis are equivalent in spanning the lattice, some basis is more preferable than others. For example, in Fig. 2, the lattice is spanned either by the basis

[TABLE]

or by the basis

[TABLE]

Apparently, the latter is more appealing as the basis vectors involved are shorter and more close to being orthogonal. We thus see that the defining basis of a lattice might be far from being optimal and generally there is room to reduce or simplify it, in the sense that $B_{2}$ is reduced compared to $B_{1}$. This is exactly what the LLL algorithm is about.

There are numerous literature about the LLL algorithm. In particular, the original paper is very readable.LLL However, here for our purpose, we just need to know that the LLL algorithm, when fed with an original basis $\{\vec{a}_{1},\ldots,\vec{a}_{n}\}$, will return us in polynomial time a new basis $\{\vec{b}_{1},\ldots,\vec{b}_{n}\}$ such that

[TABLE]

That the new basis is somehow reduced can be glimpsed by comparing (23) with (20). The determinant $d(L)$ and the basis vector lengths swapped position.

The way the LLL algorithm is used for simultaneous Diophantine approximation is as follows. Let $\{\alpha_{1},\ldots,\alpha_{n}\}$ be the $n$ real numbers. Construct the $(n+1)\times(n+1)$ matrix

[TABLE]

where $Q$ is some large integer to be chosen according to the precision wanted. The $(n+1)$ rows are the basis of a lattice. The determinant of the lattice is apparently $d(L)=Q^{n}$. The first vector $\vec{b}_{1}$ in the reduced basis is of the form

[TABLE]

where $q$ and $p_{i}$ ($1\leq i\leq n$) are integers. By (23), we have

[TABLE]

Let $Q=2^{n(n+1)/4}C^{n+1}$, we have

[TABLE]

Hence, we have found a $q$ such that (15) is satisfied with $\epsilon=1/C$. This $q$ is not bounded by $1/\epsilon^{n}$ but is on the same order.

In practice, the LLL algorithm is generally implemented with the components of the basis vectors being rational numbers. Hence, a small modification is necessary.hanrot Instead of (28), we construct the matrix

[TABLE]

where the rounding function $\lfloor\cdot\rceil$ is implemented in Mathematica with Round. In this case,

[TABLE]

By (23), we have

[TABLE]

As for $|q\alpha_{i}-p_{i}|$, we have

[TABLE]

Here in the second line, we used the inequality $|\alpha_{i}-\lfloor\alpha_{i}\rceil|\leq 1/2$, and in the third line, we used the Cauchy-Schwarz inequality. Taking $Q=(\frac{\sqrt{5}}{2})^{n+1}2^{n(n+1)/4}C^{n+1}$, we have similar to (31),

[TABLE]

The scaling law is unchanged. This is the very algorithm we shall use to search the exact recurrence times.

IV A case study

We now demonstrate the theory described above with a concrete model, which is illustrated in Fig. 3. The $N$ point-masses are each of mass $m$ and the springs are all of stiffness $k=m\omega^{2}$ (we shall set $m=k=\omega=1$ when it comes to numerics). The Hamiltonian is

[TABLE]

Here $\mathbf{x}\equiv(x_{1},x_{2},\ldots,x_{N})^{T}$ denotes the displacements of the point-masses with respect to their equilibrium positions and $\mathbf{p}\equiv(p_{1},p_{2},\ldots,p_{N})^{T}$ denotes their momenta. The $N\times N$ matrix $\mathbf{M}$ is tridiagonal as

[TABLE]

It is symmetric and thus can be diagonalized. That is, we have

[TABLE]

with $\mathbf{D}$ being a diagonal matrix and $\mathbf{U}$ an orthogonal matrix. The explicit forms of $\mathbf{D}$ and $\mathbf{U}$ are

[TABLE]

with $1\leq i,j\leq N$. Substituting (48) and (49) into (41), we get

[TABLE]

Here we have introduced the new coordinates $\mathbf{X}\equiv(X_{1},X_{2},\ldots,X_{N})^{T}$ and momenta $\mathbf{P}\equiv(P_{1},P_{2},\ldots,P_{N})^{T}$,

[TABLE]

The Hamiltonian is thus diagonalized as the sum of $N$ independent harmonic oscillators, with the frequency of the $i$th oscillator being

[TABLE]

The evolution of the system is then simply that

[TABLE]

where the amplitude $A_{i}$ and the phase $\phi_{i}$ are determined by the initial condition. We thus have the $N$-torus formalism and the algorithm can be applied directly.

The procedure to get a recurrence time is then as follows. Let $\alpha_{i}=\omega_{i}/\omega_{N}$, $1\leq i\leq N-1$. Choose a large integer $Q$, and construct the $B$ matrix as in (36). Feed it into the function LatticeReduce, and get a reduced basis. The corresponding code is

RB = LatticeReduce(B);

The new basis vectors are the rows of RB, and the integer $q$ we want is the first component of the first basis vector, i.e.,

q = RB[[1,1]];

Once the $q$ is obtained, the recurrence time is

[TABLE]

Here the subscript means Poincaré recurrence.

How good the algorithm works can be checked with some concrete numbers. Let $N=15$, and $Q=10^{35}$, we get

[TABLE]

and the recurrence error

[TABLE]

The error is sufficiently small, hence an arbitrary initial state will recur after time $T_{\text{pr}}=2\pi q/\omega_{N}$ to a good precision. Consider the initial state defined as

[TABLE]

Note that we have deliberately chosen an artificial state, in which only the $k$th point-mass is displaced and has a nonzero velocity. In Fig. 4, four snapshots of the coordinates and momenta of the point-masses are shown, with the corresponding times being $t=T_{\text{pr}}-200$, $t=T_{\text{pr}}$, and $t=T_{\text{pr}}\pm 3$. We see that at $t=T_{\text{pr}}-200$, almost all the point-masses are significantly displaced and have gained significant momenta. This is typical of the configuration of the system at an arbitrary time. The initially localized motion has spread over the whole system and the system remains so for a long time. However, around $t=T_{\text{pr}}$, the system exhibits a recurrence. At $t=T_{\text{pr}}-3$, the tendency is clear, the motion regathers around the $k$th point-mass, and the point-masses far away have already returned to their initial positions and come to rest. Finally, at $t=T_{\text{pr}}$, the system restores itself to its initial configuration to such good extent that the difference between the current values and the initial values of $(x_{n},p_{n})$ is hardly visible. Afterwards, as shown by the snapshots at $t=T_{\text{pr}}+3$, the motion disperses again and we have to wait yet another long time to witness a second recurrence.

Now for each $Q$ we can get a $q$ and a corresponding $error$. By taking $Q$ larger and larger, we get a series of pairs of $(q,error)$. It is necessary and straightforward to check the scaling relation between them. This we do in Fig. 5. In Fig. 5(a), we have $N=15$, and $q$ scales with $error$ as $q\propto error^{-14}$. However, in Fig. 5(b), we have $N=5$, and $q$ scales with $error$ not as $q\propto error^{-4}$ but as $q\propto error^{-3}$. The reason is simply that for $N=5$, the natural frequencies (52) of the normal modes are not linearly independent over the integers. To be specific, as can be easily verified, we have the equality,

[TABLE]

Hence, for an arbitrary integer $q$, $q=q\alpha_{1}+q\alpha_{3}$, which means, whenever $\langle q\alpha_{1}\rangle$ is less than $\epsilon$, so is $\langle q\alpha_{3}\rangle$ automatically. Therefore, the problem of looking for a $q$ such that

[TABLE]

is actually equivalent to the problem of looking for a $q$ such that

[TABLE]

This is why we have the $q\propto error^{-3}$ law.

V Discussions

We have had the privilege to witness a Poincaré recurrence, which is a very improbable event—much rarer than the recurrence of Halley’s comet for any realistic system actually. While this seems a hopeless task for a generic system, it is possible in the case of a completely integrable system, thanks to the fact that the problem can be reduced to the classic Diophantine approximation problem, for which mathematicians have already prepared many theorems, and in particular, algorithms, for us. We have thus here a case in which number theory plays a vital role in solving physical problems. Note that while analysis and algebra are routinely used in physics nowadays, number theory presents itself rarely in physics.

But frankly speaking, we have achieved only partial success. While the LLL algorithm is powerful, its insufficiency is also apparent. For given $\epsilon$ and for a generic set of numbers $\{\alpha_{1\leq i\leq n}\}$, there are infinitely many $q$ satisfying (15). Actually, collecting the appropriate $q$’s in a sequence $\{q_{1},q_{2},\ldots,\}$, Weyl’s theorem states thatstein

[TABLE]

That is, on average, the distance between two adjacent $q$’s is $1/\epsilon^{n}$. Physically, it means that recurrences with a precision on the order of $\epsilon$ will appear with a frequency on the order of $\epsilon^{n}$. However, the algorithm can get only one $q$ for each $\epsilon$. How to capture all the $q$’s is thus a problem of interest. As far as we know, this problem has been completely solved only for $n=1$.kovacs For larger $n$, only partial success has been achieved. Moreover, the $q$ returned by the algorithm is often not the smallest one, i.e., it is not smaller than $1/\epsilon^{n}$, the bound given by the Dirichlet theorem, although it is on the same order.

Historically, Poincaré recurrence was discovered in the context of classical mechanics. However, it can be generalized to quantum mechanics straightforwardly.loinger ; schulman ; wallace For the sake of simplicity, let us assume that the quantum system has a finite-dimensional Hilbert space, and the eigenstates and eigenvalues are $\{|m\rangle\}$ and $\{E_{m}\}$, respectively. The quantum Poincaré recurrence theorem then states that for any initial state

[TABLE]

the system, evolving as

[TABLE]

will come back close to $|\psi_{0}\rangle$ up to an arbitrary precision infinitely many times. Specifically, for any $t>0$ and $\epsilon>0$, there exists a $t_{r}>t$ such that

[TABLE]

This is not surprising. As (63) shows, the wave function also evolves on an $N$-torus, and hence the formalism developed for a completely integrable system applies too. Quantitatively, we have

[TABLE]

For (64) to be satisfied for a generic initial state $|\psi_{0}\rangle$, a sufficient condition is

[TABLE]

for all $1\leq m\leq N$. Let $t=2\pi q/E_{N}$ with $q$ a whole number, so that $\langle E_{N}t/2\pi\rangle=0$ and (66) is satisfied already for $m=N$, we get the problem defined in (14) again, with $\alpha_{i}=E_{i}/E_{N}$. Therefore, the same algorithm applies for a quantum system too.

ACKNOWLEDGMENTS

We are grateful to Yang Ji and Hai-feng Song for helpful comments. This work is supported by the Fujian Provincial Science Foundation under Grant No. 2016J05004 and the Science Challenge Project under Grant No. JCKY2016212A502.

Appendix

Here, we append the Mathematica code used to find the particular number in (55):

num = 15; omega = 2 * Sin[Range[num]Pi/2/(num + 1)]; alpha = Take[omega, {1, num - 1}]/omega[[num]]; Q = 10^35; B = Table[0, {i, 1, num}, {j, 1, num}]; B[[1, 1]] = 1; For[i = 1, i <= num - 1, i++, B[[1, i + 1]] = Round[Qalpha[[i]]];]; For[i = 2, i <= num, i++, B[[i, i]] = Q; ]; RB = LatticeReduce[B]; q = Abs[RB[[1, 1]]] error = Abs[N[qalpha - Round[q alpha], 200]];

Once the most crucial number $q$ is obtained, all other calculations can be done straightforwardly.

The readers can modify the code a little bit to reproduce the particular number in Eq. (2).