As discussed in Section 4, in order to prove optimality of the 600-cell, it suffices to find an Hermite interpolating polynomial for the kernel function, agreeing with it at the scalar products of the 600-cell, and show that it is positive definite. It is done here by showing that the Jacobi coefficients of the polynomial $h$ are positive.

Instead of minimizing the p-frame energy on $\mathbb S^3$, we symmetrize measures under consideration and minimize $$ \iint_{\mathbb R\mathbb P^3} \left(\frac{\tau(x,y) +1}2 \right)^{p/2} d\mu(x)d\mu(y) $$

Recall that here $\tau(x,y) = \cos(\theta(x,y))$, and $\theta$ denotes the geodesic distance renormalized to $[0,\pi]$.

In [1]:

```
a = 1/2
b = -1/2
scalar_prods = [1,(sqrt(5)-1)/4, -1/2, -(sqrt(5)+1)/4, -1]
```

In [2]:

```
A = [var('h_%d'%i) for i in (0..8)]
p,t = var('p t')
```

$h$ is a polynomial spanned by Jacobi poynomials $ C_n^{(\frac12,-\frac12)}(t), $ $ n=0,\ldots,8; n\neq 6 $: $$ h = \left(\sum_{n=0}^5 + \sum_{n=7}^8\right) h_n\, C_n^{\left(\frac12,-\frac12\right)}(t). $$

In [3]:

```
h = symbolic_expression(0)
for i in (0..5):
h += A[i]*jacobi_P(i,a,b,t)
for i in (7..8):
h += A[i]*jacobi_P(i,a,b,t)
```

Derivative of h:

In [4]:

```
hprime = diff(h,t)
```

Let $f$ be the kernel of the symmetrized problem

In [5]:

```
f = ((t+1)/2)^(p/2)
fprime = diff(f, t)
```

Jacobi coefficients of the interpolating polynomial are found from the Hermite interpolation conditions. f and h must agree at all the scalar products, and their derivatives must be equal at all the inner products except the endpoints $\pm1$:

In [6]:

```
inter0 = [f.subs(t=s)==h.subs(t=s) for s in scalar_prods]
inter1 = [fprime.subs(t=s)==hprime.subs(t=s) for s in scalar_prods[1:-1]]
interpolate = inter0 + inter1
```

In [7]:

```
coeffs_sol = solve(interpolate, A)[0]
```

Coefficients $h_0,\ldots h_8$ as functions of $p$:

In [8]:

```
coeffs = [c.rhs() for c in coeffs_sol]
```

We shall carry out all the non-symbolic computations in the interval arithmetic. First, we set the format for the output of a computation as an interval.

In [9]:

```
sage.rings.real_mpfi.printing_style = 'brackets'
```

We shall need a bound on the absolute value of the derivative of $h_i$. It is obtained by expanding the expression for $h_i$ into a sum, then replacing every term in the sum by its maximal absolute value on $[8,10]$. Finally, all the absolute values are summed up using triangle inequality.

Summands of $\frac{dh_i}{dp}$ are easy to estimate by monotonicity. After the *expand* command, derivative cprime is a sum of several summands, each of which is a product of factors, monotonic on $[8,10]$. We exploit this structure by replacing every (positive) factor with its maximal value.
deriv_bound is obtained by summing up the absolute values of the operands of cprime.

In [10]:

```
deriv_bounds = []
for c in coeffs:
cprime = diff(expand(c),p)
csub = cprime.subs((1/8*sqrt(5) + 3/8)^(1/2*p)==RIF(1/8*sqrt(5) + 3/8)^(1/2*RIF(8))
).subs(4^(1/2*p)==4^(1/2*RIF(8))
).subs((-1/8*sqrt(5) + 3/8)^(1/2*p)==RIF(-1/8*sqrt(5) + 3/8)^(1/2*RIF(8))
).subs(p=RIF(10))
deriv_bound = symbolic_expression(0)
for o in csub.operands():
deriv_bound += abs(o)
deriv_bounds.append(RIF(deriv_bound))
```

In [11]:

```
deriv_bounds
```

Out[11]:

[[0.0087005417083128467 .. 0.0087005417083128728], [0.046057770601501762 .. 0.046057770601502013], [0.55126143476111155 .. 0.55126143476111378], [0.35931366896653027 .. 0.35931366896653167], [0.11947185504173116 .. 0.11947185504173164], [0.049788504505656708 .. 0.049788504505657035], [0.0000000000000000 .. -0.0000000000000000], [0.11682913291810448 .. 0.11682913291810493], [0.57379946418760030 .. 0.57379946418760253]]

It follows that derivatives of all the coefficients are uniformly bounded on $[8,10]$ by e.g. 0.6.

Notice that, since we are using the default precision of 53 bits, here and in what follows the last 3-4 digits are not significant. This can be also seen by switching to the default printing style using

`sage.rings.real_mpfi.printing_style = 'question'`

Check that the coefficients $h_0,\ldots,h_4$ are positive at $p=8$:

In [12]:

```
[RIF(coeffs[i].subs(p=RIF(8))) for i in (0..4)]
```

Out[12]:

[[0.054687499999999861 .. 0.054687500000000174], [0.21874999999999936 .. 0.21875000000000092], [0.20833333333332870 .. 0.20833333333333793], [0.087499999999996913 .. 0.087500000000003131], [0.014285714285713013 .. 0.014285714285715512]]

The last three $h_5,h_7,h_8$ are equal to zero at $p=8$, and so will requre somputing the second derivative.

In [13]:

```
[RIF(coeffs[i].subs(p=RIF(8))) for i in [5,7,8]]
```

Out[13]:

[[-1.8041124150158794e-16 .. 1.7347234759768071e-16], [-1.1379786002407855e-15 .. 1.0824674490095277e-15], [-4.3298697960381106e-15 .. 4.3298697960381106e-15]]

For the first set of coefficients we proceed as follows. Since the first derivative of $h_n$ is bounded by $0.6$, if for some $p_0$, $h_n(p_0) > 0$, then the same applies to $h_n(p_0 + h_n(p_0))$, and also $$ h_n(p) >0, \qquad p \in [p_0, p_0 + h_n(p_0)]. $$ The loop below iterates this argument, and stops once p has reached the value 10.

In [14]:

```
for n in range(5):
c = coeffs[n]
p_it = RIF(8)
numit = 0
while p_it <= RIF(10):
p_it = p_it + RIF(c.subs(p = p_it))
numit += 1
# Since we got outside the loop, our coefficient must be positive.
print("Coefficient h_%d:" % n)
show(expand(c))
print("is positive for p in [8,10].\n\n")
# print(numit)
```

Coefficient h_0:

is positive for p in [8,10]. Coefficient h_1:

is positive for p in [8,10]. Coefficient h_2:

is positive for p in [8,10]. Coefficient h_3:

is positive for p in [8,10]. Coefficient h_4:

is positive for p in [8,10].

Symbolic verification that the last three coefficients turn to 0 at $p=8$:

In [15]:

```
[coeffs[i].subs(p=8).expand() for i in [5,7,8]]
```

Out[15]:

[0, 0, 0]

The same applies to $h_7$ and $h_8$ at $p=10$:

In [16]:

```
[coeffs[i].subs(p=10).expand() for i in [7,8]]
```

Out[16]:

[0, 0]

For $h_5$ it will suffice to verify that $dh_5/dp > 0 $ on $[8,10]$. This will be done similarly to the verification of positivity of the coefficients above, using the second derivative $ \frac{d^2\,h_5}{dp^2} $. For the same reason we shall need 2nd derivatives of $h_7, h_8$, so we compute them as well.

In [17]:

```
coeffs2 = [coeffs[5]] + coeffs[7:9]
```

In [18]:

```
deriv2_bounds = []
for c in coeffs2:
c2prime = diff(diff(expand(c),p),p)
c2sub = c2prime.subs((1/8*sqrt(5) + 3/8)^(1/2*p)==(1/8*sqrt(5) + 3/8)^(1/2*RIF(8))
).subs((-1/8*sqrt(5) + 3/8)^(1/2*p)==(-1/8*sqrt(5) + 3/8)^(1/2*RIF(8))
).subs(4^(1/2*p)==4^(1/2*RIF(8))
).subs(p=RIF(10))
deriv2_bound = symbolic_expression(0)
for o in c2sub.operands():
deriv2_bound += abs(o)
deriv2_bounds.append(RIF(deriv2_bound))
```

In [19]:

```
deriv2_bounds
```

Out[19]:

[[0.017230522423253111 .. 0.017230522423253195], [0.034153895491620095 .. 0.034153895491620262], [0.16177123014881653 .. 0.16177123014881740]]

Hence a uniform bound for $h_n'', n = 5,7,8$ is e.g. 0.2. Just as above, if for some $p_0$, $h_n'(p_0) > 0$, then the same applies to $h_n'(p_0 + 5h_n(p_0))$, and also $$ h_n(p) >0, \qquad p \in [p_0, p_0 + 5h_n(p_0)]. $$ We iterate this argument for $h_5$.

In [20]:

```
cprime = diff(expand(coeffs[5]),p)
p_it = RIF(8)
numit = 0
while p_it <= RIF(10):
p_it = p_it + RIF(5*cprime.subs(p = p_it))
numit += 1
# print(numit)
print("Derivative of the coefficient h_%d:" % 5)
show(cprime)
print("is positive for p in [8,10].\n\n")
```

Derivative of the coefficient h_5:

is positive for p in [8,10].

This gives the desired positivity of $h_5$ on $[8,10]$. To verify positivity of $h_n, n=7,8$, we show

- positivity of $h_n'$ on $[8,8.5]$;
- negativity of $h_n'$ on $[9.5, 10]$;
- positivity of $h_n$ on $[8.5,9.5]$.

In [21]:

```
R = RealIntervalField(100)
```

Instead of performing steps with variable length as we did above, we shall make the step size fixed, in order to avoid accumulation of error. To prevent making steps that are too long, the fixed step is compared to $5h_n'(p_0)$. We also increase the precision of interval arithmetic from the default 53 bits to 100 bits.

In [22]:

```
for n in (7..8):
cprime = diff(expand(coeffs[n]),p)
p_it = R(8)
step = R(0.00002)
numit = 0
flag = True
while p_it <= R(8.5):
if (R(cprime.subs(p = p_it)) > step):
p_it = p_it + 5*step
numit += 1
else:
flag = False
break
if (flag):
print("Derivative of the coefficient h_%d:" % n)
show(cprime)
print("is positive for p in [8,8.5].\n\n")
```

Derivative of the coefficient h_7:

is positive for p in [8,8.5]. Derivative of the coefficient h_8:

is positive for p in [8,8.5].

Negativity of $h_n'(p)$, for $p$ in $[9.5,10]$:

In [23]:

```
for n in (7..8):
cprime = diff(expand(coeffs[n]),p)
p_it = R(9.5)
step = R(0.00001)
numit = 0
flag = True
while p_it <= R(10):
if (R(cprime.subs(p = p_it)) < -step):
p_it = p_it + 5*step
numit += 1
else:
flag = False
break
if (flag):
print("Derivative of the coefficient h_%d:" % n)
show(cprime)
print("is negative for p in [9.5,10].\n\n")
```

Derivative of the coefficient h_7:

is negative for p in [9.5,10]. Derivative of the coefficient h_8:

is negative for p in [9.5,10].

It remains to justify positivity of $h_n$ itself in the interval $[8.5,9.5]$. Since $h_n$ is positive at the endpoints of this interval, the strategy used for the first five coefficients applies here as well.

In [24]:

```
for n in (7..8):
c = expand(coeffs[n])
p_it = R(8.5)
step = R(0.00001)
numit = 0
flag = True
while p_it <= R(9.5):
if (R(c.subs(p = p_it)) > step):
p_it = p_it + step
numit += 1
else:
flag = False
break
if(flag):
print("Coefficient h_%d:" % n)
show(c)
print("is positive for p in [8.5,9.5].\n\n")
```

Coefficient h_7:

is positive for p in [8.5,9.5]. Coefficient h_8:

is positive for p in [8.5,9.5].