Question on using_auto_fermion

I know that using_auto_fermion is experimental (A naiove query on using_auto_fermion), but I still want to get more understanding of how it works so that I can do some tests on it or implement a similar version myself. After reading related codes( NDTensors/src/blockspare/linearalgebra.jl, NDTensors/src/blockspare/contract_generic.jl ,src/fermions.jl) and paper SciPost: SciPost Phys. 18, 012 (2025) - Fermionic tensor network methods, I have a basic question:

Is the contraction of the fermionic tensor in ITensor just the fermionic contraction, i.e, the \mathcal{C} in that paper, which satisfys \mathcal{C}(|j\rangle \langle i|)=(-1)^{|i||j|}\mathcal{C}(\langle i||j\rangle)=(-1)^{|i|}\delta_{ij} where |i|=0,1 for fermionic even/odd subspace.
Otherwise, does it automatically include the parity tensor \mathcal{P}=(-1)^{|i|} \delta_{ij} |i\rangle\langle j| when the direction of contraction is different from the direction of the legs to be contracted?

I have this problem because I find that the square of the norm of a fermionic tensor c is different from contracting c and c^\dagger, which suggests that the \mathcal{P} is not included, and this contraction is not an inner product: as mentioned in that paper, \langle A| B\rangle \neq \mathcal{C} (A^\dagger B). Instead, a parity tensor should be placed and \langle A| B\rangle = \mathcal{C} (A^\dagger \mathcal{P}(B) ).

using ITensors: op, siteind, norm, dag, enable_auto_fermion
enable_auto_fermion()
C = op(siteind("Fermion";conserve_qns=true), "C")
norm(C) # = 1.0
Array(dag(C)*C)# = -1.0
Array(C*dag(C))# = 1.0

This seems consistent with c=|0\rangle\langle 1| and \mathcal{C}(|1\rangle\langle 0||0\rangle\langle 1|) = \mathcal{C}(|1\rangle\langle 1|)=-1, \mathcal{C}(|0\rangle\langle 1||1\rangle\langle 0|)= \mathcal{C}(|0\rangle\langle 0|)=1 and \mathcal{P} is not included.

However, when I look at the compute_alpha function in src/fermions/fermions.jl, it looks like that in
α = alpha1 * alpha2 * alphaR * alpha_arrows the alpha_arrows already takes the \mathcal{P} into consideration. Besides, when I read the MPS inner product routine _log_or_not_dot in abstractmps.jl, I find that the inner product is calculated just by O = M1dag[1] * M2[1] , e.t.c.

It seems that there is a contradiction, or I misunderstood how the fermionic tensor works.

Thanks for any possible comments and suggestions!

I guess the problem comes from the fact that, in that paper, the tensors are assumed to have zero total flux, but in ITensor a tensor can have nonzero flux. I find that for MPS that has total parity even inner(M1,M2)=conj(inner(M2,M1)) as expected, while for MPS with total parity odd inner(M1,M2) = -conj(inner(M2,M1)).

For your question about what formalism we are using, I’d refer you to the following slides I made for a group meeting at Flatiron:
https://itensor.org/miles/AutoFermion_ITensor_TensorMtg.pdf
Please let me know if the slides don’t answer your questions about the formalism.

About your second post, yes MPS with overall odd parity can have behaviors like the one you mentioned (e.g. inner(M1,M2) = - conj(inner(M2,M1))). Actually it can be even worse than that: the result can depend on how the code is written, somewhat analogous to the code not being “gauge invariant”.
So I am considering actually forbidding odd-parity tensors, which would also forbid odd-parity MPS. While their behaviors are mathematically well-defined, they are just extremely difficult to use correctly especially within the ITensor system where index ordering is meant to be unknown to the user / programmer.