WDVV equations: symbolic computations of Hamiltonian operators

Abstract

We describe software for symbolic computations that we developed in order to find Hamiltonian operators for Witten–Dijkgraaf–Verlinde–Verlinde (WDVV) equations, and verify their compatibility. The computation involves nonlocal (integro-differential) operators, for which specific canonical forms and algorithms have been used.

1 Introduction

The Witten–Dijkgraaf–Verlinde–Verlinde (WDVV) equations, also known as associativity equations, originated in 2D-topological field theory, but they acquired mathematical interest from the work of B. Dubrovin [6], who showed the deep connections between solutions of WDVV equations and integrable systems of PDEs, in particular, bi-Hamiltonian hierarchies.

Recently [25], we proved that, in low dimensions, the WDVV equations themselves are bi-Hamiltonian systems of PDEs. The proof required sophisticated symbolic computational tools that have been only recently made available [4, 5].

In this paper, we would like to describe the algorithms and symbolic software that we developed in order to carry out the computations in [25]. To do that, let us first describe the mathematical problem that we solved. Following [7], in ${\mathbb {R}}^N$ we are to find a function $F=F(t^1,\ldots ,t^N)$ such that

• 1. $\displaystyle \frac{\partial ^3 F}{\partial t^1\partial t^\alpha \partial t^\beta } = \eta _{\alpha \beta }$ is a constant symmetric nondegenerate matrix;

• 2. $c^\gamma _{\alpha \beta } = \eta ^{\gamma \epsilon } \displaystyle \frac{\partial ^3 F}{\partial t^\epsilon \partial t^\alpha \partial t^\beta }$ are the structure constants of an associative algebra;

• 3. F is quasihomogeneous: $F(c^{d_1}t^1,\ldots ,c^{d_N}t^N) = c^{d_F}F(t^1,\ldots ,t^N)$.

The conditions of associativity, or WDVV equations, are the system of PDEs

$$\begin{aligned} \eta ^{\mu \lambda }\frac{\partial ^{3}F}{\partial t^{\lambda }\partial t^{\alpha }\partial t^{\beta }}\frac{\partial ^{3}F}{\partial t^{\nu }\partial t^{\mu }\partial t^{\gamma }}=\eta ^{\mu \lambda }\frac{\partial ^{3}F}{\partial t^{\nu }\partial t^{\alpha }\partial t^{\mu }}\frac{\partial ^{3}F}{\partial t^{\lambda }\partial t^{\beta }\partial t^{\gamma }}. \end{aligned}$$

(1)

Item 1 implies that F can be rewritten as

$$\begin{aligned} F=\frac{1}{6}\eta _{11}(t^1)^3 + \frac{1}{2}\sum _{k>1}\eta _{1k}t^k(t^1)^2 + \frac{1}{2}\sum _{k,s>1} \eta _{sk}t^st^kt^1 + f(t^2,\ldots ,t^N), \end{aligned}$$

(2)

up to a second degree polynomial of the field variables. This implies that the WDVV system is an overdetermined system in one unknown function f of $N-1$ independent variables.

A technique introduced in [17] (see also [10]) allows to rewrite the WDVV system as $N-2$ commuting quasilinear first-order systems of conservation laws of the form

$$\begin{aligned} u^i_t = (V^i({\mathbf {u}}))_x,\qquad i=1,\ldots ,n. \end{aligned}$$

(3)

For example, if $N=3$ and

$$\begin{aligned} \eta _{\alpha \beta }= \delta _{\alpha +\beta ,N+1}= \begin{pmatrix} 0 &{}\quad 0 &{}\quad 1 \\ 0 &{}\quad 1 &{}\quad 0 \\ 1 &{}\quad 0 &{}\quad 0 \end{pmatrix}, \end{aligned}$$

(4)

there is only one WDVV equation:

$$\begin{aligned} f_{ttt}=f_{xxt}^2 - f_{xxx}f_{xtt}. \end{aligned}$$

(5)

Introducing the new coordinates $a=f_{xxx}$, $b=f_{xxt}$, $c=f_{xtt}$ we have the compatibility conditions

$$\begin{aligned} \left\{ \begin{array}{l} a_t = b_x,\\ b_t = c_x,\\ c_t = (b^2 - ac)_x. \end{array} \right. \end{aligned}$$

(6)

The above representation allowed to look for a bi-Hamiltonian formalism. That was found in [9]: the right-hand side of (6) was rewritten in two ways

$$\begin{aligned} u_{t}^{i}=A_{1}^{ij}\mathchoice{\frac{\delta H_2}{\delta u^j}}{\delta H_2/\delta u^j}{\delta H_2/\delta u^j}{\delta H_2/\delta u^j}=A_{2}^{ij}\mathchoice{\frac{\delta H_{1}}{\delta u^j}}{\delta H_{1}/\delta u^j}{\delta H_{1}/\delta u^j}{\delta H_{1}/\delta u^j}, \end{aligned}$$

(7)

where $\mathchoice{\frac{\delta }{\delta u^i}}{\delta /\delta u^i}{\delta /\delta u^i}{\delta /\delta u^i}$ stands for the variational derivative, $H_1$ and $H_2$ are two Hamiltonian densities and $A_{1}$, $A_{2}$ are two compatible local Hamiltonian operators. Here, ‘compatible’ means that the pencil $A_1+\lambda A_2$ is a Hamiltonian operator for every value of the parameter $\lambda$. We recall that the Hamiltonian property for a differential operator is equivalent to differential conditions on its coefficients, to be discussed later. Such a property is crucial for the integrability of the system of PDEs.

Only a few other examples of bi-Hamiltonian WDVV systems were known until recently. In [25] we were able to prove that, for $N=3$, all WDVV systems admit a bi-Hamiltonian formulation, with Hamiltonian operators of the form

$$\begin{aligned} A_1^{ij} =&g^{ij}\partial _x^{} + \Gamma ^{ij}_{k}u^k_x + \alpha V^i_qu^q_x\partial _x^{-1} V^j_pu^p_x \nonumber \\&+\beta \left( V^i_qu^q_x\partial _x^{-1}u^j_x + u^i_x\partial _x^{-1}V^j_qu^q_x \right) +\gamma u^i_x\partial _x^{-1}u^j_x, \end{aligned}$$

(8)

$$\begin{aligned} A_2 =&\partial _x(h^{ij}\partial _x + c^{ij}_k u^k_x)\partial _x, \end{aligned}$$

(9)

with suitable Hamiltonian densities $H_1$ and $H_2$. We observe that $A_1$ is a nonlocal first-order homogeneous Hamiltonian operator of Ferapontov type [8]. That is a member of a more general family of nonlocal operators, namely weakly nonlocal operators, that is discussed in [5]. In this context, ‘nonlocal’ means integro-differential operator (derivatives are taken with respect to the variable x). $A_2$ is a local third-order homogeneous Hamiltonian operator in canonical form [11, 12]. See more on these operators in next Section.

One of the computations that is needed in order to prove the bi-Hamiltonian property is the compatibility of the operators. That is made difficult by the fact that one of the operators is weakly-nonlocal (the Ferapontov operator $A_1$). The computation has been made possible by the recent computer algebra packages developed in [4].

Here we will describe how we used the Reduce computer algebra system, see https://reduce-algebra.sourceforge.io/, and its CDE package in particular (see [14] for a general description and [4] for calculations with weakly-nonlocal operators) in the nontrivial calculations that led to the results presented in [25]. We also checked some of the computations in Maple using the package jacobi.mpl from [4]. Maple was also used in a significant part of the computation in the case $N=5$.

All files that are described in the text can be found in a GitHub repository. The files are also downloadable as a single .zip file, see [26].

2 Computing the first-order Hamiltonian operator $A_1$

Below we will focus on finding a first-order, weakly nonlocal Hamiltonian operator of Ferapontov type.

In the expression (8) $(g^{ij})$ is a non-degenerate matrix of functions of the field variables $(u^i)$, whose inverse $(g_{ij})$ can be interpreted as a covariant 2-tensor. The matrix $V^i_q$ is just the Jacobian of the vector function $(V^i)$ of the fluxes of the system (3), and $\alpha$, $\beta$, $\gamma$ are three constants.

The Hamiltonian property of $A_1$ is equivalent to the following conditions [8]:

• $g^{ij}$ is symmetric;

• $\Gamma ^{i}_{jk} = - g_{js}\Gamma ^{si}_k$ are the Christoffel symbols of $g_{ij}$, regarded as a metric;

• the identities:

$$\begin{aligned} g^{ik}V^j_k&= g^{jk}V^i_k, \end{aligned}$$

(10a)

$$\begin{aligned} \nabla _kV^i_j&= \nabla _j V^i_k, \end{aligned}$$

(10b)

$$\begin{aligned} R^{ij}_{kl}&= \alpha \big ( V^i_k V^j_l - V^i_l V^j_k \big ) \nonumber \\&\quad +\beta \big ( V^i_k \delta ^j_l - V^j_k \delta ^i_l - V^i_l \delta ^j_k + V^j_l \delta ^i_k \big ) +\gamma (\delta ^i_k\delta ^j_l - \delta ^i_l\delta ^j_k) \end{aligned}$$

(10c)

hold, where $\nabla$ denotes the covariant derivative with respect to the Levi–Civita connection of $g_{ij}$ and $R^{ij}_{kl}=g^{is}R^j_{skl}$ is the Riemannian curvature tensor. The above conditions also imply that the operator is a Hamiltonian operator for the system of PDEs (3) [24, 27].

Thus, finding an operator (8) actually reduces to finding the constants $\alpha , \beta$ and $\gamma$ and the metric $g^{ij}$. However, finding the metric in general is by no means an easy computational task. To this end, we shall use a theorem from [2] which states that the metric of a first-order Hamiltonian operator for a non-diagonalizable quasilinear first-order system in $n=3$ unknown functions is proportional to a contraction of the square of the Haantjes tensor:

$$\begin{aligned} g_{ij} = f\, H^\alpha _{i\beta }H^{\beta }_{j\alpha }=f\,H_{ij}, \qquad f=f({\mathbf {u}}). \end{aligned}$$

(11)

We recall that the Haantjes tensor is obtained from $(V^i)$ and its derivatives by means of a straightforward formula (see [25]).

After this simplification we are left with a much less demanding computational problem – we need to find one unknown function $f({\mathbf {u}})$ and the above constants. To this end, we will use a computer algebra system to determine our unknowns.

2.1 The computation

The computation of a nonlocal first-order operator is divided into three smaller steps. Firstly, we will write down the candidate metric $g_{ij} = f\,H_{ij}$ and check (10a). Next we will solve the Eq. (10b) which will yield an explicit form of the metric $g^{ij}$ and as a final step, we find the constants $\alpha , \beta$ and $\gamma$ after plugging the found metric into the condition (10c). We will concentrate on a case from Mokhov-Pavlenko’s classification [18], where the determining matrix $\eta$ has the following form:

$$\begin{aligned} \eta ^4= \begin{pmatrix} 1 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad \lambda &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad \mu \end{pmatrix},\ \lambda ^2=1, \ \mu ^2=1, \end{aligned}$$

(12)

since it is not as computationally demanding and can be solved in full generality with respect to the constants $\lambda ,\mu$. The WDVV equation takes the form:

$$\begin{aligned} \mu f_{xxt}^2-\mu f_{xxx}f_{xtt}+\lambda f_{xtt}^2-\lambda f_{xxt}f_{ttt} -1= 0. \end{aligned}$$

(13)

After introducing the coordinates $u^1=f_{xxx}$, $u^2=f_{xxt}$, $u^3=f_{xtt}$ we obtain the first-order WDVV system in conservative form:

$$\begin{aligned} \begin{aligned}&u^1_t = u^2_x,\\&u^2_t = u^3_x,\\&u^3_t = \left( \frac{\mu ((u^2)^2- u^1u^3)+\lambda (u^3)^2-1}{\lambda u^2}\right) _x. \end{aligned} \end{aligned}$$

(14)

All files used for computation are available at a GitHub repository [26].

2.1.1 The first step

The following computation can be found in the file WDVV-3c-Eta4\dne3_lho2. red. In Reduce, we start by loading the package cde and initialization of our environment by the following sequence:

followed by initialization of the right hand-side of the system:

where the last element is just the WDVV equation. Next, define the velocity matrix $V^i_j$ (and its operator) of our system:

We also need to load a riemann4.red library in order to generate the Nijenhuis and Haantjes tensors by the commands:

And finally the contraction of the square of Haantjes tensor is obtained by:

As a last step we check the symmetry condition $H_{ih}V^h_j = H_{jh}V^h_i$:

and save the result hmet for later computations. We remark that, by construction, hmet is a rational function of the field variables u1, u2, u3 and the parameters lam, mu.

2.1.2 The second step

The following computation can be found in the file WDVV-3c-Eta4\dne3_lho3. red. We, again, start by initializing the environment as in previous step. After constructing the velocity matrix $V^i_j$ and loading the previous result we proceed to defining the metric $g_{ij}$ as a functional factor of the contraction of the square of the Haantjes tensor:

Besides the metric itself, we will need its inverse $g^{ij}$ and Christoffel symbols of the third kind $\Gamma ^{ij}_k$ (which are constructed inside the package riemann3.red by means of Christoffel symbols of the first and second kind):

To find the coefficients of the metric we will use the condition (10b) which is in our case equivalent to:

$$\begin{aligned} V^k_j \Gamma ^{si}_k = V^s_k \Gamma ^{ki}_j. \end{aligned}$$

Hence, we first construct the left and right hand-side respectively:

and we assemble the system:

Using the package crack, a solver of overdetermined systems of PDEs [28, 29], we obtain a solution for our unknown function $f({\mathbf {u}})$ by the following sequence:

Then by a mere substitution, have the explicit form of the metric g:

We again export this result in order to finish the final step - finding the nonlocal part of the first-order operator.

2.1.3 The third step

The following computation can be found in the file WDVV-3c-Eta4\dne3_lho4. red. At last, we need to find the nonlocal tail of the operator $A_1$. We simply need to plug everything we found into the condition (1) and see for which values of constants $\alpha , \beta , \gamma$ it is satisfied.

As usual, we initialize the environment. As an addition to that we load the results from the previous file and then construct the metric and generate all Christoffel symbols $\Gamma ^{ij}_k$:

Let us also define the identity matrix which will represent the Kronecker $\delta ^i_j$:

and finally assemble the condition (10c):

The system is a rational expression that should vanish for any value of the field variables u1, u2, u3. That yields an overdetermined linear system on the coefficients alpha, beta, gamma that can be solved only if alp=mu, bet=0, gam=lam. The substitution:

yields a list of non-zero equations in terms of the parameters mu, lam. However, we should check only for all the possible combinations of the values of constants $\mu , \lambda$:

which finally yields four lists of zeroes.

The metric $g^{ij}$ that defines $A_1$, according with the expression (8), is:

$$\begin{aligned} g^{ij}= \frac{\lambda }{\mu } \begin{pmatrix} -a^2\mu -b^2\lambda -4\mu \lambda &{} -b( a\mu +c\lambda ) &{} -b^2\mu - c^2\lambda - 1 \\ -b( a\mu +c\lambda ) &{} -b^2\mu - c^2\lambda - 1 &{} \frac{c( ac\mu -2b^2\mu -c^2\lambda + 1)}{b} \\ -b^2\mu - c^2\lambda - 1 &{} \frac{c(ac\mu -2b^2\mu -c^2\lambda + 1)}{b} &{} \frac{\delta }{b^2} \end{pmatrix}, \end{aligned}$$

(15)

where $\delta = 2ab^2c\lambda -a^2c^2\lambda +2ac^3\mu -2ac\mu \lambda -b^4\lambda -3b^2c^2\mu -2b^2\mu \lambda -c^4\lambda +2c^2 - \lambda$ and the values of constants in (8) are $\alpha =\mu ,\beta =0, \gamma =\lambda$, with $\lambda ,\mu = \pm 1$.

3 Computing the third-order Hamiltonian operator $A_2$

To find a bi-Hamiltonian structure for the system (3) we need a second Hamiltonian operator compatible with $A_1$. It was proved in [25] that in the case $N=3$, 4 and 5 the WDVV systems admit a third-order homogeneous Hamiltonian operator in the canonical form

$$\begin{aligned} A_2 = \partial _x(h^{ij}\partial _x + c^{ij}_ku^k_x)\partial _x \end{aligned}$$

(16)

(we always require $\det \,(h^{ij})\ne 0$; $(h_{ij})$ denotes the inverse matrix). We recall that the Hamiltonian property of $A_2$ is equivalent to the conditions

$$\begin{aligned}&c_{ijk}=\frac{1}{3}(h_{ik,j} - h_{ij,k}). \end{aligned}$$

(17a)

$$\begin{aligned}&h_{mk,s} + h_{ks,m} + h_{ms,k}=0, \end{aligned}$$

(17b)

$$\begin{aligned}&c_{msk,l}= - h^{pq}c_{pml}c_{qsk}. \end{aligned}$$

(17c)

and that the Eq. (17b) implies that $h_{ij}$ is a Monge metric [11, 12]. For computational purposes, that implies that the coefficients $h_{ij}$ are quadratic polynomials of the field variables.

Then, it is also known [13] that the conditions under which a third-order Hamiltonian operator is the Hamiltonian operator for a quasi-linear system of first-order conservation laws are:

$$\begin{aligned}&h_{im}V^{m}_{j}=h_{jm}V^m_{i},\end{aligned}$$

(18a)

$$\begin{aligned}&c_{mkl}V^m_{i}+c_{mik}V^m_{l}+c_{mli}V^m_{k}=0,\end{aligned}$$

(18b)

$$\begin{aligned}&h_{ks}V^k_{ij}=c_{smj}V^m_{i} + c_{smi}V^{m}_{j}. \end{aligned}$$

(18c)

Thus, with respect to the problem of finding a third-order operator for WDVV systems, knowing the vector of fluxes $(V^i({\mathbf {u}}))$ reduces the above system of equations to a linear system of algebraic equations in the coefficients of the quadratic polynomials $h_{ij}$ that can be solved in short time for small values of N.

For the particular example that we discussed in the previous section (system (14)) we can easily solve the system (18) and verify that the unique solution $h_{ij}$ that we get fulfills the conditions (17b) and (17c). This implies that $h_{ij}$ generates a third-order homogeneous Hamiltonian operator (16) using the formula (17a). We have the expression:

$$\begin{aligned} h_{ij}= \begin{pmatrix} b^2+\mu &{}\quad b\mu (\lambda c-\mu a) &{}\quad -\mu \lambda b^2 \\ b\mu (\lambda c-\mu a) &{}\quad \lambda +a^2-\lambda c(2\mu a-\lambda c) &{}\quad \lambda b(\mu a-\lambda c) \\ -\mu \lambda b^2 &{}\quad \lambda b(\mu a-\lambda c)) &{}\quad b^2 \\ \end{pmatrix}. \end{aligned}$$

(19)

As this computation is not using any special package, there is no need to go into details in the code itself. The computation is provided in the file WDVV-3c-Eta4 \wdvv_3ord_op_eta4.red

4 Compatibility of the Hamiltonian operators $A_1$, $A_2$

A bi-Hamiltonian system is just a system of PDEs which is Hamiltonian with respect to two Hamiltonian operators, $A_1$ and $A_2$ in our case, which are compatible. The compatibility condition can be either checked by means of the requirement that the pencil $A_1+\lambda A_2$ is a Hamiltonian operator for every value of the parameter, or, equivalently, by computing the Schouten bracket $[A_1,A_2]$ (see [5]).

However, tensorial conditions which are equivalent to the compatibility condition of $A_1$ and $A_2$ as in our case are not available. So, the calculation should be made through the definition, which makes a pen-and-paper approach to the problem almost completely unfeasible. This is also due to the fact that an algorithm to do the calculation in the nonlocal case was missing until recently. A computational algorithm [5] and software packages [4] became available a short time ago. Through this software we checked the compatibility of the operators $A_1$ and $A_2$.

All of the computations described below are contained in the files wdvv_comp_ nl1_eta4.red and wdvv_comp_nl2_eta4.red. For convenience they are split into two separate files where in the first one we check the bracket [$A_1,A_1$] and then in the second file we check the compatibility condition [$A_1,A_2$]. Let us focus on the second computation.

We start as usual by initializing the environment:

From now on, we initialize the Ferapontov operator $A_1$. The metric of the first-order operator is loaded from another file:

To construct the local part we further need to initialize the Christoffel symbols $\Gamma ^{ij}_k$ in the same way as before. We construct the local part of the operator by the following lines:

where b represents the contraction $\Gamma ^{ij}_k({\mathbf {u}})u^k_x$.

The other, nonlocal, part of our operator consists of a nonlocal tail which is defined by the previously found values of constants $\alpha , \beta , \gamma$:

Here the constants are represented by a symmetric matrix $c^{ij}$; $w^k_i$ on the other hand represents the nonlocal terms.

Now, we import the metric (in lower indices) of the third-order operator, stored in gl3, and we generate the constants $c_{ijk}= \frac{1}{3}(h_{ki,j} - h_{ji,k})$:

then we raise two indices:

and finally we need the contraction $c^{ij}_ku^k_x$:

We assemble the local part:

and even though there is no nonlocal part, we need to introduce a zero tail:

and put the two parts together:

All that remains to do is to configure the nonlocal variables. Here we have two distinct operators so we need three nonlocal variables and their names will be different this time:

Finally we prepare the jet space

and run the computation:

We have to plug in all the possible combinations of values $\mu , \lambda = \pm 1$ to obtain the desired result:

the result is always zero.

Thus we have proved by a direct computation that previously found operators $A_1$ and $A_2$ are compatible Hamiltonian operators.

5 A computation in Maple using jacobi

In this section we describe another tool which might be used to check if a weakly nonlocal operator is Hamiltonian and possibly compatible with another one: the package jacobi [4]. Since that package is written in Maple some readers might find it more convenient, and we thought that it is worth to explain how to achieve the same result with a different computer algebra system.

The main structural difference with the computation in the previous section is that jacobi uses the language of distributions, which is more popular among Theoretical Physicists.

The plan is the same as before, we need the metrics of both operators to generate their local parts and the three constants for the nonlocal part. However the notation is slightly different and that will be the main point of this section.

For the following computation we also need the external Maple library jacobi.mpl which can be found on [23]. The computation itself is contained in the fileWDVV-3c-Compat_Example_Eta4.mw. For simplicity, we consider $\lambda =\mu =1$ for this case.

We start by constructing $A_1$. The metric is declared as follows:

where by u[1, x, 0] we denote in our case $u^1_{0x}$ (with the convention $u_{0x} = u)$. Other coefficients are defined in a similar way. After the matrix, we need the Christoffel symbols which can be computed in Maple or, if we have them available from previous computations, input them as a three-dimensional field:

and so on. Next, we input the two tail vectors W1X and W1Y which in this notation differ only by using a different symbol for the independent variable. The W1X is introduced as follows:

and the same goes for W1Y with y instead of x. The first-order nonlocal operator $A_1$ is then assembled by:

Next for the third-order operator we need its metric $h_{ij}$ as well as constants $c^{ij}_k$. We input the metric, denoted by g3, as before and, similarly to Christoffel symbols for the previous metric, the constants $c^{ij}_k$ are handled as a three-dimensional field:

and so on. For the local operator the tail matrix is just a zero matrix:

Finally, the operator $A_3$ is assembled by:

We can then proceed to computing the Schouten bracket itself. We load the package jacobi by:

The bracket is calculated by the following command and saved into an array T3:

We can view the result just by printing the contents of T3 (further simplification might be needed):

The same command is used for computing the Schouten bracket of $[A_1, A_3]$, proving the compatibility of the two operators:

6 Large scale computations

While doing the above computations for $N=3$ is not a time demanding task, the situation quickly changes for $N=4$, where just finding the third-order operator usually takes around 5 h on an average PC. To obtain the results for $N=5$ in a similar way one would have to have an access to a huge amount of computational resources.

In this section, we describe the way in which we found the third-order operators in the case $N=5$ (results in [25]). Since we did not even try to find a first-order operator when $N=5$, we will not deal with related problems. Only partial results are known in the case $N=4$ [20].

There are two ways how to approach this problem. First, we could find a powerful enough machine to handle this computation. However, we are talking about hundreds of GB of RAM memory and possibly a week of computing time.

The other option is to optimize the computation itself, so that it not excessively resource-demanding. It can be done in a few ways, but the most effective at the beginning is to reduce the number of equations we try to solve at once.

As both approaches were used to obtain the results for $N=5$ let us explore those below in more detail. All files used below are in a separate folder \WDVV-5c-Large_scale.

6.1 Using a supercomputer

If you have an access to a supercomputer it could be beneficial to try to use a rather “brute force” approach. Given the fact we cannot say for sure how many equations need to be solved in order to obtain the solution, it is by far the best first try.

To obtain the result in case of $N=5$ for a canonical choice of $\eta$ we have used the whole super-computing grid provided by a virtual organization CESNET – MetaCentrum. The command to start the computation on the grid can be found in the file starter.sh.

It is also important to mention that we need to provide the file with the generated system in the conservative form as the algorithm has not been transferred into Maple yet. The algorithm itself is contained inside the file w10_hydro_system_gen.red whose output file (w10_eta2_eq.red) is then loaded into Maple. Finding the quasi-linear system in the conservative form using an algorithm was described in detail in [25] and is an interesting computational problem. However, since the program literally follows the steps described there, we will omit repeating the details. Note that some syntax adjustments are required before importing.

After we obtain the result from a remote computation it is important to have an independent check. We will use Reduce, file w10_ham1_eta2_check.red. We initialize as usual:

where in the files w10_eta2_eq.red and w10_eta_gl3.red we saved the general WDVV system and the result of the computation—the metric $h_{ij}$ respectively which again needs to be transformed into the proper syntax.

Then we need to confirm that the Eq. (18) holds for the metric $h_{ij}$ already in place (in form of the constants $c_{ijk}$) to confirm the compatibility with the first-order system. Lastly we need to check the Hamiltonian condition which we assemble by:

which also passes in our case and confirms the result obtained by a super-computation.

6.2 Simplifying the computation

If we look into the algorithm we can quickly realize that it is massively inefficient in the sense of computational resources. The unknowns are exactly the coefficients of the quadratic polynomials in the metric of the third-order operator. In the case of $N=5$ we need to find 1540 unknown constants.

On the other hand, the conditions (18) generate 2100 equations but those equations are polynomial with respect to the field variables. We need to collect the coefficients of each monomial out. There will be an immense amount of those collected equations. Just for comparison, in the case $N=4$ we have 231 unknown coefficients, slightly more than 450 equations polynomial in the field variables and we generate around 350,000 of such linear equations.

It is now obvious that solving all of the collected equations is completely unnecessary and in this case a waste of resources. We could try to make an algorithm that will solve only a few of the equations which are polynomial in the field variables at once, plug the result into the rest and repeat the process batch by batch until all the equations are solved. This process could work, however, it is very hard to predict how much resources it will consume and what is the optimal size of the mentioned batch. The two extreme cases being the batch is too small and the number of iterations will project into the time of the computation or the batch being too big literally clogging the memory and maybe never finish.

There is certainly enough space to optimize the whole process of solving those massively over-determined system of linear equations. However, by experiment we have found out that in our case there are around 60 equations out of the original 2100 polynomial ones needed to completely solve the posed problem. They originate from the system (20) below.

Also, saving partial results in the whole process is an absolute must. While the setup of the problem does not consume much time, relative to solving it, it is much easier to plug in a partial result from earlier and continue the computation.

More specifically, instead of solving the whole set of equations (18) we can simply restrict the system to

$$\begin{aligned} h_{ik}V^{k}_{j}=h_{jk}V^k_{i}, \qquad i,j,k=1,\dots ,4. \end{aligned}$$

(20)

This computation takes around 8 hours on an average PC and the result is checked in the same way as in the previous case.

7 Conclusions

In this paper we described the computations that we performed in [25] in order to support the conjecture that all WDVV equations are endowed with bi-Hamiltonian formalism. These are symbolic software calculations that are nontrivial under both the viewpoint of the mathematical algorithm and the viewpoint of the symbolic programming that is needed to achieve the goal.

However, we think that further computational goals might be achieved to support and broaden the above conjecture.

• First of all, in the $N=4$ case the bi-Hamiltonian property has been proved only for one of the two canonical forms of $\eta$, $\eta ^{(1)}$. This was a result found in [20]. For the system arising from the second canonical form $\eta ^{(2)}$ we found a third-order operator but we haven’t even tried to find a first-order operator. It should be said that the method that we used in the case $N=3$, coming from a characterization in [2], is not applicable in higher dimensions, hence the calculation might be a real challenge. Of course, the case $N=5$ is even harder. A way to find a first-order operator might be to restrict the search to matrices $g^{ij}$ of rational functions of a certain degree. This common pattern can be observed in our examples (see also [25]) and might greatly simplify the search. More generally, it should be possible to program a cde or jacobi computation of all first-order nonlocal operators that are compatible with a given third-order operator. This is a task for a future research project.

• Then, there is a well-known generalization of the WDVV equations, the oriented associativity equation [15]. These equations and the associated geometric structures are called F-manifolds with compatible flat structures [16], or simply flat F-manifolds. They share many properties with WDVV equations and Dubrovin-Frobenius manifolds including the existence of an associated integrable dispersive hierarchy (see [1] for details). We observe, in particular, that the oriented associativity equation has an infinite hierarchy of nonlocal symmetries [22], a first-order local Hamiltonian operator of the same type as $A_1$ [19] and a third-order nonlocal Hamiltonian operator which is the straightforward generalization of $A_2$ [3, 21]. Until now, the results on Hamiltonian operators are known only for one of the simplest cases of oriented associativity equation; more calculations are needed in order to support a conjecture on the bi-Hamiltonianity of the F-manifold equation. Indeed, the compatibility of $A_1$ and $A_2$ is still an open question even in the simplest case, see [21] for a discussion.

The above problems are on our schedule.