Kuroiwa and Nakagawa analyse the two penalty forms already in use and separate them: one "works properly in that eigenstates obtained by the VQE with the penalty term reside in the desired symmetry sector", and they "further give a convenient formula to determine the magnitude of the penalty term, which may lead to the faster convergence of the VQE". The other is not a weaker option but a wrong one — they "prove that the other type of penalty terms does not work for obtaining the target state with the desired symmetry in a rigorous sense and even gives completely wrong results in some cases". Two conditions follow. The target symmetry sector has to be known in advance, since it is what the penalty is written against; and this route reaches a symmetry-resolved spectrum rather than the k-th state as such, so it answers "the lowest state with these quantum numbers" and not "the third state overall". Their own framing is that the work "lay[s] the theoretical foundation for the use of the VQE with the penalty terms", validated by numerical simulation.

A Hermitian H reachable as a sum of terms that can be measured separately; a statement of which state is wanted — an index k, a symmetry sector, or a target energy to sit nearest; for most of the methods here, the ground state already solved, and for the deflation route every lower state as well; a target additive error and a shot budget.

1. Hamiltonian whose excited state is wanted → Parameterised circuit familyChoose a parameterised trial state

   The constraint is not built into the family here. This route prepares an ordinary $|\psi(\boldsymbol{\theta})\rangle = V(\boldsymbol{\theta})|\psi_0\rangle$ with $V(\boldsymbol{\theta}) = V_l(\theta_l)\cdots V_2(\theta_2)V_1(\theta_1)$ and pays for the symmetry in the objective instead, because a symmetry-adapted ansatz "strongly depends on the heuristic knowledge of the symmetry and is not always readily available for generic symmetries". The guarantees are then stated over the whole space — not as a consequence of that choice, but because the purpose is to explore the performance of the cost functions themselves — assumption: the ansatz can represent an arbitrary $n$-qubit state, so $|\psi(\boldsymbol{\theta})\rangle = \sum_{i=0}^{2^n-1} a_i|i\rangle$ with $\sum_i |a_i|^2 = 1$ which reduces the parameters to populations $|a_i|^2$ on the joint eigenbasis of $\hat{H}$ and $\hat{C}$. The simulations use finite circuits: a hardware-efficient ansatz of depth $D = 4$ and $D = 12$ for $\mathrm{H}_2$ from $R_Y = e^{i\theta Y/2}$, $R_Z = e^{i\theta Z/2}$, and a real-valued symmetry-preserving depth-12 ansatz for $\mathrm{H}_4$ that conserves the electron number of $|\psi_0\rangle = |00001111\rangle$, so no $\hat{N}_e$ penalty is set.

   assumption

   Open this step

2. Parameterised circuit family → State you can prepareMinimise the objective over the parameters

   The objective here is not the energy. This route minimises $F^{(1)}(\boldsymbol{\theta}) = L(\boldsymbol{\theta}) + \sum_l \mu_{C^{(l)}}\langle\psi(\boldsymbol{\theta})|(\hat{C}^{(l)}-c^{(l)})^2|\psi(\boldsymbol{\theta})\rangle$ — one penalty per conserved quantity commuting with $\hat{H}$ — with $L$ either $L_{\mathrm{VQE}}$ or $L_{\mathrm{VQD}}^{(k)}$ assumption: the analysis is stated for the lowest state of a sector; VQD enters as the ground-state problem for $\hat{H}' = \hat{H} + \sum_{i=0}^{k-1}\beta_i|\psi_i\rangle\langle\psi_i|$, so the levels below are $\hat{H}'$'s. For a single $\hat{C}$ the populations give $F^{(1)} = \sum_{i=0}^{2^n-1}|a_i|^2(E_i + \mu_C(C_i-c)^2)$, the spectrum with every out-of-sector level lifted. Ordering $E_0 \leq \cdots \leq E_{2^n-1}$ and taking $i_0$ as the target index, so that $C_i \neq c$ for all $i < i_0$, any $\mu_C > 0$ with $\mu_C \geq \max_{i<i_0}(E_{i_0}-E_i)/(C_i-c)^2$ returns $E_{i_0}$; relaxed, $\mu_C^{(\mathrm{simple})} = (E_{i_0}-E_0)/C_{\min}^2$ over $C_{\min}$, the smallest gap between distinct eigenvalues of $\hat{C}$ ($3/4$ for $\hat{S}^2$, $1/2$ for $\hat{S}_z$). Since $E_{i_0}$ is the unknown, $E_{i_0}-E_0$ is either estimated by a cheaper classical method ($\mu_C^{(\mathrm{ce})}$) or bounded, $\mu_C^{(\mathrm{rough})} = 2\sum_j|c_j|/C_{\min}^2$ from $E_{i_0}-E_0 \leq 2\|\hat{H}\| \leq 2\sum_j|c_j|$: overestimating still works and only slows convergence, underestimating can move the minimum off the target. The rival $F^{(2)} = L + \mu_C(\langle\hat{C}\rangle-c)^2$ is quadratic in the populations where the energy is linear, so it penalises the deviation less; on the tangent edge $E-E_{i_0} = \alpha(C-c)$ of the convex envelope its global minimum sits at energy $E_{i_0}-\alpha^2/(2\mu_C)$ approximation: the target is recovered exactly only as $\mu_C\to\infty$, and only in case (A) with $i_0 \neq 0$, where $(c,E_{i_0})$ is a boundary point of that envelope; in case (B) no $\mu_C$ suffices.

   approximationassumption

   Open this step

3. State you can prepare → Number with an error barEstimate an observable

   Every evaluation of the objective is bought in Pauli measurements, and Eq. (5) makes this route's readout heavier than plain VQE's: beside $\hat{H} = \sum_j c_j P_j$ it must also estimate $\langle(\hat{C}-c)^2\rangle$, an operator that "consists of more Pauli operators than the original $\hat{C}$ does" — given as the main reason $F^{(2)}$, needing only $\langle\hat{C}\rangle$, is cheaper per evaluation. Cost is charged as (cost-function evaluations) $\times$ (Pauli operators per evaluation): on 8-qubit $\mathrm{H}_4$'s triplet, with BFGS and averaged over ten initial parameter sets, $F^{(1)}$ takes $1.28 \times 10^7$ measurements at $\mu_C = 0.01$ and $1.18 \times 10^8$ at $\mu_C = 100$, where the residual also degrades from $0.000000$ to $0.015919$ Ha. The same estimate doubles as a diagnostic: $\langle\psi(\boldsymbol{\theta})|(\hat{C}-c)^2|\psi(\boldsymbol{\theta})\rangle$ is tracked during the optimisation, and a value that "deviates from 0 drastically" after convergence marks an undershot $\mu_C$ and a missed sector. assumption: the simulations use exact expectation values and model neither circuit noise nor finite sampling Noise is treated separately: under a global depolarizing channel $\Delta_p(\rho) = (1-p)\rho + pI/2^n$, $\tilde{F}^{(1)} = (1-p)F^{(1)} + p\,\mathrm{tr}(\hat{H}+\mu_C(\hat{C}-c)^2)$ — a positive rescaling plus a $\boldsymbol{\theta}$-independent constant — so its minimiser is unchanged, while $F^{(2)}$'s shifts at $O(p)$.

   assumption

   Open this step

A scalar estimate — of that eigenvalue, or of the gap between it and the ground state — together with the run budget it consumed. Some methods return a whole low-energy subspace at once and others return one state per run, and that difference is a cost, not a convenience.

given  Hamiltonian H on n qubits and a conserved-quantity observable C with
       [H, C] = 0 (particle number N, S^2, S_z, ...)                     (Sec. II C)
       target symmetry sector: an eigenvalue c of C, fixed before the run --
           c determines the sector, it is an input here, not an output   (Sec. II C)
       ordinary ansatz  |psi(theta)> = V(theta)|ref>,
           V(theta) = V_l(theta_l) ... V_1(theta_1)                      (Eq. 1-2)
           # |ref> is the paper's |psi_0> of Eq. 1, "some reference state" -- one
           #   FIXED state (|00001111> in the H4 runs, Sec. IV) reused unchanged by
           #   every round. Renamed here because Sec. II B writes |psi_0> for a
           #   different object: the converged ground state found earlier.
           # no symmetry structure is built into V(theta); the guarantees below
           #   are then derived assuming this ansatz can reach an ARBITRARY
           #   n-qubit state, |psi(theta)> = sum_i a_i |i>, because the paper's
           #   purpose is to study the two cost functions themselves, not a
           #   specific circuit's reach                                  (Sec. III, Eq. 8)
       for an excited state: eigenstates |psi_0>, ..., |psi_{k-1}> already
           obtained by earlier rounds of this same procedure             (Sec. II B)

requires  the target sector known in advance as the scalar c -- this route
    resolves a symmetry sector, so it answers "the lowest state with C = c",
    not "the k-th state overall" without one; denote that state's energy
    E_{i0}, meaning C_i != c for every i < i0 under E_0 <= ... <= E_{2^n-1}  (Sec. III)

# --- cost function: penalize deviation from the sector, don't build it in ---------
if excited state k:
    L(theta) = LVQD(theta) = <psi(theta)|H|psi(theta)>
                   + sum_{i=0}^{k-1} beta_i |<psi_i|psi(theta)>|^2        (Eq. 4)
        # equivalently plain VQE against H' = H + sum_i beta_i |psi_i><psi_i|,
        #   beta_i > 0 chosen only to enforce orthogonality against the
        #   already-found lower states                                   (Sec. II B)
else:
    L(theta) = LVQE(theta) = <psi(theta)|H|psi(theta)>                   (Eq. 3)

F(theta) = L(theta) + mu_C * <psi(theta)|(C - c)^2|psi(theta)>  = F^(1)   (Eq. 5)
    # the rival form F^(2) = L + mu_C*(<psi|C|psi> - c)^2 is NOT used here:
    #   proved to miss the target by O(1/mu_C) at every finite mu_C even in
    #   the best case, and when (c, E_{i0}) sits interior to the (C,E) convex
    #   envelope, no mu_C -- however large -- recovers it at all   (Sec. III B, Eq. 18-23)

# --- set the penalty weight --------------------------------------------------------
# proved sufficient, but needs the full spectrum so is not directly usable:
#   mu_C >= max_{i<i0} (E_{i0} - E_i) / (C_i - c)^2                       (Eq. 11)
C_min = smallest gap between distinct eigenvalues of C
    # 1 for the particle-number operator N; 3/4 for S^2; 1/2 for S_z      (Sec. III A)
mu_C_simple = (E_{i0} - E_0) / C_min^2
    # a looser but simpler bound than Eq. 11, obtained from E_i >= E_0
    #   for all i < i0 -- still proved sufficient                        (Eq. 12-13)
# E_{i0} - E_0 is itself unknown before the run; two ways to fill it in:
mu_C_ce    = (E_{i0}_classical - E_0_classical) / C_min^2   # e.g. Hartree-Fock
    #   or Moeller-Plesset estimates of the two energies                 (Eq. 14)
mu_C_rough = (2 / C_min^2) * sum_j |c_j|     # H = sum_j c_j P_j (Pauli decomposition)
    #   from the rigorous bound E_{i0} - E_0 <= 2||H|| <= 2 sum_j |c_j|  (Eq. 15)
mu_C := mu_C_rough      # pick exactly ONE of the three; F^(1) takes a single scalar
    # overestimating E_{i0} - E_0 still finds the target, only slower;
    #   UNDERestimating can move the minimum off the target -- so
    #   mu_C_rough (never an underestimate) is the safe default          (Sec. III A)
    # but the paper also says mu_C_rough "may be too large for the fast
    #   convergence", and that a large mu_C slows the optimizer, so
    #   mu_C := mu_C_ce is the faster choice where the classical estimate
    #   can be trusted not to underestimate the gap             (Sec. III A, Eq. 15)
    # several conserved quantities at once: one term and one mu_C^(l) per
    #   C^(l), each set by the same formula with its own C_min^(l)       (Eq. 16-17)

# --- optimize -----------------------------------------------------------------------
run a classical optimizer on F^(1)(theta) with that mu_C, theta -> theta*
    # the paper's own numerics use BFGS on noiseless simulation with exact
    #   expectation values -- no shot noise, no device noise              (Sec. IV)
while optimizing, track <psi(theta)|(C - c)^2|psi(theta)>
    # if this has not fallen near 0 by convergence, mu_C was too small and
    #   the run missed the sector -- raise mu_C and rerun                 (Sec. III A)

return  |psi(theta*)>, F^(1)(theta*) as the estimate of the sector's E_{i0}
# under a global depolarizing channel this minimizer is UNCHANGED -- the noisy
#   F^(1) equals (1-p) times F^(1) plus a theta-independent constant, so the
#   same theta* minimizes it; F^(2)'s minimizer instead shifts at O(p)  (Sec. III C, Eq. 26)

Kuroiwa and Nakagawa price penalty-constrained VQE in Pauli measurements — (cost-function evaluations) $\times$ (Pauli operators per evaluation) — and bound neither factor: the working penalty $F^{(1)}$ measures $(\hat{C}-c)^2$ once per conserved quantity, more Pauli operators than $\hat{C}$ alone but never counted, plus $k$ overlap terms for the $k$-th state by VQD. Their figures are measured, not proved — 8-qubit $\mathrm{H}_4$ triplet, depth-12 symmetry-preserving ansatz, BFGS, noiseless, averaged over ten initial parameter sets. $F^{(1)}$ takes $1.28 \times 10^7$ Pauli measurements at $\mu_C = 0.01$, rising to $1.18 \times 10^8$ at $\mu_C = 100$ where the residual degrades from $0.000000$ to $0.015919$ Ha. That cheapest accurate point still exceeds the $9.78 \times 10^6$ the unguaranteed $F^{(2)}$ needs at $\mu_C = 100$, though the ordering reverses in their $k = 4$ VQD run. On $\mathrm{H}_4$ with BFGS the classically estimated $\mu_C$ needs $17{,}604$ evaluations against $147{,}063$ for the loose bound, which converged less accurately. Every guarantee assumes an ansatz spanning arbitrary $n$-qubit states.

• Penalty-based excited-state VQE · Qiskit

  From the repository — run, not written up from a paper · unsupported

  About

  Orthogonality or symmetry penalties augment the energy objective to exclude previously identified sectors.

  Methods

  None found yet.

  Data

  None found yet.

  Code

  Qiskit

  Results

  Literature-backed method record; algorithmic scope and evidence boundary reviewed, with no benchmark run claimed.

• Penalty-based excited-state VQE

  Orthogonality or symmetry penalties augment the energy objective to exclude previously identified sectors.