Resource Estimation for Spectroscopy Applications

Simulating Spectroscopy on a Quantum Computer

First challenge is translating a physical spectroscopy experiment into a quantum circuit.

Regardless of the specific technique, whether it is X-ray Absorption Spectroscopy(XAS) [1], Vibrational Spectroscopy [3], or Electron Energy Loss Spectroscopy [2], the core goal is the same: calculating the time-domain correlation function, \(\tilde{G}(t)\).

To measure this property on a quantum computer, we rely on a standard algorithmic template called the Hadamard Test.

The dominant cost for this circuit comes from the controlled time evolution block, the efficiency of which thus dictates the resource requirements for our algorithms. This cost is dictated by the intersection of the spectroscopic domain and our implementation strategy.

While the physical problem of interest dictates the type of Hamiltonian (e.g., electronic, vibrational etc.), we have significant freedom to optimize the implementation. We can select specific Hamiltonian representations and pair them with the time evolution algorithm of choice, i.e. Trotterization or Qubitization.

We begin with the example of X-ray Absorption Spectroscopy (XAS) to demonstrate how to use the PennyLane resource estimator to generate logical resource counts.

X-Ray Absorption Spectroscopy

XAS, is a critical spectroscopic method used to study the electronic and local structural environment of specific elements within a material, by probing core-level electron transitions. For this simulation, we follow the algorithm established in Fomichev et al. (2025) [1], which utilizes the Compressed Double Factorization (CDF) representation of the electronic Hamiltonian combined with Trotterization.

To benchmark this approach, we focus on Lithium Manganese (LiMn) oxide clusters, which are widely studied as critical cathode materials for next-generation batteries.

The first step to determining the resources for this simulation is to define the Hamiltonian. Here, we must note that the resource requirements for the simulation depend on the structural parameters of the Hamiltonian, specifically the number of orbitals and the number of fragments, rather than the exact integral values. We leverage this by utilizing the specialized compact Hamiltonian representation feature offered by PennyLane, skipping the expensive Hamiltonian construction while retaining the exact cost topology required for analysis.

import pennylane.estimator as qre

active_spaces = [11, 14, 16, 18]
limno_ham = [qre.CDFHamiltonian(num_orbitals=i, num_fragments=i) for i in active_spaces]

import numpy as np

eta = 0.05  # Lorentzian width (experimental resolution) in Hartree
jmax = 100  # Number of time points (determines resolution limit)
Hnorm = 2.0  # Maximum final state eigenvalue used to determine tau.

tau = np.pi / (2 * Hnorm)  # Sampling interval
trotter_error = 1  # Hartree
delta = np.sqrt(eta / trotter_error)  # Trotter step size
num_trotter_steps = int(np.ceil(2 * jmax + 1) * tau / delta)  # Number of Trotter steps

def xas_circuit(hamiltonian, num_trotter_steps, measure_imaginary=False, num_slaters=1e4):

    # State preparation
    num_qubits = int(np.ceil(np.log2(num_slaters)))
    qre.QROMStatePreparation(
        num_state_qubits=num_qubits,
        positive_and_real=False,
        select_swap_depths=1,
    )
    qre.QROM(num_bitstrings=2**num_qubits, size_bitstring=num_qubits, select_swap_depth=1)

    # Hadamard and S gates
    qre.Hadamard()

    if measure_imaginary:
        qre.Adjoint(qre.S())

    # Controlled time evolution
    qre.Controlled(
        qre.TrotterCDF(hamiltonian, num_trotter_steps, order=2), num_ctrl_wires=1, num_zero_ctrl=0
    )

    qre.Hadamard()

    # Uncompute state preparation
    qre.Adjoint(
        qre.QROMStatePreparation(
            num_state_qubits=num_qubits,
            positive_and_real=False,
            select_swap_depths=1,
        )
    )
    qre.Adjoint(
        qre.QROM(num_bitstrings=2**num_qubits, size_bitstring=num_qubits, select_swap_depth=1)
    )

Resource Estimation

With the circuit fully defined, we turn to cost estimation. Before generating counts, however, we must select an implementation strategy for the rotation gates, as this choice heavily influences the final resource overhead.

PennyLane’s default compiler synthesizes rotation gates using **Repeat-Until-Success circuits [#Alex2014], which decompose rotations into sequences of probabilistic T-gates. While effective for general circuits, the algorithm we are implementing specifically calls for the Phase Gradient Trick proposed by Craig Gidney [6],

The phase gradient trick is algorithmically superior for this application because it allows us to implement rotations with deterministic cost using arithmetic, rather than relying on probabilistic synthesis sequences. This results in better scaling as the precision requirements increase.

To adopt this more efficient strategy, we configure the estimator to use the adder-based synthesis by leveraging ResourceConfig.

def single_qubit_rotation(precision=None):
    """Gidney-Adder based decomposition for single qubit rotations"""
    num_bits = int(np.ceil(np.log2(1 / precision)))
    return [qre.GateCount(qre.resource_rep(qre.SemiAdder, {"max_register_size": num_bits + 1}))]

We can now set up the resource estimation with this custom decomposition using the set_decomp function and also set the targeted precision:

cfg = qre.ResourceConfig()
cfg.set_decomp(qre.RX, single_qubit_rotation)
cfg.set_decomp(qre.RY, single_qubit_rotation)
cfg.set_decomp(qre.RZ, single_qubit_rotation)
cfg.set_single_qubit_rot_precision(1e-3)

With the configuration set, we run the resource estimation by sweeping over the different active-spaces of our Hamiltonian.

xas_resources = []
toffolis = []
qubits = []
for ham in limno_ham:
    resource_counts = qre.estimate(xas_circuit, config=cfg)(
        hamiltonian=ham, num_trotter_steps=num_trotter_steps, measure_imaginary=False
    )
    xas_resources.append(resource_counts)
    toffolis.append(resource_counts.gate_counts["Toffoli"])
    qubits.append(resource_counts.total_wires)

Let’s visualize how these initial estimates scale with the system size.

These results highlight that while qubit requirements (~90-100) are feasible for early fault-tolerant devices, the gate complexity approaches \(10^9\) Toffolis.

Crucially, these estimates represent a single shot; the total cost for the full algorithm will scale linearly with the number of samples required. This gate overhead is therefore the primary bottleneck we must address. Let’s see how we can optimize these gate counts further.

Optimizing the Estimates

In order to optimize our algorithm, instead of accepting the standard implementations, we can inject specific high-performance subroutines directly into the resource estimator.

Optimization 1: Efficient Basis Rotations¶

First, we target the basis rotation, and replace the generic standard decomposition with the specialized, lower-cost circuit described in Kivlichan et al. (2018).

def basis_rotation_cost(dim):
    """Custom basis rotation decomposition from Kivlichan et al. (2018)"""

    counts = dim * (dim - 1) / 2

    ops_data = [
        (qre.RZ, dim, None),
        (qre.RX, counts, None),
        (qre.RY, counts, None),
        (qre.Hadamard, 4 * counts, None),
        (qre.S, counts, None),
        (qre.Adjoint, counts, {"base_cmpr_op": qre.resource_rep(qre.S)}),
        (qre.CNOT, 2 * counts + dim, None),
    ]

    return [
        qre.GateCount(qre.resource_rep(op, params or {}), count) for op, count, params in ops_data
    ]

cfg.set_decomp(qre.BasisRotation, basis_rotation_cost)

toffolis_opt1 = []
for ham in limno_ham:
    res = qre.estimate(xas_circuit, config=cfg)(hamiltonian=ham, num_trotter_steps=num_trotter_steps, measure_imaginary=False)
    toffolis_opt1.append(res.gate_counts["Toffoli"])

Optimization 2: Double Phase Trick¶

Second, we implement the double phase trick for the Controlled-RZ (CRZ) gates. As described in Section III A of Fomichev et al. (2025) [1], this optimization reduces the cost of the controlled rotations inside the Trotter steps by combining phase shifts.

def custom_CRZ_decomposition(precision):
    """Decomposition of CRZ gate using double phase trick"""
    rz = qre.resource_rep(qre.RZ, {"precision": precision})  # resource representation of RZ
    cnot = qre.resource_rep(qre.CNOT)

    return [qre.GateCount(cnot, 2), qre.GateCount(rz, 1)]

cfg.set_decomp(qre.CRZ, custom_CRZ_decomposition)

toffolis_final = []
for ham in limno_ham:
    res = qre.estimate(xas_circuit, config=cfg)(hamiltonian=ham, num_trotter_steps=num_trotter_steps, measure_imaginary=False)
    toffolis_final.append(res.gate_counts["Toffoli"])

# Let's visualize the cumulative impact of our optimizations:
import matplotlib.pyplot as plt

plt.plot(active_spaces, toffolis, "o-", label="Baseline", color="gray", linewidth=2.5)
plt.plot(active_spaces, toffolis_opt1, "^-", label="Optimization 1 (Basis Rot.)", color="goldenrod", linewidth=2.5)
plt.plot(active_spaces, toffolis_final, "*-", label="Fully Optimized", color="fuchsia", linewidth=2.5, markersize=8)

plt.xlabel("Number of Orbitals")
plt.ylabel("Toffoli Gate Count")
plt.title("XAS Resource Estimation Comparison")
plt.legend()
plt.show()

The plot illustrates the effectiveness of our optimization strategy. The specialized basis rotation (gold) delivers the most significant reduction, visibly correcting the slope of the cost scaling compared to the baseline. The double phase trick then provides a final constant-factor improvement (pink), bringing the total gate count even further down. This stepwise reduction validates the importance of matching the gate synthesis strategy to specific algorithmic requirements, effectively “engineering” a lower simulation cost by targeting dominant bottlenecks.

Photodynamic Therapy Applications

Having validated our resource estimation workflow against the XAS benchmarks, let’s see how we can apply these tools to a different application, i.e. Photodynamic Therapy (PDT), where we need to estimate the cumulative absorption rates for a transition-metal photosensitizer. [7]

This application requires a fundamental shift in algorithmic strategy. While XAS simulates the system’s time evolution to observe dynamic changes (Trotterization), PDT employs a spectral filtering approach. Here, rather than resolving individual eigenstates, we use Generalized Quantum Signal Processing (GQSP) to isolate and measure the total signal within a specific therapeutic energy window (typically 700–850 nm).

To achieve this efficiently, the algorithm utilizes a walk operator constructed from the tensor hypercontracted Hamiltonian, which allows for a highly compact block encoding. Similar to the XAS example, we can use the compact Hamiltonian representation to skip the expensive Hamiltonian construction. Let’s verify our model using the 11-orbital BODIPY system from the reference.

bodipy_ham = qre.THCHamiltonian(num_orbitals=11, tensor_rank=22, one_norm=6.48)

We now construct the walk Operator from the Hamiltonian using the ~.pennylane.estimator.templates.QubitizeTHC template. For comprehensive details on how to construct and configure this operator, we recommend the Qubit and gate trade-offs in Qubitized Quantum Phase Estimation demo. Let’s define the precision parameters based on the error budget from reference and construct the walk operator accordingly:

error = 0.0016  # Error budget from Zhou et al. (2025)
n_coeff = int(np.ceil(2.5 + np.log2(10 * bodipy_ham.one_norm / error)))  # Coeff precision
n_angle = int(
    np.ceil(5.652 + np.log2(10 * bodipy_ham.one_norm * 2 * bodipy_ham.num_orbitals / error))
)  # Rotation angle precision

prep_op = qre.PrepTHC(bodipy_ham, coeff_precision=n_coeff, select_swap_depth=4)
select_op = qre.SelectTHC(bodipy_ham, rotation_precision=n_angle, num_batches=5)
walk_op = qre.QubitizeTHC(bodipy_ham, prep_op=prep_op, select_op=select_op)

Next, we need to set up the GQSP parameters to construct the spectral filter. The filter is defined by a polynomial whose degree determines the sharpness of the filter. Following Zhou et al. (2025) [7], we set the polynomial degree using Figure 7 from the paper, which relates the degree to the desired spectral resolution.

def polynomial_degree(one_norm):
    """Calculate polynomial degree parameters from Zhou et al. (2025)"""
    e_hi = 0.0701  # Ha
    e_lo = 0.0507  # Ha
    degree_hi = int(np.ceil(4.7571 * (one_norm + e_hi) / 0.01 + 321.2051))
    degree_low = int(np.ceil(4.7571 * (one_norm + e_lo) / 0.01 + 321.2051))

    return degree_hi, degree_low

We can now set up the resource estimation for the PDT algorithm by defining the circuit as shown in Figure 2.

def pdt_circuit(walk_op, poly_degree_hi, poly_degree_low, num_slaters=1e4):
    num_qubits = int(np.ceil(np.log2(num_slaters)))
    qre.QROMStatePreparation(
        num_state_qubits=num_qubits,
        positive_and_real=False,
    )
    qre.QROM(num_bitstrings=2**num_qubits, size_bitstring=num_qubits, select_swap_depth=1)

    # Hadamard
    qre.Hadamard()

    # GQSP Spectral Filter
    qre.GQSP(walk_op, d_plus=poly_degree_hi)
    qre.GQSP(walk_op, d_plus=poly_degree_low)

    # Multi-Controlled-X
    qre.MultiControlledX(
        num_ctrl_wires=2,
        num_zero_ctrl=1,
    )

    # Hadamard
    qre.Hadamard()

    # Uncompute state preparation
    qre.Adjoint(
        qre.QROMStatePreparation(
            num_state_qubits=num_qubits,
            positive_and_real=False,
        )
    )
    qre.Adjoint(
        qre.QROM(num_bitstrings=2**num_qubits, size_bitstring=num_qubits, select_swap_depth=1)
    )

Now, that we have all the components set up, we can run the resource estimation for the PDT algorithm to calculate cumulative absorption:

qubits_pdt = []
toffolis_pdt = []

poly_deg_hi, poly_deg_low = polynomial_degree(bodipy_ham.one_norm)  # Calculate polynomial degree
resource_counts = qre.estimate(pdt_circuit, config=cfg)(walk_op, poly_deg_hi, poly_deg_low)
print(resource_counts)

--- Resources: ---
 Total wires: 174
   algorithmic wires: 78
   allocated wires: 96
     zero state: 96
     any state: 0
 Total gates : 2.298E+8
   'Toffoli': 1.965E+7,
   'CNOT': 1.491E+8,
   'X': 7.846E+6,
   'Z': 1.786E+5,
   'S': 3.022E+5,
   'Hadamard': 5.280E+7

The estimated resources of 174 qubits and \(1.96 \times 10^7\) Toffoli gates align with the reference values of 177 qubits and \(2.72 \times 10^7\) Toffoli gates. The small differences can be attributed to different tunable parameters being used. We encourage users to explore further by testing other systems from the reference or analyzing how the resources scale with different error budgets, using the parameter tuning techniques detailed in our qubitization demo.

Conclusion

In this demo, we successfully performed end-to-end resource estimation for two distinct spectroscopic paradigms: the time-domain simulation of X-ray Absorption (XAS) and the spectral filtering approach for Photodynamic Therapy (PDT).

Beyond providing specific resource counts, this demo highlights the flexibility of the resource estimator as a programmable design tool. We showed how to move beyond “black box” standard implementations by actively prototyping advanced algorithmic choices—such as the phase gradient trick and specialized basis rotations. This flexibility empowers researchers to “engineer” their quantum algorithms, identifying cost bottlenecks and seamlessly swapping out subroutines to optimize performance for future fault-tolerant hardware, regardless of the underlying spectroscopic technique.

References