# Bond order correlation function

I am trying to compute the bond order correlation function, $$g_6$$. It is defined based on the bond order parameter: $$\psi_6(x_i) = \frac{1}{N_i}\sum_{i=1}^{N_i}{\exp(i6\theta_i^j)}$$ where $$\theta_i^j$$ is defined as the bond angle between closest neighbors to the particle located at $$x_i$$. Then, we have: $$g_6(x=|x_i - x_j|) =\langle\psi_6(x_i) \psi_6^*(x_j)\rangle .$$ I am calculating this parameter for my simulation results; the problem is that the $$g_6$$ function should give real values, in all papers that I have seen, they plot it in a log-log diagram. But it is obvious to me that $$\psi_6(x_i) \psi_6^*(x_j)$$ does not yield a real positive number necessarily, since they are for two different particles and we are producting two different complex numbers. It is not like $$|a|=a.a^*$$. And I am getting complex values for my results. Could you help me understand where I am making a mistake? If you want to know more about this function you can take a look at this paper : http://dspace.nbuv.gov.ua/bitstream/handle/123456789/32106/04-Brodin.pdf

In the context in which this usually appears, the intention is to define this function to be symmetric to the exchange of indices $$i$$ and $$j$$. In other words, it doesn't matter which of the two terms comes first, and which second, in your expression. This has the same effect as taking the complex conjugate. So the $$ij\leftrightarrow ji$$ symmetrized version would be real.
The ensemble average would normally take care of this symmetry, and in any case there is often a sum over $$i$$ and $$j$$ pairs in which both the $$ij$$ and $$ji$$ terms appear explicitly. But if you are worried about the imaginary parts of the expression, you could simply symmetrize it (take the real part) before averaging in the simulation. This should make no difference, as we can be confident that the imaginary part will average to zero.