Exact Splitting Methods for Kinetic and Schrödinger Equations

Abstract

In (Bernier in Exact splitting methods for semigroups generated by inhomogeneous quadratic differential operators. arXiv:1912.13219, (2019)), some exact splittings are proposed for inhomogeneous quadratic differential equations including, for example, transport equations, Fokker–Planck equations, and Schrödinger type equations with an angular momentum rotation term. In this work, these exact splittings are used combined with pseudo-spectral methods in space. High accuracy and efficiency of exact splitting methods are illustrated and comparison are performed with the numerical methods in literature. We show that our methods can be used to improve significantly some classical splitting methods for some nonlinear or non-quadratic equations.

Appendix

1.1 Backward Error Analysis of the Strang Splittings for Rotations

We consider the ODE associated with a rotation in ${\mathbb {R}}^n$

$$\begin{aligned} {\dot{y}} = By \end{aligned}$$

(37)

where $B\in {\mathfrak {so}}_n({\mathbb {R}})$ is a real skew-symmetric matrix of size n. For $t\in {\mathbb {R}}$, the Strang splitting naturally associated with the decomposition of this rotation in shear transforms is

$$\begin{aligned} P_t = \left( \prod _{j=1}^{n-1} \mathrm {I}_n + \frac{t}{2} (e_j\otimes e_j) B \right) \left( \mathrm {I}_n + t (e_n\otimes e_n) B \right) \left( \prod _{j=1}^{n-1} \mathrm {I}_n + \frac{t}{2} (e_{n-j}\otimes e_{n-j}) B \right) . \end{aligned}$$

where $e_1,\dots ,e_n$ is the canonical basis of ${\mathbb {R}}^n$. It is indeed a Strang splitting because since $(e_j\otimes e_j) B$ is nilpotent of order 1 we have

$$\begin{aligned} e^{t (e_j\otimes e_j) B} = \mathrm {I}_n + t (e_j\otimes e_j) B. \end{aligned}$$

The following proposition states that, up to a change of coordinate close to the identity, $P_t$ is a rotation of which angles coincide with those of $\exp (tB)$ up to an error of order $t^3$.

Proposition 2

There exists $t_0>0$ and two analytic functions $t\in (-t_0,t_0) \mapsto (B_t,Q_t) \in \mathfrak {so}_n({\mathbb {R}}) \times \mathrm {S}_n({\mathbb {R}})\cap \mathrm {GL}_n({\mathbb {R}})$, where $B_t$ is real skew symmetric and $Q_t$ is real symmetric and invertible such that $B_0 = B$, $Q_0 = \mathrm {I}_n$ and

$$\begin{aligned} \forall t \in (-t_0,t_0) , \ P_t = Q_{t^2} e^{t B_{t^2}} Q_{t^2}^{-1}. \end{aligned}$$

Proof

Observing that since

$$\begin{aligned} \mathfrak {sp}(B) = \{ SB \ | \ S \in \mathrm {S}_{n}({\mathbb {R}}) \} \ \mathrm {is \ a \ Lie \ algebra}, \end{aligned}$$

it follows from the Bakker-Campbell-Hausdorff formula that there exist $t_1>0$ and an analytic function $t \in (-t_1,t_1)\mapsto S_t \in \mathrm {S}_{n}({\mathbb {R}})$ such that $S_0= \mathrm {I}_n$ and

$$\begin{aligned} \forall t\in (-t_1,t_1), \ P_t = e^{tS_tB}. \end{aligned}$$

Since the Strang splitting is reversible, we have $P_t = P_{-t}^{-1}$. Thus, since the exponential map is injective in a neighborhood of the identity, we deduce that $t\mapsto S_t$ is an even function. Consequently, there exists an analytic function $t\mapsto K_t$ such that $K_{t^2} = S_t$. Furthermore, since $K_0= \mathrm {I}_n$, there exists $t_0\in (0,t_1)$ such that for all $t\in (-t_0,t_0)$, $K_t$ is positive-definite. Finally, we deduce that,

$$\begin{aligned} \forall t\in (-t_0,t_0), \ P_t = e^{tS_tB} = \sqrt{K_{t^2}}e^{t\sqrt{K_{t^2}}B\sqrt{K_{t^2}}}\sqrt{K_{t^2}}^{-1}= Q_{t^2} e^{t B_{t^2}} Q_{t^2}^{-1} \end{aligned}$$

where $Q_t = \sqrt{K_{t}}$ and $B_t = \sqrt{K_{t}}B\sqrt{K_{t}}$. $\square$

1.2 2D Magnetic Schrödinger Equation

$$\begin{aligned} i\epsilon \partial _t \psi (\mathbf{x},t) = -\frac{\epsilon ^2}{2} \Delta \psi (\mathbf{x},t) + i\epsilon \mathbf{A} \cdot \nabla \psi (\mathbf{x},t) + \frac{1}{2}{|\mathbf{A}|}^2\psi (\mathbf{x},t), \end{aligned}$$

(38)

where $\mathbf{x} = (x_1,x_2) \in {\mathbb {R}}^2$, $\mathbf{A} = \frac{1}{2}(A_1, A_2)$, $A_1= -x_2$, $A_2 = x_1$. The above system can be split into three systems:

$$\begin{aligned}&i\epsilon \partial _t \psi (\mathbf{x},t) = -\frac{\epsilon ^2}{2} \Delta \psi (\mathbf{x},t), \end{aligned}$$

(39)

$$\begin{aligned}&\partial _t \psi (\mathbf{x},t) = \mathbf{A} \cdot \nabla \psi (\mathbf{x},t), \end{aligned}$$

(40)

$$\begin{aligned}&i\epsilon \partial _t \psi (\mathbf{x},t) = \frac{1}{2}{|\mathbf{A}|}^2\psi (\mathbf{x},t), \end{aligned}$$

(41)

The solutions of the above three subsystems can be obtained by operators $e^{it\frac{\epsilon }{2} \Delta }$, $e^{t\text {Rot}}$, and $e^{tV}$ respectively. Since the second is nothing but a 2D rotation, we call the associated solution $e^{t\text {Rot}}$. Then we have the following second order splitting method

$$\begin{aligned} \psi ^{n+1} = e^{\frac{\Delta t}{2}V} e^{i \Delta t\frac{\epsilon }{4} \Delta } e^{\Delta t\text {Rot}} e^{i\Delta t\frac{\epsilon }{4} \Delta } e^{\frac{\Delta t}{2}V}, \end{aligned}$$

(42)

from which we derive two variants according to the treatment of $e^{\Delta t\text {Rot}}$. Indeed, ESR denotes the splitting method (42) when $e^{\Delta t \text {Rot}}$ is solved by exact splittings for transport equation in Proposition 1. Strang denotes (42) when $e^{\Delta t \text {Rot}}$ is approximated by Strang directional splitting.

1.3 2D Rotating Gross–Pitaevskii Equation

The rotating Gross-Pitaevskii equation (GPE) [9, 11] is

$$\begin{aligned} \begin{aligned}&{i}\partial _t \psi (\mathbf{x}, t) = -\frac{1}{2}\Delta \psi (\mathbf{x},t) + V(\mathbf{x}) \psi (\mathbf{x},t) + \beta | \psi |^2 \psi (\mathbf{x},t) - \Omega L_{x_3} \psi (\mathbf{x},t),\ {\mathbf{x }} \in {{\mathbb {R}}}^2, \end{aligned} \end{aligned}$$

(43)

where $\psi (\mathbf{x}, t)$ is the macroscopic wave function, $\mathbf{x} = (x_1, x_2)$, $L_{x_3} = - i (x_1\partial _{x_2} - x_2 \partial _{x_1})$. Two operator splittings are presented to approximate (43).

1.3.1 BW Method

Here we recall the splitting method introduced in [11] to approximate (43). We will call it BW in the sequel. BW splitting for rotating GPE (43) is based on the following two-steps splitting

$$\begin{aligned}&i\partial _t \psi (\mathbf{x}, t) = -\frac{1}{2}\Delta \psi (\mathbf{x}, t) - \Omega L_{x_3} \psi (\mathbf{x}, t), \end{aligned}$$

(44)

$$\begin{aligned}&\partial _t \psi (\mathbf{x}, t) = V(\mathbf{x})\psi (\mathbf{x}, t) + \beta |\psi (\mathbf{x}, t)|^2\psi (\mathbf{x}, t). \end{aligned}$$

(45)

Then, the authors in [11] noticed that (44) can be split further as

$$\begin{aligned}&i\partial _t \psi (\mathbf{x}, t) = -\frac{1}{2}\partial ^2_{x_1}\psi (\mathbf{x}, t) -i\Omega x_2 \partial _{x_1} \psi (\mathbf{x}, t), \end{aligned}$$

(46)

$$\begin{aligned}&i\partial _t \psi (\mathbf{x}, t) = -\frac{1}{2}\partial ^2_{x_2}\psi (\mathbf{x}, t) + i\Omega x_1 \partial _{x_2} \psi (\mathbf{x}, t). \end{aligned}$$

(47)

The solutions of subsystems (45), (46) and (47) can be obtained by operators $e^{tN}, e^{tX}$ and $e^{tY}$ respectively, the second order BW method is then derived from the following composition

$$\begin{aligned} \begin{aligned} \psi ^n(\mathbf{x})&= \left( e^{\Delta t/2 \, Y}e^{\Delta t/2\,X}e^{\Delta t\, N}e^{\Delta t/2\,X}e^{\Delta t/2\, Y} \right) ^n \psi _0(\mathbf{x}),\\&= e^{\Delta t/2\, Y} (e^{\Delta t/2\, X}e^{\Delta t\, N}e^{\Delta t/2\, X}e^{\Delta t\, Y})^{n-1}e^{\Delta t/2\, X}e^{\Delta t\, N}e^{\Delta t/2\, X}e^{\Delta t/2\, Y} \psi _0(\mathbf{x}). \end{aligned} \end{aligned}$$

(48)

Combined with Fourier pseudo-spectral method in space, we can see that in each time step, we need six calls to FFT.

1.3.2 Lagrangian Method

Here we recall the main step of the method introduced in [10] to approximate (43). We will call it Lagrangian in the sequel. First, a change of coordinates is considered

$$\begin{aligned} \phi (\tilde{\mathbf{x}}, t) := \psi (\mathbf{x}, t) = \psi ( e^{Jt} \mathbf{x}, t), \end{aligned}$$

where J is the $2\times 2$ symplectic matrix. Then, $\phi$ satisfies the following equation

$$\begin{aligned} \partial _t \phi (\tilde{\mathbf{x}}, t) = - \frac{1}{2}\Delta \phi (\tilde{\mathbf{x}},t) + W(\tilde{\mathbf{x}}, t) \phi (\tilde{\mathbf{x}},t) + \beta | \phi (\tilde{\mathbf{x}},t) |^2 \phi (\tilde{\mathbf{x}},t),\ \tilde{\mathbf{x }} \in {{\mathbb {R}}}^2, \end{aligned}$$

(49)

with the initial condition $\phi (\tilde{\mathbf{x}}, 0)=\psi _0(\mathbf{x})$ and where W is defined by

$$\begin{aligned} W(\tilde{\mathbf{x}}, t) = V(e^{Jt}\tilde{\mathbf{x}}), \end{aligned}$$

so that a time dependency is created. However, if V is isotropic, $W(\tilde{\mathbf{x}}, t) = V(\tilde{\mathbf{x}})$ is time-independent. The main advantage of (49) does not involve the angular momentum rotation term and the following splitting is used

$$\begin{aligned}&i\partial _t \phi (\tilde{\mathbf{x}}, t) = -\frac{1}{2}\Delta \phi (\tilde{\mathbf{x}}, t), \end{aligned}$$

(50)

$$\begin{aligned}&\partial _t \phi (\tilde{\mathbf{x}}, t) = W(\tilde{\mathbf{x}}, t)\phi (\tilde{\mathbf{x}}, t) + \beta |\phi (\tilde{\mathbf{x}}, t)|^2\phi (\tilde{\mathbf{x}}, t). \end{aligned}$$

(51)

For harmonic potential V, each step can be solved exactly (which is not the case for general potential V). Combined with a Strang splitting and with Fourier pseudo-spectral method in space, we can see that in each time step, this method needs four calls to FFT.

1.4 3D Time-Periodic Quadratic Linear Schrödinger Equation

For (23) with $f=0$ and B and V are specified in (32) and (33), we consider two numerical methods: ESQM and a standard Strang operator splitting.

1.4.1 Exact Splitting

The coefficients for ESQM (21) are given by

$$\begin{aligned}&A_{\Delta t} \simeq \begin{pmatrix} 0.503369336514750 &{}0.09260872887966 &{}-0.086577853155386 \\ 0.092608728879667 &{} 0.499175997238123&{}0.090475411725230 \\ -0.086577853155386 &{}0.090475411725230 &{}0.482430618251455 \end{pmatrix}, \\&V_{\Delta t}^{(\ell )} \simeq \begin{pmatrix} 1.838313777101704 &{} 0&{} 0 \\ 0 &{}1.405233579215994 &{}0 \\ 0 &{}0 &{}2.416160688906186 \end{pmatrix}, \\&V_{\Delta t}^{(r)} \simeq \begin{pmatrix} 0.765638127548775 &{}0.097739062052903 &{} -0.244124321719139 \\ 0.097739062052903 &{}1.408683914880933 &{} 0.141925135897144 \\ -0.244124321719139 &{}0.14192513589714 &{} 3.535113753227984 \end{pmatrix}, \\&L_{\Delta t} \simeq \begin{pmatrix} 0 &{} 0&{}0 \\ 0.957867410476376 &{}0 &{} 0 \\ -0.917880413070041 &{}1.133563918623215 &{} 0 \end{pmatrix}, \\&U_{\Delta t} \simeq \begin{pmatrix} 0 &{} -1.132325985517193&{} 0.915677911046419 \\ 0 &{} 0 &{} -0.957661219232001 \\ 0 &{} 0&{} 0 \end{pmatrix}. \end{aligned}$$

1.4.2 Strang Method

Classically, we use the following operator splitting

$$\begin{aligned}&i \frac{\partial \psi (\mathbf{x}, t)}{\partial t} = -\frac{1}{2}\Delta \psi (\mathbf{x},t),\\&\frac{\partial \psi (\mathbf{x}, t)}{\partial t} =- ({ B} \mathbf{x})\cdot \nabla \psi (\mathbf{x},t),\\&i \frac{\partial \psi (\mathbf{x}, t)}{\partial t} = V(\mathbf{x})\psi (\mathbf{x},t). \end{aligned}$$

The solutions of the above three subsystems can be obtained by operators $e^{-it\frac{1}{2}\Delta }$, $e^{t \text {Rot}}$, and $e^{-it \text {V}}$ respectively so that we have the following second order splitting method

$$\begin{aligned} \psi ^{n+1}({\mathbf{x }}) = e^{-i\frac{\Delta t}{2} \text {V}} e^{-i\Delta t\frac{1}{4}\Delta } e^{\Delta t \text {Rot}} e^{-i\Delta t\frac{1}{4}\Delta }e^{-i\frac{\Delta t}{2} \text {V}} \psi ^n({\mathbf{x }}). \end{aligned}$$

(52)

Strang denotes (52) when $e^{\Delta t \text {Rot}}$ is also approximated by a Strang directional splitting.

1.5 3D Magnetic Schrödinger Equation

From (35), where $\mathbf{A}(\mathbf{x}) = \mathbf{x} \times \mathbf{B}$, and $V_{nq}$ is given by (36), we can use the following operator splitting

$$\begin{aligned}&i \frac{\partial \psi (\mathbf{x}, t)}{\partial t} = -\frac{1}{2}\Delta \psi (\mathbf{x},t),\\&\frac{\partial \psi (\mathbf{x}, t)}{\partial t} = {\mathbf{A }} (\mathbf{x})\cdot \nabla \psi (\mathbf{x},t), \\&i \frac{\partial \psi (\mathbf{x}, t)}{\partial t} = \frac{1}{2}|\mathbf{A}(\mathbf{x})|^2\psi (\mathbf{x},t) + V_{nq}(\mathbf{x})\psi (\mathbf{x},t), \end{aligned}$$

The solutions of the above three subsystems can be obtained by operators $e^{-it\frac{1}{2}\Delta }$, $e^{t \text {Rot}}$, and $e^{t \text {VA}}$ respectively and we can derive a second order splitting method:

$$\begin{aligned} \psi ^{n+1}({\mathbf{x }}) = e^{\frac{\Delta t}{2} \text {VA}} e^{-i\Delta t\frac{1}{4}\Delta } e^{\Delta t \text {Rot}} e^{-i\Delta t\frac{1}{4}\Delta }e^{\frac{\Delta t}{2} \text {VA}} \psi ^n({\mathbf{x }}). \end{aligned}$$

(53)

ESR denotes the splitting method (53) when $e^{\Delta t \text {Rot}}$ is solved by exact splittings for transport equation in Proposition 1. Strang denotes (53) when $e^{\Delta t \text {Rot}}$ is approximated by Strang directional splitting.

The coefficients when $\Delta t = 0.1$ for ESQM (25) are as follows

$$\begin{aligned}&A_{\Delta t} \simeq \begin{pmatrix} 0.506160069704187 &{}0.098840554692409&{}0.001683128724191 \\ 0.098840554692409 &{} 0.508317718832811&{}0.050167780151672 \\ 0.001683128724191 &{}0.050167780151672 &{}0.501715861437068 \end{pmatrix}, \\&V_{\Delta t}^{(\ell )} \simeq \begin{pmatrix} 2.025343613765655 &{} 0&{} 0 \\ 0 &{}0.508168767491105 &{}0 \\ 0 &{}0 &{}0.000099459606977 \end{pmatrix}, \\&V_{\Delta t}^{(r)} \simeq \begin{pmatrix} 0.072891278447532 &{}0.242556937819776 &{} -1.026420948565178 \\ 0.242556937819776 &{}1.959142247295385 &{} -0.046535665904951 \\ -1.026420948565178 &{}-0.046535665904951 &{} 0.508102737430800 \end{pmatrix}, \\&L_{\Delta t} \simeq \begin{pmatrix} 0 &{} 0&{}0 \\ 2.003434507092443 &{}0 &{} 0 \\ -0.099043028107977 &{}1.016569585390557&{} 0 \end{pmatrix}, \\&U_{\Delta t} \simeq \begin{pmatrix} 0 &{} -1.963756896350695&{} -0.099988990937417 \\ 0 &{} 0 &{} -1.006635420674690 \\ 0 &{} 0&{} 0 \end{pmatrix}. \end{aligned}$$

1.6 Proof of the Period 360

Lemma 1

The function $t\mapsto U_t = e^{ i t ( \Delta /2 - V(x)) - t Bx \cdot \nabla }$, where V and B are given by (32) satisfies

$$\begin{aligned} \forall t\in {\mathbb {R}}, \ U_{t+180} = -U_t. \end{aligned}$$

Proof

Since $t\mapsto U_t$ is a group, we just have to prove that

$$\begin{aligned} U_{180}=-I_{L^2({\mathbb {R}}^3)}. \end{aligned}$$

We recall that by construction, we have $U_t = e^{-t q_{(\hbox {QM})}^w}$ where

$$\begin{aligned} q_{(\hbox {QM})}^w = i\frac{|\varvec{\xi }|^2}{2} + iB \mathbf{x } \cdot {\varvec{\xi }} + iV(\mathbf{x }) \end{aligned}$$

Step 1: To conjugate $q_{(\hbox {QMS})}^w$ to a sum of harmonic oscillators. We are going to prove that there exists $V\in {\mathcal {U}}(L^2({\mathbb {R}}^n))$ such that

$$\begin{aligned} U_t = V \exp (-it \sum _{j=1}^3 \omega _j (x_j^2 -\partial _{x_j}^2) ) V^* \end{aligned}$$

(54)

where $(\omega _1,\omega _2,\omega _3) = \frac{\pi }{180}(20,75,132)$. Assuming first this decomposition, we deduce that

$$\begin{aligned} U_{180} = V \exp (-20 i \pi (x_1^2 -\partial _{x_1}^2) ) \exp (-75 i \pi (x_2^2 -\partial _{x_2}^2) ) \exp (-132 i \pi (x_3^2 -\partial _{x_3}^2) )V^*. \end{aligned}$$

But, in dimension 1, the eigenvalues of the harmonic oscillator $x^2-\partial _x^2$ being the odd positive integers, we know that $\exp (i\pi (x^2-\partial _x^2) ) = -I_{L^2({\mathbb {R}}^3)}$. Thus, we deduce that

$$\begin{aligned} U_{180} = V I_{L^2({\mathbb {R}}^3)} (-I_{L^2({\mathbb {R}}^3)}) I_{L^2({\mathbb {R}}^3)} V^* = -I_{L^2({\mathbb {R}}^3)}. \end{aligned}$$

In order to prove (54) we are going to apply the following theorem due to Hörmander.

Theorem 2

(Hörmander, Theorem 21.5.3 in [27]) Let $Q\in S_{2n}^{++}({\mathbb {R}})$ be a real symmetric positive matrix of size 2n. There exists a real symplectic matrix $P\in \mathrm {Sp}_{2n}({\mathbb {R}})$ of size 2n such that and some positive numbers $\omega _1,\dots ,\omega _n$ such that

where $D(\omega ) = \mathrm {diag}(\omega _1,\dots ,\omega _n,\omega _1,\dots ,\omega _n)$ is the diagonal matrix such that, for $j=1,\dots ,n$, $D(\omega )_{j,j} = D(\omega )_{j+n,j+n} = \omega _j$ .

Indeed, here, it can be checked that $Q_{(\hbox {QM})}$ (the matrix of $q_{(\hbox {QM})}$) is a symmetric positive matrix (computing, for example, an approximation of its eigenvalues). Thus, applying Theorem 2, we get a symplectic matrix P and some positive numbers $\omega _1<\omega _2<\omega _3$ such that

(55)

Consequently, since P is symplectic, we have

$$\begin{aligned} \exp (2tJQ_{(\hbox {QM})}) = P \exp (2tJD(\omega )) P^{-1}, \end{aligned}$$

where J is the symplectic matrix of ${\mathbb {R}}^{2n}$. Now, applying the monoid morphism (Theorem 3.1 in [12]) introduced also by Hörmander in [26], we get a function $t\mapsto \sigma _t\in \{\pm 1\}$ such that

$$\begin{aligned} \forall t\in {\mathbb {R}}, \ U_t = e^{-itq_{(\hbox {QM})}^w} = \sigma _t V \exp (-it \sum _{j=1}^3 \omega _j (x_j^2 -\partial _{x_j}^2) ) V^* \end{aligned}$$

where $\pm V$ is the Fourier Integral Operator associated with P. Note that V is unitary. Furthermore, by a straighforward argument of continuity we deduce that $\sigma _t=1$ for all $t\in {\mathbb {R}}$. Thus, to conclude, we just have to prove that $(\omega _1,\omega _2,\omega _3) = \frac{\pi }{180}(20,75,132)$.

Step 2: To determine $\omega$. First, we observe that the matrices $JQ_{(\hbox {QM})}$ and $JD(\omega )$ are similar. Indeed, since $P\in \mathrm {Sp}_6({\mathbb {R}})$, we have and applying (55) we deduce that

A fortiori, $JQ_{(\hbox {QM})}$ and $JD(\omega )$ have the same eigenvalues. Thus, the eigenvalues of $JQ_{(\hbox {QM})}$ are

$$\begin{aligned} \sigma (JQ_{(\hbox {QM})}) = \sigma (JD(\omega )) = \{ i \omega _1,-i \omega _1,i\omega _2,-i\omega _2,i\omega _3,-i\omega _3\}. \end{aligned}$$

(56)

Consequently, to determine $\omega$ we just have to determine the roots of the characteristic polynomial of $JQ_{(\hbox {QM})}$, denoted $\chi _{(\hbox {QM})}$. By a straightforward calculation, we observe that

$$\begin{aligned} \left( \frac{3}{\pi }\right) ^{6} \chi _{(\hbox {QM})}(\frac{\pi X}{3}) = X^6 + \frac{\lambda _1+\lambda _2+\lambda _3+3}{2} X^4 + \frac{\lambda _1\lambda _2+ \lambda _1\lambda _3+\lambda _2\lambda _3+9/4}{4} X^2 \ 3\frac{\lambda _1+\lambda _2+\lambda _3}{32} + \frac{\lambda _1\lambda _2+ \lambda _1\lambda _3+\lambda _2\lambda _3}{8} - \frac{\lambda _1\lambda _2\lambda _3}{8}. \end{aligned}$$

But, by construction $\lambda _1< \lambda _2 < \lambda _3$ are the roots of the polynomial

$$\begin{aligned} 7200 X^3 - 72196 X^2 + 222088 X - 216341. \end{aligned}$$

Thus, $\lambda _1+\lambda _2+\lambda _3$, $\lambda _1\lambda _2+ \lambda _1\lambda _3+\lambda _2\lambda _3$ and $\lambda _1\lambda _2\lambda _3$ are some explicit rational numbers and we deduce that

$$\begin{aligned} \left( \frac{3}{\pi }\right) ^{6} \chi _{(\hbox {QM})}(\frac{\pi X}{3}) = X^6+\frac{407}{120}X^4+\frac{123}{80}X^2-\frac{7}{384}. \end{aligned}$$

Finally, we verify by an explicit computation that

$$\begin{aligned} \chi _{(\hbox {QM})}( i\frac{\pi }{9} )= \chi _{(\hbox {QM})}( i\frac{5\pi }{12} ) =\chi _{(\hbox {QM})}(i\frac{11 \pi }{15})=0. \end{aligned}$$

So, we deduce of (56) that $(\omega _1,\omega _2,\omega _3) = \frac{\pi }{180}(20,75,132)$. $\square$