Inconsitency in matrix multiplication

Hi All,

I was playing with the Quantum optics package and I saw something weird
So the resulting matrix is inconsistent in comparison with other CAS systems (Mathematica)
Mathematica:

 # Spin operators 

Sx   =  1/Sqrt[2] { {  0, 1, 0} , { 1, 0, 1}, {0, 1, 0 } };
Sy   =   I/Sqrt[2]{{0, -1, 0}, {1, 0, -1}, {0, 1, 0}};
Sz   =   { {1, 0, 0 }, {0, 0, 0}, {0, 0, -1}};
eye =   IdentityMatrix[3];

(* Electron triplet operators *)

eSx = KroneckerProduct[Sx, eye];
eSy = KroneckerProduct[Sy, eye];
eSz = KroneckerProduct[Sz, eye];

 (* Nuclear operators *)
 
Ix = KroneckerProduct[eye, Sx];
Iy = KroneckerProduct[eye, Sy];
Iz = KroneckerProduct[eye, Sz];
 
(* Total angular momentum *)
 Jx = eSx + Ix;
Jy = eSy + Iy;
Jz = eSz + Iz;

(*J^2*)
Jx*Jx+Jy*Jy+Jz*Jz

(*gives me *)
(4	0	0	0	0	0	0	0	0
0	1	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0
0	0	0	1	0	0	0	0	0
0	0	0	0	0	0	0	0	0
0	0	0	0	0	1	0	0	0
0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	1	0
0	0	0	0	0	0	0	0	4
)

Julia:

using TensorOperations
using LinearAlgebra

Sx=1/√2 * [ [0 1 0]
            [1 0 1] 
            [0 1 0]]

Sy=1/√2 * [ [0 -1im 0]
            [1im 0 -1im] 
            [0 1im 0]]

Sz= [   [1 0 0]
        [0 0 0] 
        [0 0 -1]]

S=[Sx,Sy,Sz]
eye=[Matrix(1*I,3,3) ,Matrix(1*I,3,3) ,Matrix(1*I,3,3) ]
eS=kron.(S,eye)
IS=kron.(eye,S)
JS=eS+IS
real.(JS'*JS)

# This gives 

9×9 Matrix{Float64}:
 6.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  4.0  0.0  2.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  2.0  0.0  2.0  0.0  0.0  0.0  0.0
 0.0  2.0  0.0  4.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  2.0  0.0  4.0  0.0  2.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  4.0  0.0  2.0  0.0
 0.0  0.0  0.0  0.0  2.0  0.0  2.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  2.0  0.0  4.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  6.0

where is Julia code that results different stuff?

sorry was editing it noticed it later my bad

This:

is not equivalent to the Mathematica code. In the Julia code, JS'*JS is equal to

JS[1]'*JS[1] + JS[2]'*JS[2] + JS[3]'*JS[3]

where ' denotes the adjoint and * is matrix multiplication. I don’t know much about Mathematica, but I assume that Jx*Jx+Jy*Jy+Jz*Jz in Mathematica uses elementwise multiplication and no adjoint.

If you use elementwise multiplication and no adjoint in Julia, you get the same result as in Mathematica:

julia> real(JS[1] .* JS[1] + JS[2].*JS[2] + JS[3].*JS[3])
9×9 Matrix{Float64}:
 4.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  4.0

Thank you for the clarification.
my bad …