Flaky test_prepare_gaussian_state due to missing qubit_order in final_state_vector #1454

Issue Title:

[Bug] Flaky test_prepare_gaussian_state due to missing qubit_order in final_state_vector

Issue Description:

Bug Summary

test_prepare_gaussian_state in src/openfermion/circuits/primitives/state_preparation_test.py fails probabilistically under certain random seeds (e.g. --randomly-seed=1921347655).

The failure manifests as a SciPy sparse matrix multiplication dimension mismatch: ValueError: matmul: dimension mismatch with signature (n,k=16),(k=2,1?)->(n,1?)

Environment Details

• OS: macOS / Linux
• Python: 3.12+
• Packages: cirq, scipy, pytest-randomly

Steps to Reproduce

Run pytest using the specific random seed:

bash

pytest src/openfermion/circuits/primitives/state_preparation_test.py --randomly-seed=1921347655

Failure Log / Traceback

text

__________________ test_prepare_gaussian_state[4-True-None-2] __________________

n_qubits = 4, conserves_particle_number = True, occupied_orbitals = None

initial_state = 2, atol = 1e-05

# Get the state using a circuit simulation

circuit = cirq.Circuit(

    prepare_gaussian_state(qubits, quad_ham, occupied_orbitals, initial_state=initial_state)

)

if isinstance(initial_state, list):

    initial_state = sum(1 << (n_qubits - 1 - i) for i in initial_state)

state = circuit.final_state_vector(initial_state=initial_state)

# Check that the result is an eigenstate with the correct eigenvalue

numpy.testing.assert_allclose(quad_ham_sparse.dot(state), energy * state, atol=atol)

E ValueError: matmul: dimension mismatch with signature (n,k=16),(k=2,1?)->(n,1?)

/Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/site-packages/scipy/sparse/_base.py:654: ValueError

Root Cause Analysis

In state_preparation_test.py:

python

state = circuit.final_state_vector(initial_state=initial_state)

1. In Cirq, calling circuit.final_state_vector() without specifying qubit_order automatically infers the qubit order and active qubit set from circuit.all_qubits() (only the qubits that have operations applied to them).
2. For specific random Hamiltonians (e.g. under seed 1921347655), the state preparation circuit leaves 3 out of 4 qubits idle, applying operations to only 1 qubit.
3. Cirq therefore returns a 2-element state vector ($2^1 = 2$) instead of a 16-element state vector ($2^4 = 16$).
4. Multiplying quad_ham_sparse (dimension $16 \times 16$) by state (dimension $2$) raises ValueError: matmul: dimension mismatch.

I'd be happy to open a quick PR to fix this if no one else is currently working on it