About the orthogonality of columns in Zadoff-Chu matrix As we know that Zadoff-Chu matrix is similar to Walsh-Hadamard matrix where every columns in those matrices is orthogonal with the any other column.

For Walsh-Hadmard matrix, the orthogonality is clear and can be demonstrated easily. But regarding the Zadoff-Chu matrix, I don't see that orthogonality is right as shown below: 

a = (sqrt(2) + 1j*sqrt(2))/2;
W = [a 1 a -1; -1 a 1 a ; a -1 a 1; 1 a -1 a];  %%Zadoff-Chu matrix 

C1 = W(:,1);     %% Take the first column

orth_results = []; 
for j = 1 : size(W,2)
    results = (abs(sum(W(:,j).*C1)).^2);  %Project the first column to every column in matrix W
    orth_results = [orth_results results];
end 


Following the condition of orthogonality between columns, orth_results should be equals to 0 except in the fist value, but what I get is orth_results = [8.0000         0    8.0000      0];

What's the problem of that? Is there any issue in the matrix itself ?
For complex matrices, the concept of orthogonality is replaced by unitarity. If $^H$ denotes the  conjugate transpose of a matrix, we should check that $W$ is a unitary matrix (up to a scaling factor) via equations:
$$W^*  W = WW^*  =\lambda I$$
with $\lambda\neq 0$, which you can test positively with:
W'*W
W*W'

What is missing in your loop, and especially in the complex scalar product, is the conjugation on one of the complex vectors. This is explained for instance in Wikipedia/Dot product for Complex vectors.
You could for instance replace C1 by conj(C1).