arXiv CSopen access

Stacking the Deck: Tunable Trainability in Stacked LCUs

Stacking the Deck: Tunable Trainability in Stacked LCUs Nikhil Khatri,1 Stefan Zohren,1 and Gabriel Matos2 1

Machine Learning Research Group, University of Oxford, United Kingdom 2 Quantinuum, London, United Kingdom (Dated: 27 July 2026)

arXiv:2607.24686v1 [quant-ph] 27 Jul 2026

Variational quantum circuits have been central to many proposed near-term applications of quantum computing, but a growing body of evidence suggests that trainability and quantum advantage are fundamentally at odds: ansätze expressive enough to resist efficient classical simulation tend to exhibit barren plateaus, while structures that provably rule out barren plateaus typically render them classically simulable. We propose a stacked linear combination of unitaries (S-LCU) as a variational ansatz which provides a tunable trade-off between barren plateaus and classical simulability. Using a diagrammatic analysis, we bound the loss-landscape variance of the Free Fermion S-LCU, whose elements are fermionic Gaussian unitaries. We prove a variance lower bound of Ω(1/(nk3l )), with a simulation cost of O(k2l n3 ) using the best known classical algorithm, compared to a quantum gate complexity of only O(lkn2 ). The number of layers l serves as a single dial that trades computational complexity against the rate of cost concentration. This offers practitioners a systematic method for constructing ansätze with a complexity-trainability trade-off that best suits their application and hardware.

I.

INTRODUCTION

Variational quantum algorithms optimise a parametrised quantum circuit against a cost function, and are an important tool in many proposed near-term applications of quantum computing, from quantum machine learning [1, 2] to ground state preparation [3–5] and combinatorial optimisation [6–8]. Their utility as quantum routines relies on a balance between two demands: a parameterised circuit must be efficiently trainable, yet be difficult to simulate classically. These goals are increasingly understood to be in tension. Models with sufficiently high expressivity tend to exhibit barren plateaus, where the variance of the gradient vanishes exponentially with system size, rendering training from a random parameter initialisation infeasible [9–11]. Conversely, the structure that provably rules out barren plateaus typically also enables efficient classical simulation [12]. One route to bridging this divide builds on the linear combination of unitaries (LCU) primitive [13]. A number of variational ansätze incorporate LCUs, optimising the underlying parameterised unitaries, the LCU coefficients, or both [14–18]. Recent work established that an LCU of barren plateau-free circuits is also free of barren plateaus, while increasing expressivity and potentially increasing classical simulation cost [19]. While these results are general, a concrete example considered is that of fermionic Gaussian unitaries and free fermion simulation. In that case, the difference between the simulation cost using the best known classical method and the quantum gate count is subquadratic in the model parameters [20]. This leaves open the question of whether there is a general variational ansatz construction that, starting from unitaries that are easy to simulate classically, leads to a larger classical-quantum time-complexity dif-

ference while keeping the loss landscape from flattening exponentially. To tackle this question, we introduce a stacked LCU variational ansatz (S-LCU): l sequentially composed LCUs, each a weighted sum of k unitaries on n qubits. Expanding the product over its l layers, the S-LCU is a linear combination of k l unitaries, implemented with only O(lkn2 ) gates, before postselection overhead. We denote the specialisation of the constituent unitaries to fermionic Gaussian unitaries by Free Fermion S-LCU (FF-S-LCU), an ansatz whose explicit decomposition is a linear combination of k l fermionic Gaussian unitaries that is implementable using a polynomial number of quantum gates. This count upper-bounds the fermionic Gaussian rank of the system and sets the cost of direct classical simulation [20]. Because it grows exponentially with the number of layers l, the model becomes increasingly costly for this classical method while remaining efficient to prepare with a polynomial-sized quantum circuit. Our analysis rests on a diagrammatic technique for computing the first and second moments of the S-LCU expectation value, requiring only that the unitaries be drawn from a continuous finite-dimensional unitary representation of a compact group G. Its central object is the Gram matrix of the commutant of each S-LCU layer. This provides a systematic means of constructing ansätze tailored to the complexity and trainability requirements of a particular application and hardware platform. Specialising to the FF-S-LCU, we derive a variance bound of Ω(1/(n k 3l )) for traceless quadratic observables. For a fixed number of layers, this rules out exponential cost concentration in n and k. More generally, Appendix F gives the variance for an arbitrary initial state and observable. The number of layers l is therefore a tunable knob, yielding an arbitrarily high polynomial difference between the direct classical simulation cost and the quantum gate count, at a corresponding cost in trainability.

2 We begin by describing the S-LCU construction in Section II and the notation we use. We then analyse the model’s loss landscape by characterising the first and second moments (Sections III A and III B) of the cost function when the unitaries are drawn from a representation of a compact group. In Section IV, we specialise to the FF-S-LCU, analyse its trainability and classical simulation complexity, and provide asymptotic bounds on cost concentration. We close by discussing algorithmic precedents for sequential LCUs and staged strategies for training the S-LCU.

III.

THE S-LCU LOSS LANDSCAPE

Barren plateaus were originally defined as an exponential decay of the variance of the gradient with system size [9]. It is now common practice to reason instead about the variance of the cost function itself, since exponential cost concentration is known, under the conditions of Ref. [21], to accompany a barren plateau. We therefore use the variance of the expectation value m as a trainability diagnostic: Var[m] = Var[tr(ρO)] = E[tr(ρO)2 ] − E[tr(ρO)]2 . (7)

II.

S-LCU DESCRIPTION

The linear combination of unitaries (LCU) [13] is a construction used to prepare a weighted sum of unitary operations, L=

k X

ai Ui ,

k X

ai ≥ 0,

i=1

ai = 1.

(1)

i=1

The S-LCU is obtained by the sequential composition of multiple independent LCUs, and prepares the state †

ρ := Aρ0 A , A := Ll · · · L2 L1 ,

Lj :=

k X

aji Uij ,

i.i.d.

(8)

i.i.d.

(9)

Uij ∼ Haar(G),

(2)

aj ∼ Dir(1, . . . , 1),

(3)

with the two sampled families mutually independent. In the following sections, we derive expressions for the first and second moments of the S-LCU expectation value in terms of the commutants of the group from which the unitaries are drawn. We specialise to fermionic Gaussian unitaries in Section IV.

i=1

where ρ0 is a normalised input state, aji ≥ 0 with Pk j j i=1 ai = 1 for every j, and Ui is a unitary acting on n qubits. Since an LCU is not generally unitary, ρ may be subnormalised, with a postselection success probability given by ps := tr(ρ). For a Hermitian observable O, we analyse the unnormalised expectation value m := tr(ρO).

Calculating these expectations requires defining the distributions over which the moments are taken. For the unitaries, we follow the standard practice of using the Haar distribution over the relevant group [9, 10]. The LCU coefficients are drawn independently for each layer from the uniform Dirichlet distribution. These choices are maximum entropy distributions, uniform over the parameterised spaces, and model bias-free initialisation strategies used by practitioners without further priors on the target circuit. Thus,

A.

First Moment

(4) The first moment of the expectation value m is

When deriving our results, we make extensive use of diagrammatic notation, an introduction to which can be found in Appendix A. We define Vij := aji Uij . A single LCU may then be denoted as Vij i

:=

k X

Vij ,

(5)

i=1

where the shaded box indicates a sum, indexed by i. Using this notation, the state prepared by the S-LCU may be represented as ρ=

† Vil i

...

† Vi1 i

ρ0

Vi1 i

...

Vil

.

(6)

i

Note that each layer of the S-LCU gets its own shaded box, indicating that the sum is taken independently over each layer.


 E[m] = E[tr(ρO)] = E tr Aρ0 A† O .

(10)

Utilising the diagrammatic representation of the S-LCU from the previous section, the first moment may be restated as # " † 1 l l† 1 . . . . . . ρ 0 V O Vi Vi Vi i E . i

i

i

i

(11) Defining the Bell state |Φ⟩ = can use the identity

P

x∈{0,1}n |x⟩ ⊗ |x⟩, we

tr Aρ0 A† O = ⟨Φ|(O ⊗ I)(A ⊗ A∗ )(ρ0 ⊗ I)|Φ⟩

(12)

to push linear operators around the trace. See Appendix A for more details. This yields the following dia-

3 gram, equivalent to Eq. (11),  ρ0 Vi1 . . . Vil   i i E  ∗ ∗ Vi1 . . . Vil i

 O

  . 

(13)

i

The Haar distribution is invariant under multiplication by group elements, which forces any term in which a unitary U appears without a matching U ∗ to vanish under the expectation. Lemma 1. Let U ∼ Haar(G) and a ̸= b be nonnegative integers. If the representation of the compact group G contains a scalar element ωI satisfying ω a−b ̸= 1, then h i E U ⊗a ⊗ U ∗ ⊗b = 0. (14) A proof is provided in Appendix C. Note that the global phase of each LCU term is physically relevant, and controls interference between terms. It is possible to augment each term with an independent U (1) phase group, i.e. each unitary can be sampled from the Haar(U (1)×G) distribution instead. The balanced moments of the augmented ensemble coincide with those of Haar(G), since the phases cancel, whereas all unbalanced moments vanish by the lemma above. Apart from compactness, the calculations below do not assume a particular group. We now apply Lemma 1 to the unitary expansion represented by Eq. (13), yielding   1 l ρ0 Vi . . . Vi O     E (15) ,   ∗ ∗ Vi1 . . . Vil i

U∗

To evaluate this remaining Haar average, we use the vectorised moment operator. It is the orthogonal projector onto the tensors invariant under the group action, so expressing it in terms of a basis of the commutant converts the unitary integral directly into overlaps with the initial state and observable. This is closely related to the Weingarten calculus construction discussed in Appendix B. The t-th vectorised moment operator is given by M

(t)

:=

h

E

U ∼Haar(G)

U

⊗t

⊗U

∗ ⊗t

i

=

dt X

|Pit ⟩⟩⟨⟨Pit |, (20)

i=1

where {Pit }i is a Hermitian Hilbert-Schmidt-orthonormal basis for the t-th order commutant of G, of dimension dt [22]. The vectorisation operation |·⟩⟩ is defined in Appendix A. Applying this identity at first order, and the formula for the moment of the Dirichlet-distributed co2 efficients, E[a2i ] = k(k+1) (Appendix G), gives  E[tr(ρO)] =  =

2 k+1

l X d1

2 k+1

l X d1

ρ0

Pi1

Pi1

O

i=1

tr(ρ0 Pi1 ) tr(OPi1 ).

(21)

i=1

i

where the single shaded box per layer indicates shared indexing of the sum for both matrices in each column. The coefficients and unitaries are independent. Factoring out the common coefficient moment from the sum therefore yields   1 l . . . ρ 0 U O U i i     E[a2i ]l · E  (16) .   1∗ . . . l∗ Ui Ui i

i

The product of independent Haar-distributed unitaries is itself Haar-distributed. Indeed, for any fixed g ∈ G, gUl Ul−1 · · · U1

Each of the k l ways of choosing one unitary per layer therefore yields the same single-unitary Haar average, giving   ρ0 U O   E[a2i ]l · k l · E  (19) .

Haar

=

Ul Ul−1 · · · U1 ,

(17) Haar

by the left-invariance of the Haar distribution (gUl = Ul ). This in turn allows us to reduce the expectation under the Haar distribution of multiple layers to the expectation of a single layer as Haar

Ul Ul−1 · · · U1 = U.

(18)

Result 1: The S-LCU first moment For a compact unitary group G satisfying Lemma 1, with observable O, and initial state ρ0 , the first moment is  E[tr(ρO)] =

2 k+1

l X d1

tr ρ0 Pi1 tr OPi1 ,

i=1

(22) a sum over an orthonormal basis for the first-order commutant {Pi1 } of G. The dependence on ρ0 , O, and the group is therefore identical to the Haar average of a single unitary from G: stacking changes only the overall factor (2/(k + 1))l . The second moment is more involved, and is tackled below.

B.

Second Moment

Computing the variance of the cost function requires computing its second moment, given in the S-LCU case

4 by h
2 i . E[tr(ρO)2 ] = E tr Aρ0 A† O

(23)

distinct unitary label must occur equally often in unconjugated and conjugated form. This leaves exactly three index-identification patterns,

Restating diagrammatically,    E  

O

Vil

†

. . . V 1† i

i

O

Vil

ρ0

i

†

Vil

i

. . . V 1† i

i

...

Vi1

ρ0

i

...

Vi1

i

Vil

i

     

∗

,

∗

      E      

Vi1

...

Vil

i

. . . Vil∗

i

ρ0

O

i

∗ Vi1

i

Vi1

...

Vil

i

O

i

∗ Vi1

...

∗ Vil

i

k X

      .      

i ̸= i′ ,

(30)

respectively. The paired-slice picture above may be restated as below, by permuting wires appropriately, (25) ∗

∗

. ∗

i

Vij1 ⊗ (Vij2 )∗ ⊗ Vij3 ⊗ (Vij4 )∗ .

(26)

i1 ,i2 ,i3 ,i4 =1

∗

∗

∗

(31) Recall that these diagrams contain not just the unitaries, but also the coefficients drawn from the uniform Dirichlet distribution. The slice average is E[Sj ] = k · E[a4i ]Mπ1 + k(k − 1) · E[a2i a2i′ ]Mπ2

Independence between layers gives (27)

ρ0

ρ0

       E      

Vij i

Vij

∗

i

Vij i

Vij

∗

Mπ 1 =

l O

             

(32)

Note that the counts here are order-sensitive. The associated moment operators can be written compactly as

Consequently, 

i ̸= i′ .

+ k(k − 1) · E[a2i a2i′ ]Mπ3 ,

E[Sl · · · S1 ] = E[Sl ] · · · E[S1 ] = E[Sj ]l .

E[tr(ρO)2 ] =

∗

π1 : Ui ⊗ Ui∗ ⊗ Ui ⊗ Ui∗ , π2 : Ui ⊗ Ui∗ ⊗ Ui′ ⊗ Ui∗′ , π3 : Ui ⊗ Ui∗′ ⊗ Ui′ ⊗ Ui∗ ,



This diagram consists of l vertical slices, independently and identically distributed, since they correspond to different LCUs. We denote the operator formed by layer j by Sj :=

∗

(29)

where the shaded boxes enclose like unitaries and their corresponding conjugates. These represent the pairings

By pushing operators around the trace, as in the firstmoment calculation, we equivalently obtain ρ0

∗

i

(24)



∗

d2 X

|Piπ1 ⟩⟩⟨⟨Piπ1 |,

(33)

i=1

Mπq =

d1 X

π

π

|Pijq ⟩⟩⟨⟨Pijq |,

q ∈ {2, 3},

(34)

i,j=1

(28) O

i

The expansion (26) introduces four indices: one choice of unitary and one choice of conjugate unitary in each of the two copies of the trace. All four are drawn from the same j-th LCU, but the sums choose their indices independently before the Haar average imposes pairings. Applying Lemma 1 to the four factors in a slice, every

where each expression is a sum over vectorised orthonormal basis elements,

|Piπ1 ⟩⟩ :=

Pi2

Pi1

|Pijπ2 ⟩⟩ :=

Pj1

(35)

Pi1

|Pijπ3 ⟩⟩ :=

Pj1

.

(36)

5 We concatenate these orthonormal basis elements into an ordered, overcomplete frame for the commutant of a single S-LCU layer,     (37) P := [Piπ1 ]i || Pijπ2 ij || Pijπ3 ij ,

{cµ , cν } = 2δµν I. Up to a global phase, a paritypreserving fermionic Gaussian unitary has the form

of D := d2 + 2 · d21 elements. Similarly, the associated coefficients are arranged in a vector p. The ordering need only be consistent between the frame and the coefficient vector. The expectation from Eq. (28) may then be reexpressed as  l ρ0 O D X    E[tr(ρO)2 ] = pi Pi Pi  .  ρ0   O i=1

(43)

2n

U = e−iH ,

H=

V :=

D X

4

|i⟩⟨⟨Pi |

∈ CD×d

(39)

∈ RD×D .

(40)

i=1

Db := Diag(p) =

D X

pi · |i⟩⟨i|

i=1

Defining |A, A⟩⟩ := |A⟩⟩ ⊗ |A⟩⟩, we have †

2

l

E[tr(ρO) ] =⟨⟨O, O|(V Db V ) |ρ0 , ρ0 ⟩⟩ †

† l−1

=⟨⟨O, O|V (Db V V )

h ∈ so(2n),

where so(2n) := {h ∈ R2n×2n : hT = −h}

U † cµ U =

2n X

Rµν cν ,

R ∈ SO(2n).

Db V |ρ0 , ρ0 ⟩⟩. (42)

Observe that G := V V † , with entries Gij = ⟨⟨Pi |Pj ⟩⟩, is the Gram matrix containing all inner products between elements of the ordered frame. Result 2: The S-LCU second moment For a compact unitary group G satisfying Lemma 1, with observable O, and initial state ρ0 , the S-LCU second moment takes the form E[tr(ρO)2 ] = ⟨⟨O, O|V † (Db G)l−1 Db V |ρ0 , ρ0 ⟩⟩, consisting of the following elements,

(45)

ν=1

The matrix R records how U rotates the 2n Majorana operators, whereas U acts on the 2n -dimensional Fock space. The unitaries U and −U induce the same rotation, so the correspondence is two-to-one. The group that retains this sign is Spin(2n), the double cover of SO(2n). Relative global phases between LCU terms affect their interference. Accordingly, the FF-S-LCU uses the Haar distribution on GFF := U (1) × Spin(2n),

(41)

(44)

is the Lie algebra of SO(2n), consisting of real antisymmetric matrices. The conjugation action of U is linear,

(38) We define two matrices to aid in the calculation of this expression,

i X hµν cµ cν , 4 µ,ν=1

(46)

where the U (1) factor enforces Lemma 1. Appendix D collects the remaining free fermion definitions and identities. Fermionic Gaussian unitaries can be classically simulated in time polynomial in the number of qubits [23, 24]. Previous work has connected their polynomial trainability to their dynamical Lie algebra [12, 25–29]. They are a natural choice for the S-LCU: a single linear combination already increases their expressivity at a quantifiable cost in trainability [19], and stacking allows that trade-off to be tuned further. In this section we analyse concretely the variance of the loss landscape of the FF-S-LCU: an S-LCU composed of fermionic Gaussian unitaries. We begin by reviewing the first and second order commutants of the underlying group.

• the Gram matrix G = V V † ; • the initial state projections V |ρ0 , ρ0 ⟩⟩; and • the observable projections V |O, O⟩⟩. In the following section, we consider fermionic Gaussian unitaries and derive the Gram matrix G explicitly. IV.

FERMIONIC GAUSSIAN UNITARIES

1.

First-Order Commutant

The first-order commutant of the group of fermionic Gaussian unitaries is well-studied [26, 30, 31]. It is spanned by the orthonormal basis     I P 1 1 √ ,√ , d = 2n . = √ · ,√ · P d d d d (47) Here P is the parity operator, defined as

For n fermionic modes, encoded on n qubits, let c1 , . . . , c2n be Hermitian Majorana operators satisfying

P := Z ⊗n = (−i)n · c1 c2 . . . c2n .

(48)

6 2.

We plot only the elements on or below the principal diagonal, with G† = G specifying the rest.

Second-Order Commutant

We use the definition of the second-order commutant from Diaz et al. [26], which has also been characterised in Refs. [30, 31]: )2n

(

)2n

(

Q0κ

Q1κ

∪

,

κ=0

(49)

κ=0

where Q0κ

X

:= Nκ

cs ⊗ cs ,

:= iκ mod 2 ·

Result 2 contracts the layer operator with the initial state and observable. The initial state boundary vector needed for this contraction has entries 1h tr(ρ0 )2 , tr(P ρ0 ) tr(ρ0 ), tr(P ρ0 ) tr(ρ0 ), d tr(P ρ0 )2 , tr(ρ20 ), tr(P ρ20 ), tr(P ρ20 ), i eκ (ρ0 ), Ceκ (ρ0 ) tr(P ρ0 P ρ0 ), P (54) where

Q0κ

P

,

(51)

[2n]

where κ := {s ⊆ [2n] : |s| = κ} and Nκ :=  q −1
2n 2n . Algebraically, the second line is Q1κ = κ

eκ (ρ0 ) := tr(ρ⊗2 Q0 ), P κ 0

A.

3.

Gram Matrix

The full Gram matrix G for the FF-S-LCU, expressed in the ordered frame P defined in Eq. (37), is given below. We define the following shorthand: q

2n κ

, σκ := (−1)⌊κ/2⌋ , d ϵκ := (−1)κ(κ+1)/2 .

(52) (53)

where d = 2n . Derivations for each scalar element of the matrix are provided in Appendix E. Q0κ

Q0κ

Q1κ

P

Q1κ

P P

P

P

P P

P

δκ,κ′

·

·

·

·

·

·

·

·

·

0

δκ,κ′

·

·

·

·

·

·

·

·

δ0,κ

0

1

·

·

·

·

·

·

·

0

(−1)n δ2n,κ

0

1

·

·

·

·

·

·

0

δ0,κ

0

0

1

·

·

·

·

·

(−1)n δ2n,κ

0

0

0

0

1

·

·

·

·

σκ νκ

0

1 d

0

0

1 d

1

·

·

·

0

iκ mod 2 σκ νκ

0

1 d

1 d

0

0

1

·

·

0

κ mod 2

0

1 d

1 d

0

0

0

1

·

1 d

0

0

1 d

0

0

0

1

P

P P P

P

P

i

ϵκ ν κ

P P

ϵκ ν κ

0

1 Ceκ (ρ0 ) := tr(ρ⊗2 0 Qκ ).

(55)

The observable boundary vector V |O, O⟩⟩, which closes the other side of the same contraction, has the same form with ρ0 replaced by O.

iκ mod 2 (I ⊗ P )Q0κ .

νκ :=

Projections of Initial State and Observable

(50)

s∈([2n] κ )

Q1κ

4.

FF-S-LCU Loss Landscape

Using the Gram matrix G, coefficient matrix Db , and the initial state and observable vectors derived above, the variance Var[m] of the FF-S-LCU may be expressed as a function of the number of qubits n, LCU term count k, and depth l. Theorem 1. For a traceless quadratic observable O with tr(O2 ) = d, ρ0 = |0⟩⟨0|⊗n , n ≥ 3, with the variance taken over the joint distribution Haar(GFF )⊗kl × Dir(1, . . . , 1)⊗l , Var[m] ≥

1 2n − 1



24 (k + 1)(k + 2)(k + 3)

l

  1 ∈Ω . n k 3l

A proof of this bound is provided in Appendix F (via Theorem 2). The mean vanishes because tr(O) = tr(P O) = 0, so the variance equals the second moment used in the proof. Appendix F explains how the same methods extend to arbitrary initial states and observables. For a fixed number of layers, the lower bound decreases polynomially in both the number of qubits and the number of terms in each LCU, ruling out exponential cost concentration in either variable. This conclusion does not apply if l itself grows with the problem size. In Figure 1 we plot the exact variance against the number of qubits for different values of l.

7 l=1

100

l=2

l=3 k=1 k=5 k = 10 k = 15 k = 20 k = 25

10−2

Var[m]

10−4

k = 30 k = 35 k = 40 k = 45 k = 50

10−6 10−8 10−10 10−12 2

6

10

14 18 n (qubits)

22

26

30

2

6

10

14 18 n (qubits)

22

26

30

2

6

10

14 18 n (qubits)

22

26

30

FIG. 1: Variance of the unnormalised expectation value Var[m] as a function of the number of qubits n, for l ∈ {1, 2, 3} and varying LCU term count k. Observable O = Z1 , initial state ρ0 = |0⟩⟨0|⊗n . Constant simulation cost (k l = 64) k = 64, l = 1 k = 8, l = 2 k = 4, l = 3 k = 2, l = 6 k = 16, l = 1 k = 24, l = 1 k = 32, l = 1

Var[m]

10−3

10

−4

10−5

value is O(r2 n3 ), since each of the r2 cross-terms requires a phased overlap computation in O(n3 ) classical operations. For the S-LCU model, expanding the product of l layers yields a superposition of k l fermionic Gaussian states. The resulting state has fermionic Gaussian rank at most k l , and direct pairwise evaluation of these overlaps yields a complexity

10−6 5

10

15 n (qubits)

20

25

30

FIG. 2: Variance for FF-S-LCU configurations with constant expanded term count k l = 64. Solid lines: (k, l) pairs with k l = 64. Dashed lines: single-layer models with comparable k. Distributing a fixed term count over more layers produces a larger variance in these examples. Figure 2 plots the variance as a function of the number of qubits for various (k, l) combinations with k l = 64. For the configurations shown, keeping the expanded term count (and therefore direct simulation cost) constant while increasing the number of layers reduces cost concentration.

tdirect ∈ O (l + k 2l )n3 = O(k 2l n3 )

(56)

We are not aware of a method that exploits the stacked structure to lower the exponent in k. On the quantum side, each layer requires implementing an LCU of k fermionic Gaussian unitaries, realisable in O(kn2 ) twoqubit gates [34]. Before accounting for postselection or amplitude amplification, the total gate count for l layers is therefore

Ngates ∈ O(lkn2 ). B.

(k ≥ 2).

(57)

Time Complexity: Classical and Quantum

Free fermion dynamics can be classically simulated in polynomial time [23, 24]. A single fermionic Gaussian unitary on n qubits is fully characterised, up to phase, by an SO(2n) rotation matrix, and overlaps between fermionic Gaussian states can be computed in O(n3 ) time, provided the relative phase is tracked. Recent work has extended efficient classical simulation to states of bounded Gaussian rank, where the rank is the minimum number of fermionic Gaussian states in a superposition representing the state [20, 32, 33]. Given an explicit decomposition into r Gaussian states, the direct cost of computing a fixed-complexity local expectation

This runtime analysis does not account for the postselection cost associated with implementing non-unitary operations. Depending on the application, the runtime for obtaining successful shots may scale with the inverse success probability. We comment on this in Section V. The variance calculation, however, does account for the postselection probability, since the subnormalisation of ρ = Aρ0 A† is precisely this probability and m = tr(ρO) is weighted by it. No separate correction is therefore required for the variance bound.

8 Result 3: The FF-S-LCU tradeoff For a traceless quadratic observable O with tr(O2 ) = d, ρ0 = |0⟩⟨0|⊗n , and n ≥ 3, the unnormalised cost variance satisfies   1 Var[m] ∈ Ω . (58) n k 3l For fixed l, direct classical simulation using term expansion has the upper bound
tdirect ∈ O (l + k 2l )n3 , (59) while the quantum circuit’s gate count remains linear in k before postselection overhead, scaling as
Ngates ∈ O lkn2 . (60)

cording to a gradient or residual criterion, assigned a small nonzero coefficient, and then optimised together with the existing terms. This is analogous in spirit to adaptive variational methods that grow an ansatz one operator at a time [41], but here the added resource is an LCU term; the expanded term count upper-bounds the fermionic Gaussian rank. These warm-start strategies keep the effective depth and term count small when initialising. They do not, however, constitute a guarantee during optimisation: our present bounds concern independent uniform initialisation (as captured by the Haar/Dirichlet distributions), and analysing structured, correlated initialisations and the gradients encountered throughout this staged optimisation remains an open problem.

VI.

CONCLUSION

Sequentially composed LCUs already occur in several quantum algorithms. Truncated-Taylor Hamiltonian simulation divides an evolution into short-time segments and implements every segment by an LCU followed by robust oblivious amplitude amplification [35]. The truncated-Dyson algorithm uses the same segmented structure for time-dependent Hamiltonians [36], while multiproduct formula simulation implements each shorttime step as an LCU of product formulas and repeats the amplified step [13, 37]. More abstractly, LCU constructions are a standard route to block encodings, and block encodings can be composed to encode matrix products [38, 39]. These examples are not variational ansätze, but they establish that a product of LCU-derived maps is a natural circuit primitive. They also illustrate an important implementation condition: the cited algorithms avoid postselecting each individual LCU. Instead, they keep the construction coherent and control the normalisation, or make each segment approximately unitary and amplify it before applying the next segment. The layered parameterisation suggests alternatives to the fully random initialisation analysed above. One may first optimise a single LCU, append a new layer initialised near the identity, and then resume optimisation. Incrementally increasing circuit depth has been observed to retain larger gradients than training the full circuit from scratch [40]. The same idea can be applied term-wise within a layer. Candidate unitaries can be added ac-

Our primary contribution is the S-LCU, a stacked linear combination of unitaries variational ansatz, together with a diagrammatic framework for characterising its loss landscape when the unitaries are drawn from a compact group G. This gives practitioners a systematic way to tailor the complexity-trainability balance to their application and hardware. Combining this framework with known results on the linear and quadratic commutants of fermionic Gaussian unitaries, we bound the rate of decay of the loss-landscape variance of the FF-S-LCU for traceless quadratic observables. For a fixed number of layers, the variance is bounded from below by an inverse polynomial in both the number of terms per layer and the number of qubits. We additionally derive an exact expression for arbitrary normalised initial states and Hermitian observables in Appendix F. Increasing the number of layers yields an arbitrarily high polynomial separation (in the per-layer term count k) between the simulation cost O(k 2l n3 ) of the best classical method known to us, and the coherent gate count O(lkn2 ), while the variance lower bound decreases at the corresponding rate k −3l . Note, however, that while the postselection probability is accounted for in our variance analysis, it needs to be controlled in practical applications. Further work would be necessary to realise a practical separation from classical methods, which must contend with postselection probabilities, possible classical algorithms that exploit the stacked structure, gate execution times on real quantum devices, and sources of noise such as decoherence and imprecise readout.

[1] M. Benedetti, E. Lloyd, S. Sack, and M. Fiorentini, Parameterized quantum circuits as machine learning models, Quantum Science and Technology 4, 043001 (2019).

[2] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Quantum machine learning, Nature 549, 195 (2017).

V.

DISCUSSION: PRACTICAL SETTINGS AND TRAINING STRATEGIES

9 [3] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, A variational eigenvalue solver on a photonic quantum processor, Nature Communications 5, 4213 (2014). [4] J. R. McClean, J. Romero, R. Babbush, and A. AspuruGuzik, The theory of variational hybrid quantumclassical algorithms, New Journal of Physics 18, 023023 (2016). [5] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. G. Booth, and J. Tennyson, The variational quantum eigensolver: A review of methods and best practices, Physics Reports 986, 1 (2022). [6] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm (2014), arXiv:1411.4028 [quant-ph]. [7] K. Blekos, D. Brand, A. Ceschini, C.-H. Chou, R.-H. Li, K. Pandya, and A. Summer, A review on Quantum Approximate Optimization Algorithm and its variants, Physics Reports 1068, 1 (2024). [8] A. Abbas, A. Ambainis, B. Augustino, et al., Challenges and opportunities in quantum optimization, Nature Reviews Physics 6, 718 (2024). [9] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Barren plateaus in quantum neural network training landscapes, Nature Communications 9, 4812 (2018). [10] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Cost function dependent barren plateaus in shallow parametrized quantum circuits, Nature Communications 12, 1791 (2021). [11] M. Larocca, S. Thanasilp, S. Wang, K. Sharma, J. Biamonte, P. J. Coles, L. Cincio, J. R. McClean, Z. Holmes, and M. Cerezo, Barren plateaus in variational quantum computing, Nature Reviews Physics 7, 174–189 (2025). [12] M. Cerezo, M. Larocca, D. Garcı́a-Martı́n, N. L. Diaz, P. Braccia, E. Fontana, M. S. Rudolph, P. Bermejo, A. Ijaz, S. Thanasilp, E. R. Anschuetz, and Z. Holmes, Does provable absence of barren plateaus imply classical simulability?, arXiv preprint arXiv:2312.09121 (2023). [13] A. M. Childs and N. Wiebe, Hamiltonian Simulation Using Linear Combinations of Unitary Operations, Quantum Information and Computation 12, 10.26421/qic12.11-12 (2012). [14] N. Khatri, G. Matos, L. Coopmans, and S. Clark, Quixer: A quantum transformer model, arXiv preprint arXiv:2406.04305 (2024). [15] J. Heredge, M. West, L. Hollenberg, and M. Sevior, Nonunitary quantum machine learning, Phys. Rev. Appl. 23, 044046 (2025). [16] L. Coopmans and M. Benedetti, On the sample complexity of quantum Boltzmann machine learning, Communications Physics 7, 274 (2024). [17] H. Yao, X. Liu, M. Jing, G. Li, and X. Wang, LCQNN: Linear combination of quantum neural networks, arXiv preprint arXiv:2507.02832 (2025). [18] B. Coyle, S. Raj, N. Mathur, E. A. Cherrat, N. Jain, S. Kazdaghli, and I. Kerenidis, Training-efficient density quantum machine learning, npj Quantum Information 11, 172 (2025). [19] N. Khatri, S. Zohren, and G. Matos, Trainability of parametrised linear combinations of unitaries, arXiv preprint arXiv:2506.22310 (2025). [20] B. Dias and R. König, Classical simulation of non-

Gaussian fermionic circuits, Quantum 8, 1350 (2024). [21] A. Arrasmith, Z. Holmes, M. Cerezo, and P. J. Coles, Equivalence of quantum barren plateaus to cost concentration and narrow gorges, Quantum Science and Technology 7, 045015 (2022). [22] A. A. Mele, Introduction to Haar measure tools in quantum information: A beginner’s tutorial, Quantum 8, 1340 (2024). [23] L. G. Valiant, Quantum computers that can be simulated classically in polynomial time, in Proceedings of the Thirty-Third Annual ACM Symposium on Theory of Computing (ACM, 2001) pp. 114–123. [24] B. M. Terhal and D. P. DiVincenzo, Classical simulation of noninteracting-fermion quantum circuits, Physical Review A 65, 032325 (2002). [25] G. Matos, C. N. Self, Z. Papić, K. Meichanetzidis, and H. Dreyer, Characterization of variational quantum algorithms using free fermions, Quantum 7, 966 (2023). [26] N. Diaz, D. Garcı́a-Martı́n, S. Kazi, M. Larocca, and M. Cerezo, Showcasing a barren plateau theory beyond the dynamical lie algebra, arXiv preprint arXiv:2310.11505 (2023). [27] E. Fontana, D. Herman, S. Chakrabarti, N. Kumar, R. Yalovetzky, J. Heredge, S. H. Sureshbabu, and M. Pistoia, Characterizing barren plateaus in quantum ansätze with the adjoint representation, Nature Communications 15, 7171 (2024). [28] E. Kökcü, T. Steckmann, Y. Wang, J. K. Freericks, E. F. Dumitrescu, and A. F. Kemper, Fixed depth hamiltonian simulation via cartan decomposition, Physical Review Letters 129, 070501 (2022). [29] R. Wiersema, E. Kökcü, A. F. Kemper, and B. N. Bakalov, Classification of dynamical lie algebras for translation-invariant 2-local spin systems in one dimension, npj Quantum Information 10, 110 (2024). [30] P. Sierant, X. Turkeshi, and P. S. Tarabunga, Theory of the matchgate commutant, arXiv preprint arXiv:2603.12392 (2026). [31] M. Lastres and S. Moudgalya, Geometry of free fermion commutants, arXiv preprint arXiv:2604.05031 (2026). [32] O. Reardon-Smith, M. Oszmaniec, and K. Korzekwa, Improved simulation of quantum circuits dominated by free fermionic operations, Quantum 8, 1549 (2024). [33] J. Cudby and S. Strelchuk, Gaussian decomposition of magic states for matchgate computations, arXiv preprint arXiv:2307.12654 (2024). [34] R. Babbush, C. Gidney, D. W. Berry, N. Wiebe, J. McClean, A. Paler, A. Fowler, and H. Neven, Encoding electronic spectra in quantum circuits with linear T complexity, Physical Review X 8, 041015 (2018). [35] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Simulating Hamiltonian dynamics with a truncated Taylor series, Physical Review Letters 114, 090502 (2015). [36] M. Kieferová, A. Scherer, and D. W. Berry, Simulating the dynamics of time-dependent Hamiltonians with a truncated Dyson series, Physical Review A 99, 042314 (2019). [37] G. H. Low, V. Kliuchnikov, and N. Wiebe, Wellconditioned multiproduct Hamiltonian simulation (2019), arXiv:1907.11679 [quant-ph]. [38] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, Quantum singular value transformation and beyond: Exponential improvements for quantum matrix arithmetics, in Pro-

10 ceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing (2019) pp. 193–204. [39] S. Chakraborty, A. Gilyén, and S. Jeffery, The power of block-encoded matrix powers: Improved regression techniques via faster Hamiltonian simulation, in 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), Vol. 132 (2019) pp. 33:1–33:14. [40] A. Skolik, J. R. McClean, M. Mohseni, P. van der Smagt, and M. Leib, Layerwise learning for quantum neural networks, Quantum Machine Intelligence 3, 5 (2021). [41] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, An adaptive variational algorithm for exact molecular simulations on a quantum computer, Nature Communications 10, 3007 (2019). [42] R. Penrose, Applications of negative dimensional tensors, in Combinatorial Mathematics and its Applications, edited by D. J. A. Welsh (Academic Press, 1971) pp. 221–244. [43] S. Abramsky and B. Coecke, A categorical semantics of quantum protocols, in Proceedings of the 19th Annual

IEEE Symposium on Logic in Computer Science (LICS) (IEEE, 2004) pp. 415–425. [44] B. Coecke and A. Kissinger, Picturing Quantum Processes: A First Course in Quantum Theory and Diagrammatic Reasoning (Cambridge University Press, 2017). [45] J. Biamonte and V. Bergholm, Tensor networks in a nutshell, arXiv preprint arXiv:1708.00006 (2017). [46] D. Weingarten, Asymptotic behavior of group integrals in the limit of infinite rank, Journal of Mathematical Physics 19, 999 (1978). [47] B. Collins, Moments and cumulants of polynomial random variables on unitary groups, the Itzykson-Zuber integral, and free probability, International Mathematics Research Notices 2003, 953 (2003). [48] B. Collins and P. Śniady, Integration with respect to the Haar measure on unitary, orthogonal and symplectic group, Communications in Mathematical Physics 264, 773 (2006).

Appendix A: Reading the Diagrams

All diagrams in the paper are tensor networks, with minimal additional decoration to represent sums. This follows a long tradition of diagrammatic representations in quantum computing, finding use in tensor networks and categorical quantum mechanics [42–45]. Linear algebraic expressions can be represented using diagrams composed of wires and boxes. Wires carry finitedimensional Hilbert spaces, while boxes represent linear maps. For example, A : Cn → Cm

n /

←→

A

m / .

(A1)

When the system acted upon is unambiguous, the system size is omitted from the notation: .

A

(A2)

Linear maps may be composed sequentially, BA =

A

B

,

(A3)

.

(A4)

or in parallel, A

A⊗B = B

The unnormalised maximally entangled vector is represented by a “cup”: X |Φ⟩ = |x⟩ ⊗ |x⟩ = .

(A5)

x∈{0,1}n

The transpose of a linear map is represented by bending wires as AT =

A

.

(A6)

Complex conjugation and the adjoint are not represented by any change to the diagram: the box is drawn identically, with only its label decorated, as A∗ or A† . The trace of a linear map is represented by tr(A) =

A

.

(A7)

11 The vectorisation of a linear map is A

|A⟩⟩ := (A ⊗ I) |Φ⟩ =

.

(A8)

We denote by |A, B⟩⟩ the paired vectorisation A

|A, B⟩⟩ := |A⟩⟩ ⊗ |B⟩⟩ =

B

.

(A9)

The Hilbert-Schmidt inner product of two linear maps is given by B

⟨A, B⟩HS = ⟨⟨A|B⟩⟩ = tr(A† B) =

A†

.

(A10)

Tensor network notation does not natively offer a suitable standard for sums of linear maps, for which we introduce shaded boxes: X Ai Ai = . (A11) i

i

Two separate shaded boxes in the same diagram indicate independent sums, Ai

X

Ai ⊗

i

X

i

Bj =

j

,

(A12)

Bj j

while a single box spanning multiple linear maps indicates a shared index for the sum, Ai

X i

Ai ⊗ B i =

.

(A13)

Bi i

Appendix B: Relation to the Weingarten Calculus

Our method bears resemblance to the Weingarten calculus [46–48], which also computes Haar averages from the Gram matrix of a group commutant [22]. Our method differs in two respects. The Weingarten calculus builds its Gram matrix from a basis of the commutant, whereas our diagrammatic pairings give an overcomplete frame. It then inverts this Gram matrix, whereas we instead raise the Dirichlet-weighted Gram matrix Db G to a power, one factor for each layer. Since our frame is overcomplete, its Gram matrix is singular; this does not matter here because no inverse is required. Appendix C: Vanishing of Unbalanced Haar Moments

Lemma 1. Let U ∼ Haar(G) and a ̸= b be nonnegative integers. If the representation of the compact group G contains a scalar element ωI satisfying ω a−b ̸= 1, then h i E U ⊗a ⊗ U ∗ ⊗b = 0. (14) Proof. By left-invariance of the Haar measure, for any fixed g ∈ G, h i h i ⊗b E U ⊗a ⊗ U ∗ ⊗b = E (gU )⊗a ⊗ (gU )∗ .

(C1)

12 Taking g = ωI gives E[U ⊗a ⊗ U ∗ ⊗b ] = ω a−b E[U ⊗a ⊗ U ∗ ⊗b ].

(C2)

E[U ⊗a ⊗ U ∗ ⊗b ] = 0.

(C3)

Since ω a−b ̸= 1,

Appendix D: Free Fermion Definitions and Identities

In this section we introduce the definitions and identities required to prove our results involving fermionic Gaussian unitaries. 1.

Definitions

A free fermion system on n modes is described by creation and annihilation operators a†j , aj (j = 1, . . . , n) obeying the canonical anticommutation relations (CAR), {aj , a†k } = δjk I,

{aj , ak } = 0,

(D1)

where {A, B} := AB + BA is the anticommutator. It is convenient to work instead with the 2n Hermitian Majorana operators c2j−1 := aj + a†j ,

c2j := i (a†j − aj ),

(D2)

which obey c2µ = I,

{cµ , cν } = 2δµ,ν I,

(D3)

so that each cµ is its own inverse and distinct Majoranas anticommute. A free fermion system evolves under a quadratic Majorana Hamiltonian 2n

H=

i X hµν cµ cν , 4 µ,ν=1

h ∈ R2n×2n ,

hT = −h,

(D4)

and the associated unitaries e−iHt realise Spin(2n) on Fock space, inducing rotations in SO(2n) on the Majorana operators. Matchgate computation is historically connected to Valiant’s perfect-matching identities [23]. A Majorana string cs denotes multiple sequentially composed Majorana operators, indexed by an ordered subset s = {ν1 < · · · < νκ } ⊆ [2n]: cs := cν1 cν2 · · · cνκ .

κ := |s|,

(D5)

The parity operator P is the full Majorana string c1 c2 · · · c2n , including the phase that makes it equal to Z ⊗n : P := Z ⊗n = (−i)n c1 c2 . . . c2n

(D6)

2

(D7)

P = I.

Q0κ

:= Nκ

X s∈(

Q1κ Equivalently, Q1κ = iκ mod 2 (I ⊗ P )Q0κ .

cs ⊗ cs

(D8)

Q0κ

(D9)

[2n] κ

:= iκ mod 2 ·

)

P

13 2.

Identities

We make use of the following general identities about majorana operators, the parity operator, and their interaction. cs † = (−1)κ(κ−1)/2 cs = (−1)⌊κ/2⌋ cs .

(D10)

tr(cs cs ) = (−1)⌊κ/2⌋ · tr(cs cs † ) = (−1)⌊κ/2⌋ · d

(D11)

†

tr(P c[2n] ) = (−i)n tr(c[2n] c[2n] ) = (−i)n · (−1)n·(2n−1) tr(c[2n] c[2n] ) = (−i)n · (−1)n·(2n−1) · d tr(Q0κ ) = δ0,κ · Nκ · tr(Id ⊗ Id ) = δ0,κ · cs

cs

1 2 · d = δ0,κ · d d

= tr(cs cs ) = (−1)κ(κ−1)/2 tr(cs cs † ) = (−1)κ(κ−1)/2 · d P cs = (−1)κ · cs P

(D12) (D13)

(D14) (D15)

Appendix E: Gram Matrix Elements for the Fermionic Gaussian Group 1.

First-Order Terms

We evaluate gram matrix elements associated with the first-order commutant of the fermionic Gaussian group. All terms here are in {0, 1, d1 }. G3,3 : 1 · d2

=

1 1 tr(I) tr(I) = 2 · d2 = 1 2 d d

(E1)

=

1 tr(P ) tr(I) = 0 d2

(E2)

=

1 tr(I) tr(P ) = 0 d2

(E3)

=

1 tr(P ) tr(P ) = 0 d2

(E4)

G3,4 : P

1 · d2 G3,5 : 1 · d2

P

G3,6 : P

1 · d2

P

G7,7 : 1 · d2

=

1 1 tr(I) tr(I) = 2 · d2 = 1 2 d d

(E5)

14 G7,8 : P

1 · d2

=

1 tr(P ) tr(I) = 0 d2

(E6)

1 · d2

=

1 tr(I) tr(P ) = 0 d2

(E7)

=

1 tr(P ) tr(P ) = 0 d2

(E8)

1 1 1 tr(I) = 2 · d = 2 d d d

(E9)

G7,9 :

P

G7,10 : P

1 · d2

P

G3,7 : 1 · d2

=

G3,8 : P

1 · d2

=

1 tr(P ) = 0 d2

(E10)

1 · d2

=

1 tr(P ) = 0 d2

(E11)

1 1 1 tr(P 2 ) = 2 · d = d2 d d

(E12)

G3,9 :

P

G3,10 : P

1 · d2

=

P

2.

Second-Order Terms

Here we evaluate the inner products Gα,β , for all terms involving at least one second-order commutant element. The indices follow the ordering of P. G3,1 :

1 · d

Q0κ =

1 1 · tr(Q0κ ) = δ0,κ · 2 · d2 = δ0,κ d d

(E13)

15 G3,2 : Q1κ

1 · d

=

1 · tr(Q1κ ) = 0 d

(E14)

G4,1 : P

1 · d

Q0κ

P

1 = · d

Q0κ

Nκ = d

P

X

cs cs

=0

(E15)

cs ∈([2n] κ )

G4,2 : P

Q1κ

1 · d

P

1 = · d =

Q1κ

P

iκ mod 2 = · d

P

Q0κ

(E16)

iκ mod 2 · Nκ · δ2n,κ · tr(P c[2n] )2 = δ2n,κ · (−1)n d

(E17)

G5,1 :

1 · d

Q0κ

1 · d

P

Q1κ

=

1 · d

P

1 · d

P

=

P

Q0κ

=

1 · Nκ d

cs

X P

cs ∈ [2n] κ

)

Q0κ

=

(

cs

=0

(E18)

G5,2 :

1 · d

Q1κ =

P

1 · d

P

iκ mod 2 · d

P P

iκ mod 2 · tr(Q0κ ) = δ0,κ d

(E19)

G6,1 : P

1 · d

Q0κ =

P

P

Q0κ

1 = δ2n,κ 2 tr(P c[2n] )2 = δ2n,κ · (−1)n d

Q1κ

=

(E20)

G6,2 : P

1 · d

Q1κ =

P

P

1 κ mod 2 ·i · d

P

Q0κ

=0

cs

(−1)⌊κ/2⌋ = · d

(E21)

G7,1 :

1 · d

Q0κ

cs

1 = · Nκ d

X cs ∈([2n] κ )

cs

1 = · Nκ d

X cs ∈([2n] κ )

cs

s  2n κ (E22)

16 G7,2 : cs

Q1κ

1 · d

=

1 · Nκ · iκ mod 2 d

cs

X cs ∈ [2n] κ

(

=

P

)

1 · Nκ · iκ mod 2 d

cs ∈ [2n] κ

(

cs

P

X

cs

=0

) (E23)

G8,1 : P

Q0κ

1 · d

=

1 · Nκ d

P

X

cs

cs

=0

(E24)

cs ∈ [2n] κ

(

)

G8,2 : P

P

Q1κ

1 · d

κ mod 2

=

i

Q0κ

·

d

P

cs κ mod 2

=

i

Nκ

d

cs

X cs ∈([2n] κ )

iκ mod 2 (−1)⌊κ/2⌋ = d

s  2n . κ

(E25)

G9,1 :

1 · d

Q0κ =

P

1 · Nκ d

P

X

cs

cs

=0

(E26)

cs ∈ [2n] κ

(

)

=

iκ mod 2 Nκ d

G9,2 :

1 · d

Q1κ P

Q0κ

iκ mod 2 = · d

P

P

iκ mod 2 = (−1)κ mod 2 · · Nκ d

cs

X

cs

cs ∈([2n] κ )

P

X

cs

P

cs

s

2n κ

cs ∈([2n] κ )

iκ mod 2 = (−1)κ(κ+1)/2 · d

 (E27)

G10,1 : P

1 · d

P

Q0κ =

1 · Nκ d

X cs ∈([2n] κ )

P

cs

P

cs

=

(−1)κ · Nκ d

X

cs

cs

1 · d

s  2n κ

cs ∈([2n] κ ) κ(κ+1)/2

= (−1)

(E28)

17 G10,2 : P

1 · d

P

Q1κ

κ mod 2

=

P

i

·

d

Q0κ

P

P

Q0κ

iκ mod 2 = d

P

=0

(E29)

G1,1 : Q0κ′

Q0κ

= tr(Q0κ′ Q0κ ) = δκ,κ′

(E30)

= tr(Q1κ′ Q0κ ) = 0

(E31)

= tr(Q1κ′ Q1κ ) = δκ,κ′

(E32)

G1,2 : Q0κ

Q1κ′

G2,2 : Q1κ

Q1κ′

Appendix F: Bounds for the Second Moment of the FF-S-LCU

We derive the P second moment of the FF-S-LCU for arbitrary ρ0 , O. From Sec. III B, the single-layer average Λ := E[S1 ] = i (Db )i |Pi ⟩⟩⟨⟨Pi | gives E[m2O ] = ⟨⟨O, O|Λ l |ρ0 , ρ0 ⟩⟩.

(F1)

In terms of the moment operators Mπ1 , Mπ2 , Mπ3 of Sec. III B,
Λ = k · E[a4i ] · Mπ1 + k(k − 1) · E[a2i a2j ] · Mπ2 + Mπ3 .

(F2)

Each Mπm is an orthogonal projector, and the π2 , π3 ranges are spanned by that of π1 . This may be observed by inspection of the Gram matrix in Sec. IV 3 and Appendix E, and provides the following properties of the projectors, Mπ2m = Mπm (m = 1, 2, 3),

Mπ1 Mπ2 = Mπ2 ,

Mπ1 Mπ3 = Mπ3 ,

∥Mπ2 Mπ3 ∥ = 2/d.

(F3)

We define the projector composed of the first-order commutants as R := Mπ2 + Mπ3 , which commutes with the second-order commutant projector as Mπ1 R = RMπ1 . Expanding Λl , these properties collapse the expansion to terms composed of Mπ1 , R, and a term ∆, exponentially small in the number of qubits. Λl = (k · E[a4i ])l · Mπ1 + cl · R + ∆,

∥∆∥ ≤ 4(l − 1) · cl

4 l−2 1 1+ , d d

(F4)

l where cl = k · E[a4i ] + k(k − 1) · E[a2i a2j ] − (k · E[a4i ])l , and ∆=

l   X l r=2

r

(k E[a4i ])l−r (k(k − 1)E[a2i a2j ])r R r − R .

(F5)

18 For the ∥∆∥ bound, R is a projector up to exponentially small error O(1/d). R2 = R + W with W := Mπ2 Mπ3 + Mπ3 Mπ2 , ∥W ∥ ≤ 2∥Mπ2 Mπ3 ∥ = 4/d (adjoints have equal norm). Self-adjointness gives ∥R∥2 = ∥R2 ∥ ≤ ∥R∥ + ∥W ∥, so ∥R∥ ≤ 1 + 4/d. We may decompose the matrix term appearing in the expansion of ∆ as R r − R = (I + R + · · · + Rr−2 )W,

(F6)

the operator norm of which is bounded by (1 + d4 )r−2 . ∥R r − R∥ ≤ 4(r−1) d

(F7)

Using this bound in the sum in (F5) yields the bound of (F4). The remaining terms may be contracted with the boundary vectors, beginning with the second-order commutant projector, ⟨⟨O, O|Mπ1 |ρ0 , ρ0 ⟩⟩ =

X

2n h i

X j eκ (O) P eκ (ρ0 ) + Ceκ (O) Ceκ (ρ0 ) . tr O⊗2 Qκj tr ρ⊗2 P 0 Qκ =

(F8)

κ=0

j,κ

eκ = √σκ Pκ and Ceκ = √σκ Cκ are scaled and signed versions of the κ-purity and κ-coherence of Ref. [26], with Here P Dκ
Dκ σκ = (−1)⌊κ/2⌋ and Dκ = 2n κ . The terms involving the first-order moment projectors may similarly be calculated as 2 1 tr(O) tr(ρ ) + tr(P O) tr(P ρ ) , 0 0 d2   1 ⟨⟨O, O| Mπ3 |ρ0 , ρ0 ⟩⟩ = 2 tr(O2 ) tr(ρ20 ) + 2 tr(P O2 ) tr(P ρ20 ) + tr(P OP O) tr(P ρ0 P ρ0 ) . d

⟨⟨O, O| Mπ2 |ρ0 , ρ0 ⟩⟩ =

(F9) (F10)

Combining these gives the second moment. The error term for the second moment, ⟨⟨O, O|∆|ρ0 , ρ0 ⟩⟩, is at most ∥∆∥ tr(O2 ) tr(ρ20 ) (using ∥|X, X⟩⟩∥ = tr(X 2 )); with tr(ρ20 ) ≤ 1 and cl ≤ 1, and l ≪ d, this term is of order O(l 2−n tr(O2 )). Result 4: The FF-S-LCU second moment For any normalised initial state ρ0 and any Hermitian observable O, the second moment is E[m2O ] = (k E[a4i ])l

2n h X

eκ (O) P eκ (ρ0 ) + Ceκ (O) Ceκ (ρ0 ) P

i

κ=0

2 cl h + 2 tr(O) tr(ρ0 ) + tr(P O) tr(P ρ0 ) + tr(O2 ) tr(ρ20 ) d i + 2 tr(P O2 ) tr(P ρ20 ) + tr(P OP O) tr(P ρ0 P ρ0 )
+ O l 2−n tr(O2 ) ,

where cl = (k E[a4i ] + k(k − 1)E[a2i a2j ])l − (k E[a4i ])l .

1.

Concrete bounds on the second moment

The result provided above is general for the FF-S-LCU, and makes very few assumptions beyond the norm of ρ0 . With minimal additional restrictions, we may provide a concrete lower bound for the variance of the expectation value. Let oR , sR , GR be the restrictions of o, s, G to the eight first-order frame elements (the π2 , π3 block). Theorem 2. Let tr(O) = tr(P O) = 0 (O orthogonal to the first-order commutant, so that E[m] = 0 and Var[m] = E[m2 ]) and oR , sR ≥ 0. Then Var[m] ≥ (k E[a4i ])l ⟨⟨O, O|Mπ1 |ρ0 , ρ0 ⟩⟩ = (k E[a4i ])l Var[mU ],
where mU := tr U ρ0 U † O is the expectation value for a single Fermionic Gaussian unitary.

(F11)

19 Proof. By (F3), Λl = (k E[a4i ])l Mπ1 + term,

Pl

l 4 l−r (k(k − 1)E[a2i a2j ])r R r . r=1 r (k E[ai ])

Considering only the first-order

r−1 ⟨⟨O, O|R r |ρ0 , ρ0 ⟩⟩ = o⊤ sR . R GR

(F12)

Because each R r acts only within the first-order block, this walk involves only oR , sR and GR ; the second-order overlaps enter solely through the retained Mπ1 term. The entries of GR are traces of products of I and P (Appendix E), hence lie in {0, d1 , 1}, and oR , sR ≥ 0 by hypothesis; so every r ≥ 1 term is non-negative and drops, giving (F11). The factors decouple: (k E[a4i ])l ∼ k −3l carries the stacking, Var[mU ] the state and observable (closed form for all ρ0 , O [26]); any plateau is inherited from the twirl, not created by stacking. The hypotheses hold for even-parity ρ0 and parity-preserving O with tr(P O2 ) ≥ 0. Theorem 1. For a traceless quadratic observable O with tr(O2 ) = d, ρ0 = |0⟩⟨0|⊗n , n ≥ 3, with the variance taken over the joint distribution Haar(GFF )⊗kl × Dir(1, . . . , 1)⊗l , 1 Var[m] ≥ 2n − 1



24 (k + 1)(k + 2)(k + 3)

l

  1 . ∈Ω n k 3l

Proof. A traceless quadratic O is a Majorana bilinear: it commutes with P and carries all weight at degree κ = 2. So tr(O) = tr(P O) = 0 (the cc entries of oR vanish), and via P OP = O, tr(P OP O) = tr(O2 ) = d,

tr(P O2 ) = 0

(n ≥ 3),

(F13)

the latter since O2 has degree ≤ 4 < 2n. Hence oR = (0, 0, 0, 0, 1, 0, 0, 1) and sR = d1 (1, . . . , 1) (tr(ρm 0 ) = 1), both ≥ 0, so Theorem 2 applies. Only κ = 2 survives in Var[mU ], which for tr(O2 ) = d equals 1/(2n − 1) [26]. Substituting into (F11) gives Var[m] ≥ (k E[a4i ])l /(2n − 1); with k E[a4i ] = 24/((k+1)(k+2)(k+3)) (Appendix G) this is the stated bound. The case O = Z1 holds also at n = 2, as Z12 = I keeps tr(P Z12 ) = 0. Appendix G: Dirichlet Coefficients

For (a1 , . . . , ak ) ∼ Dir(1, . . . , 1), the uniform Dirichlet distribution over the (k − 1)-simplex, the following moments hold for i ̸= j: E[ai ] =

1 , k

(G1)

E[a2i ] =

2 , k(k + 1)

(G2)

24 , k(k + 1)(k + 2)(k + 3) 1 E[ai aj ] = , k(k + 1) 4 . E[a2i a2j ] = k(k + 1)(k + 2)(k + 3) E[a4i ] =

(G3) (G4) (G5)