first update:

am now looking into the (nicely written it seems) lecture notes linked below. I hope I shall find an answer there provided I make it to the end :-)

Actually, eq.(17.51) in this document provides a first clue in that already for pure orbital angular momentum states there is a (-1)^(ML) dependence in the relation between time reversed states. If that carries through my problem is probably solved and the answer would be that whether or not a minus sign appears in the relation between Kramers conjugate states depends on the value of MJ (e.g. through J-MJ being even or odd).

However, for half integral angular momenta (odd number of electrons) things are a bit more intricate and I'll have to give it some more time.

http://bohr.physics.berkeley.edu/classes/221/9697/timerev.pdf

I am surprised that the phase of the wave function affects the energy of the state. The absolute phase should be observable only if the wave function hybridizes/overlaps with something else. Is this "something else" you representation of the crystal fields?

No, energies do not depend on the phases. But what operators look like does depend on it. And that in turn may affect the signs of coefficients describing the eigenvectors in terms of the base states.

What happened in our case was that we obtained twofold degenerate eigenvectors which to me did not look like they were Kramers conjugate. So I suspected an error in our implementation of the Hamiltonian (hence my desire to actually check Kramers conjugatedness of these states).

Many measurements don't even depend on these signs. We would get the same XA spectra and also magnetization would be independent (all these measurements being dipolar in character). The spatial distribution of (4f, in our case) charge density, however, may look very different. In case of Cerium, in a particular case it was shown that x-ray scattering with quadrupolar sensitivity (NIXS) does have the sensitivity to reveal the phases in the linear combinations of bases states forming the eigenstates.

So from a practical point of view, one might say "don't care, shut up and calculate". But then again, if the solutions don't look proper, I am left with the uncertainty of not knowing what went wrong. And here, the answer seems to be: nothing went wrong but I was uninformed about how the (implicitely defined) phases of my base vectors relate to those pertaining to time reversal.

For those interested I enclose the link to the NIXS paper below. In that case, the rotational symmetry was fourfold, while ours is threefold. This does make a difference in what the operators (matrices) look like with repect to a sign reversal of MJ.

http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.109.046401

In the Willers paper you cite the wave function in Fig. 1 contradicts the text on page 1, right column (definition of Gamma^7_1 and _2,  correct as far as I can see). The difference between _1 and _2 is the relative sign of the +/- 5/2 and -/+ 3/2 contributions. You could also write alpha and sqrt(1-alpha^2) as sin(theta) and cos(theta). _1 and _2 would then be 90deg offset in theta.

I guess you have seen this, and you question is: Are the states as given with +/- 5/2 and -/+3/2 really the proper Kramers pair? Is that correct?

Note also that as long as time reversal is completely intact and the two wave functions are completely degenerate you could chose a different basis in the Gamma^7_1 space. These two functions would then not be related by time reversal. This freedom of choice evaporates as soon as there is a magnetic field that lifts the degeneracy between the two wave functions.

(Messiah, Vol. 2 chapter XV section 22 may also help)

"and you question is: Are the states as given with +/- 5/2 and -/+3/2 really the proper Kramers pairs?"

Somehow along these lines, yes. What (I think) I can tell by now is that suppose phases in our basis set are defined such that

|5/2>  and  |-5/2>  are a "proper" Kramers pair then for |m|=3/2

|3/2> and  -|-3/2>  are a proper Kramers pair.

This implies that if CF interaction creates a mixed state

|psi> = a* |5/2> + b*|-3/2>

then the Ktamers conjugate state will be

|psi'> = a*|-5/2> - b*|3/2>

This is exactly the kind of result we got (but for looking into the case of Erbium) and seemed strange and incorrect to me because I was unaware of the above relations. Now I have semi-understood that it actually must be like this.

As you say, whether or not this result is correct depends on the phase conventions in the base vectors.  Again, in coding the problem, the phase conventions are implicitly made by the rules we apply to construct the angular momentum operators.

For me, the question is practically resolved, since I got aware of the above relations, though I haven't explicitly made a derivation for arbitrary 4f configurations...

From what I understand of Messiah's book, J-states with angular and spin contributions are pretty much a mess with the time-reversed wave function being, in general, multiplied not just with +/-1 but a complex phase factor that depends on the details of what L and S states are mixed in. Kramer' s degeneracy seems to be directly related to the fact that for a spin T^2 = -1, whereas for orbital momenta T^2=1.

So the key seems to be to understand this factor. The CF eigenstates are obtained by diagonalizing the CF Hamiltonian, and the degeneracy is found even without taking time reversal into account, as nicely shown by Lea, Leask, and Wolf.