# $S^2$ expectation value of a circuit seems wrong

I need to rewrite a following circuit from Cirq to Qiskit with $$2\varphi = \pi$$.

It's a circuit representing an initial state for VQE calculation. Considering, that qiskit_nature considers a different ordering of spin-orbitals (first all spin-up, then all spin-down), I've rewritten the circuit like this (switching 1st and 2nd qubit, indexing from 0):

In both cases the statevector equals to

array([ 0.        , -0.        ,  0.        , -0.        ,  0.        ,
-0.        ,  0.70710678, -0.        ,  0.        , -0.70710678,
0.        , -0.        ,  0.        , -0.        ,  0.        ,
-0.        ])


, i.e.

$$\left|\psi\right> = \frac{1}{\sqrt{2}}\left| 0110 \right> - \frac{1}{\sqrt{2}}\left| 1001 \right>.$$

In Cirq, the total angular momentum $$\left<\psi \middle| S^2 \middle| \psi \right> = 0$$, i.e. the state should be singlet. But when I try the same in qiskit with the help of second_q_ops(), I'm getting 1.999999....

Why does my result differ?

### My Code

#!/usr/bin/env python3
import numpy as np
import qiskit_nature
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.providers.aer import Aer
from qiskit.quantum_info import Statevector
from qiskit_nature.converters.second_quantization import QubitConverter
from qiskit_nature.drivers import Molecule
from qiskit_nature.drivers.second_quantization import ElectronicStructureMoleculeDriver, ElectronicStructureDriverType
from qiskit_nature.mappers.second_quantization import JordanWignerMapper
from qiskit_nature.problems.second_quantization import ElectronicStructureProblem
from qiskit_nature.transformers.second_quantization.electronic import ActiveSpaceTransformer

if __name__ == '__main__':

qiskit_nature.settings.dict_aux_operators = True
backend = Aer.get_backend('statevector_simulator')

geometry = [('N', [0.000000000000, 0.000000000000, 0.000000000000]),
('C', [0.000000000000, 0.000000000000, 1.498047000000]),
('H', [0.000000000000, -0.938765985000, 2.004775984000]),
('H', [0.000000000000, 0.938765985000, 2.004775984000]),
('H', [-0.744681452, -0.131307432, -0.634501434])]

# Obtaining ground state
driver = ElectronicStructureMoleculeDriver(Molecule(geometry),
basis="sto3g",
driver_type=ElectronicStructureDriverType.PSI4)

as_transformer = ActiveSpaceTransformer(num_electrons=2, num_molecular_orbitals=2)
es_problem = ElectronicStructureProblem(driver, transformers=[as_transformer])
converter = QubitConverter(JordanWignerMapper())

# Create circuit
qreg_q = QuantumRegister(4, 'q')
circuit = QuantumCircuit(qreg_q)

circuit.x(qreg_q[2])
circuit.ry(np.pi, qreg_q[0])
circuit.ch(qreg_q[0], qreg_q[3])
circuit.cx(qreg_q[3], qreg_q[2])
circuit.cx(qreg_q[3], qreg_q[0])
circuit.cx(qreg_q[0], qreg_q[1])
circuit.z(qreg_q[3])
circuit.x(qreg_q[0])
print(f'Circuit:\n{circuit}')

statevec = Statevector(circuit).data
print(f'Statevector:\n{statevec}')

# Obtain S^2 operator
ss_op = es_problem.second_q_ops()['AngularMomentum']
print(f'S^2 operator:\n{ss_op}')

print(f'<psi|S^2|psi> = {statevec @ converter.convert(ss_op).to_matrix() @ statevec}')


### The Output

Circuit:
┌───────┐          ┌───┐     ┌───┐
q_0: ┤ Ry(π) ├──■───────┤ X ├──■──┤ X ├
└───────┘  │       └─┬─┘┌─┴─┐└───┘
q_1: ───────────┼─────────┼──┤ X ├─────
┌───┐    │  ┌───┐  │  └───┘
q_2: ──┤ X ├────┼──┤ X ├──┼────────────
└───┘  ┌─┴─┐└─┬─┘  │  ┌───┐
q_3: ─────────┤ H ├──■────■──┤ Z ├─────
└───┘          └───┘
Statevector:
[ 0.00000000e+00+0.j  0.00000000e+00+0.j  0.00000000e+00+0.j
0.00000000e+00+0.j  0.00000000e+00+0.j  6.12323400e-17+0.j
7.07106781e-01+0.j  0.00000000e+00+0.j  0.00000000e+00+0.j
-7.07106781e-01+0.j  0.00000000e+00+0.j  0.00000000e+00+0.j
0.00000000e+00+0.j  0.00000000e+00+0.j  0.00000000e+00+0.j
0.00000000e+00+0.j]
S^2 operator:
Fermionic Operator
register length=4, number terms=12
(0.75+0j) * ( +_0 -_0 )
+ (0.75+0j) * ( +_1 -_1 )
+ (0.75+0j) * ( +_2 -_2 )
+ (0.75+0j) * ( +_3 -_3 )
+ (0.5+0j) * ( +_0 -_0 +_1 -_1 )
+ (-1.5+0j) * ( +_0 -_0 +_2 -_2 )
+ (-0.5+0j) * ( +_0 -_0 + ...
<psi|S^2|psi> = (1.9999999999999996+0j)