Let $d(a,b) = 1 - \frac{2\gcd(a,b)^3}{ab(a+b)}$ be a metric on natural numbers without $0$.

The metric space $X = \{x_0,x_1,\cdots,x_n\},n>2$ is isometric embeddable in $\mathbb{R}^n$ if and only if the matrix: $$M(x_0,x_1,\cdots,x_n) = (1/2 (d(x_0,x_i)^2+d(x_0,x_j)^2-d(x_i,x_j)^2))_{1 \le i,j \le n}$$ is positive semidefinite.

So my question is:

Is the matrix above for $d$ as above positive semidefinite for all choices of $x_i \in \mathbb{N}$? (Maybe it is possible to prove this using quadratic forms and then transform it to $\sum_{i} a_{ii} y_i^2$ showing then that $a_{ii}\ge 0$?

If it is so, then this would one allow to do euclidean geometry of natural numbers. For instance for three (pairwise distinct) points / natural numbers we would have:

1. a triangle

2. law of sines

3. law of cosines

4. All other theorems concerning triangles

Then in the limit three consecutive numbers / primes would build an equilateral triangle of side length $1$. Hence one could imagine primes ("in the limit") as an infinite dimensional simplex, which would be a funny thing to think about.

Thanks for your help.

Related question: https://math.stackexchange.com/questions/3385102/is-this-metric-matrix-positive-semidefinite

See Theorem 2.4 in https://books.google.de/books?id=7_DuCAAAQBAJ&printsec=frontcover&hl=de&source=gbs_ge_summary_r&cad=0#v=onepage&q&f=false for isometrically embedding of $(\mathbb{N},d)$ in a Hilbert space.

Edit: Here is some Sage code in case one wants to check this numerically for some examples:

def dABC(a,b):
    """ABC"""
    return 1- 2*gcd(a,b)**3/(a*b*(a+b))

def MM(xx,d=dABC):
    N = len(xx)
    return matrix([[1/2*(d(xx[0],xx[i])**2+d(xx[0],xx[j])**2-d(xx[i],xx[j])**2) for i in range(1,N)] for j in range(1,N)])

def skp(a,b,d=dABC):
    return 1/2*(d(a,1)**2+d(b,1)**2-d(a,b)**2)

def schur(M):
    from scipy.linalg import schur
    import numpy as np
    M_np = np.matrix(M,dtype='float64')
    A,B = schur(M_np,output="complex")
    return (matrix(np.asmatrix(A)),matrix(np.asmatrix(B)))

def createEmbedding(rr):
    M = MM(rr)
    n = len(rr)+1
    A,B = schur(M)
    E = diagonal_matrix([sqrt(x) for x in A.diagonal()])
    X = B*E
    ee = [ matrix([[i==j] for i in range(1,n-1)],ring=QQ) for j in range(1,n-1)]
    #print ee
    xx = [ X.transpose()*ee[i] for i in range(n-2)]
    return xx

N = 20
for i in primes(N):
    for j in primes(i+1,N):
        for k in primes(j+1,N):
            a = dABC(i,j)
            b = dABC(j,k)
            c = dABC(k,i)
            s = 1/2*(a+b+c)
            area = sqrt(s*(s-a)*(s-b)*(s-c)).n()
            alpha = pi.n()-arccos((skp(j,k)-skp(j,i)-skp(k,k)+skp(k,i))/(b*c))
            beta = pi.n()-arccos((skp(j,i)-skp(k,j)-skp(i,i)+skp(i,k))/(a*c))
            gamma = pi.n()-arccos((skp(j,k)-skp(k,i)-skp(j,j)+skp(j,i))/(b*a))
            print i,j,k,"area:",area, "sum:",(alpha+gamma+beta).n(),pi.n()
            print i,j,k,"sine law:",a/sin(alpha).n(),b/sin(beta).n(),c/sin(gamma).n()
            print i,j,k,"lengths:", a.n(),b.n(),c.n()
            print i,j,k,"cosine law: c", c**2.0,(a**2+b**2-2*a*b*cos(gamma)).n(),cos(gamma).n()
            print i,j,k,"cosine law: b", b**2.0,(c**2+a**2-2*c*a*cos(beta)).n(),cos(beta).n()
            print i,j,k,"cosine law: a", a**2.0,(c**2+b**2-2*c*b*cos(alpha)).n(),cos(alpha).n()
for n in range(2,101):
    print n, MM(range(1,n)).is_positive_definite()

Second Edit: Just out of curiosity: For $(a,b,c)=(1,2,2k+1)$, so $c \ge 3$ is odd, we get using the sum of angles in a triangle:

$$\alpha + \beta + \gamma = \pi$$

the following curious identity. For each odd $c \ge 3$ we have:

$$\operatorname{acos}(\frac{4 \, c^{5} + 28 \, c^{4} + 62 \, c^{3} + 2 \, c^{2} - 153 \, c - 135}{12 \, {\left(c + 2\right)}^{3} {\left(c + 1\right)} c} ) +$$ $$ \operatorname{acos}(\frac{14 \, c^{5} + 98 \, c^{4} + 226 \, c^{3} + 142 \, c^{2} - 135 \, c - 153}{18 \, {\left(c^{2} + 2 \, c - 1\right)} {\left(c + 2\right)}^{2} {\left(c + 1\right)}}) + $$ $$\operatorname{acos}(\frac{4 \, c^{6} + 24 \, c^{5} + 70 \, c^{4} + 156 \, c^{3} + 187 \, c^{2} - 18 \, c - 135}{12 \, {\left(c^{2} + 2 \, c - 1\right)} {\left(c + 2\right)} {\left(c + 1\right)}^{2} c}) = \pi$$

Third edit:

I think the main property which distinguishes $d$ for example from the Jaccard or other metrics is the proven property ( https://mathoverflow.net/a/342921/6671) :

For all $a \neq b, a\neq c$ we have:

$$d(a,b)+d(a,c) > 1$$

I have tested other metrics with this property and they also seem to embedd in Euclidean Space. On the other hand metrics who do not have this property do not seem to embedd. So I think this is the point to be taken into consideration.

If someone has an idea how to exploit this property that would be very nice!