Paper deep dive
Deep Sequence Modeling with Quantum Dynamics: Language as a Wave Function
Ahmed Nebli, Hadi Saadatdoorabi, Kevin Yam
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 89%
Last extracted: 7/20/2026, 12:48:14 PM
Summary
The paper introduces a sequence modeling framework where the latent state is a complex-valued wave function evolving under a learned, time-dependent Hamiltonian. This approach utilizes quantum interference for disambiguation, where conflicting interpretations cancel via destructive interference while compatible ones reinforce. The dynamics are strictly unitary, preserving state norm, and token probabilities are extracted via the Born rule. A key theoretical contribution is a separation theorem demonstrating that this complex unitary model solves certain disambiguation tasks with dimension N, whereas real-valued orthogonal models require dimension Ω(N²).
Entities (10)
Relation Signals (8)
Quantum Sequence Model → uses → Born Rule
confidence 95% · Token probabilities are extracted using the Born rule, a quadratic measurement operator that couples magnitudes and relative phases.
Quantum Sequence Model → uses → Hamiltonian
confidence 92% · Its evolution follows the time-dependent Schrödinger equation... where H(t) is a Hermitian operator.
Quantum Sequence Model → discretizes → Cayley Transform
confidence 90% · To implement this system, we discretize the evolution using the Cayley transform... produces an exactly unitary update
Quantum Sequence Model → comparedto → S4
confidence 85% · Our framework relates to SSMs in the following precise ways... This is a diagonal linear system identical in structure to the diagonalized form of S4
Quantum Sequence Model → comparedto → Mamba
confidence 85% · Mamba extended SSMs... This broke the time-invariance that enabled S4’s convolutional mode
Quantum Sequence Model → sharesconceptwith → Unitary Evolution RNN
confidence 85% · Our framework shares the use of unitary evolution and the Cayley transform with these models
Quantum Sequence Model → outperforms → LSTM
confidence 80% · Unlike standard recurrent architectures that rely on gating mechanisms... our framework utilizes quantum interference
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We introduce a sequence modeling framework in which the latent state is a complex-valued wave function evolving on a finite-dimensional Hilbert space under a learned, time-dependent Hamiltonian. Unlike standard recurrent architectures that rely on gating mechanisms to suppress competing hypotheses, our framework utilizes quantum interference: the Hamiltonian steers the phases of complex amplitudes so that conflicting interpretations cancel while compatible ones reinforce. The dynamics are strictly unitary, ensuring that the state norm is preserved exactly at every time step via a Cayley (Crank--Nicolson) discretization. Token probabilities are extracted using the Born rule, a quadratic measurement operator that couples magnitudes and relative phases. Our primary theoretical contribution is a separation theorem characterizing the representational advantage of this readout: we define a family of disambiguation tasks that a complex unitary model of dimension $N$ solves exactly, but which requires a state dimension of $\Omega(N^2)$ for any real-valued orthogonal model equipped with a standard affine-softmax readout. This quadratic gap arises because the Born rule implicitly lifts the $N$-dimensional state into the space of rank-one Hermitian matrices, accessing pairwise phase correlations that are inaccessible to linear projections. Finally, we derive a continuity equation for the latent probability mass, yielding conserved pairwise currents that serve as a built-in diagnostic for tracing information flow between dimensions.
Tags
Links
- Source: https://arxiv.org/abs/2602.22255v1
- Canonical: https://arxiv.org/abs/2602.22255v1
Trouble viewing inline? Open PDF directly →
Full Text
203,298 characters extracted from source content.
Expand or collapse full text
Deep Sequence Modeling with Quantum Dynamics: Language as a Wave Function Ahmed Nebli ahmed.nebli@cai-technology.ai cAI Technology GmbH Hadi Saadatdoorabi hadi.saadatdoorabi@cai-technology.ai cAI Technology GmbH Kevin Yam kevin.yam@coeo-group.ai cAI Technology GmbH Corresponding author Abstract We introduce a sequence modeling framework in which the latent state is a complex-valued wave function evolving on a finite-dimensional Hilbert space under a learned, time-dependent Hamiltonian. Unlike standard recurrent architectures that rely on gating mechanisms to suppress competing hypotheses, our framework utilizes quantum interference: the Hamiltonian steers the phases of complex amplitudes so that conflicting interpretations cancel while compatible ones reinforce. The dynamics are strictly unitary, ensuring that the state norm is preserved exactly at every time step via a Cayley (Crank–Nicolson) discretization. Token probabilities are extracted using the Born rule, a quadratic measurement operator that couples magnitudes and relative phases. Our primary theoretical contribution is a separation theorem characterizing the representational advantage of this readout: we define a family of disambiguation tasks that a complex unitary model of dimension N solves exactly, but which requires a state dimension of Ω(N2) (N^2) for any real-valued orthogonal model equipped with a standard affine-softmax readout. This quadratic gap arises because the Born rule implicitly lifts the N-dimensional state into the space of rank-one Hermitian matrices, accessing pairwise phase correlations that are inaccessible to linear projections. Finally, we derive a continuity equation for the latent probability mass, yielding conserved pairwise currents that serve as a built-in diagnostic for tracing information flow between dimensions. 1 Introduction Sequence modeling, the task of predicting the next element given a preceding context, underlies modern language models and drives applications from text generation to code synthesis. The dominant architectures for this task, including Transformers vaswani2017attention , recurrent neural networks hochreiter1997long , and state-space models gu2023mamba , differ substantially in their computational mechanisms but share a fundamental representational choice: the latent state at each time step is a vector of real numbers. In Transformers, context is stored in a cache of real-valued key-value pairs; in recurrent networks, it is a hidden vector updated multiplicatively; in state-space models, it is a vector governed by a discretized linear dynamical system. This paper examines the algebraic implications of this choice. In a real vector space, the superposition of two vectors is strictly additive. While real-valued neural networks can suppress incorrect hypotheses using learned nonlinear gating mechanisms (such as the sigmoid gates in LSTMs or the projection layers in Transformers), they lack the intrinsic geometric property of phase. In a complex vector space, superposition allows for interference: depending on the relative phase, adding two amplitudes can result in constructive reinforcement or destructive cancellation. Complex-valued representations provide a distinct mechanism for hypothesis interaction. A complex number carries both a magnitude and a phase. This structure allows competing interpretations of an ambiguous context to interact algebraically rather than requiring a dedicated gating module to suppress one against the other. Consider the prefix “The bank was …” where a model must maintain latent probability mass for both the ‘financial institution’ and ‘river edge’ interpretations. If the subsequent token is “steep,” it acts as a filter that should reinforce the river interpretation and suppress the financial one. In the framework proposed here, this suppression arises from interference: the evolution of the state rotates the phases of the latent components such that the amplitudes associated with the financial interpretation interfere destructively with the new input, while the river interpretation interferes constructively. This architectural choice is inspired by, though distinct from, the field of quantum cognition. Busemeyer and Bruza busemeyer2012quantum and Pothos and Busemeyer pothos2013quantum have documented that human probability judgments often violate classical probability axioms in ways formally described by quantum probability theory. While our goal is to model the statistical distribution of text rather than human cognitive processes, these findings suggest that the mathematical structure of complex Hilbert spaces may offer a parsimonious inductive bias for sequence disambiguation. There is also growing interest in formulating neural network dynamics in continuous time. Chen et al. chen2018neural formalized residual networks as discretizations of ordinary differential equations (Neural ODEs). Greydanus et al. greydanus2019hamiltonian introduced Hamiltonian Neural Networks to enforce conservation laws in physical simulations. haber2017stable demonstrated that constraining dynamics to be norm-preserving (orthogonal or unitary) provides strong stability guarantees, mitigating the vanishing and exploding gradient problems. Modern state-space models gu2023mamba ; gu2021efficiently leverage continuous-time linear systems for efficient sequence processing. However, existing approaches typically operate in real-valued spaces or rely on unconstrained dynamics that may lack unitary stability. None combine complex-valued states with structured Hamiltonian evolution and a quadratic, measurement-based output rule. This paper introduces a framework integrating these concepts. The latent state is a unit-norm vector in a finite-dimensional complex Hilbert space. Its evolution follows the time-dependent Schrödinger equation idt|ψ(t)⟩=H(t)|ψ(t)⟩i ddt ψ(t)=H(t) ψ(t), where H(t)H(t) is a Hermitian operator. The Hermiticity of H(t)H(t) guarantees that the time-evolution operator is strictly unitary, ensuring the state vector rotates on the unit sphere without changing its norm. We decompose H(t)=H0+Hint(t)H(t)=H_0+H_int(t), separating the dynamics into a learnable diagonal term H0H_0 that sets baseline oscillation frequencies, and an input-dependent interaction Hamiltonian Hint(t)H_int(t) generated by a neural network. This decomposition allows us to move into the interaction picture, a coordinate change that factors out the free oscillations, leaving the numerical integrator to resolve only the input-driven dynamics. To implement this system, we discretize the evolution using the Cayley transform, which corresponds to the Crank–Nicolson (implicit midpoint) scheme of numerical analysis hairer2006geometric ; crank1947practical . Unlike standard explicit integrators (e.g., Runge-Kutta), which destroy the unitary property and introduce norm drift, the Cayley transform produces an exactly unitary update for any step size. This ensures the model preserves the state norm regardless of sequence length, although we note that discretization error in the phase trajectory (integration accuracy) remains bounded by the step size. Token prediction is performed via the Born rule: the probability of token k is the squared magnitude of the inner product between the state and a learned measurement vector. This introduces a quadratic nonlinearity at the readout layer, differing fundamentally from the linear-projection-plus-softmax used in standard models. As we show, this quadratic structure allows the model to access pairwise phase correlations that are inaccessible to a linear readout. The central theoretical contribution of this paper is a separation theorem characterizing the representational capacity of this architecture. We construct a family of disambiguation tasks in which the correct output depends on the phase relationship between two context tokens. We prove that a complex-valued unitary model of dimension N can represent these tasks exactly, whereas any real-valued model restricted to orthogonal dynamics and an affine-softmax readout requires a state dimension of Ω(N2) (N^2). The gap arises because the Born rule implicitly lifts the state into the space of Hermitian matrices, accessing O(N2)O(N^2) degrees of freedom (including phase cross-terms) from an N-dimensional complex vector. Real-valued models with linear readouts must explicitly encode these pairwise interactions in their state dimensions. The lower bound (Theorem 5.10) is proved conditionally, assuming the target log-probability matrix has full row rank (a genericity-type condition on the task parameters, discussed in Section 5), and is established for the state-independent model class defined in Section 5; the extension to the full state-dependent architecture of Section 3 is left for future work. Finally, we derive a continuity equation for the probability mass in the latent dimensions. Because the dynamics are unitary, the change in occupation probability of any dimension is exactly accounted for by antisymmetric probability currents flowing between dimensions. We propose these currents as a built-in diagnostic tool for tracing the internal redistribution of semantic content. Quantum Formalism Architectural Role ML Interpretation Hilbert space ℂNC^N Complex latent space with magnitude and phase per coordinate Hidden state space Wave function |ψ(t)⟩ ψ(t) Unit-norm complex vector; superposition of interpretations Recurrent hidden state Hermitian H(t)=H†H(t)=H Constrains dynamics to be norm-preserving by construction Learned transition rule Schrödinger equation H†=H⇒H =H unitary ⇒ ‖ψ‖≡1\|ψ\|≡1 Recurrence equation Cayley transform Structure-preserving discretization; unitarity exact for any step size Discrete state update Born rule |⟨mk|ψ⟩|2|\! m_kψ\!|^2 Quadratic readout: N2N^2 features from N complex dims Token probabilities Interference (phase) Competing interpretations cancel or reinforce via relative phase Disambiguation mechanism Prob. current Jj←kJ_j← k Antisymmetric, conserved pairwise flow driven by each token Built-in interpretability Table 1: Correspondence between quantum-mechanical formalism and the sequence model. Each structure plays a specific functional role: Hermiticity guarantees unitarity, unitarity guarantees norm conservation, and norm conservation ensures the Born rule produces valid probabilities. The architecture uses this algebraic chain as an engineering constraint, not as a physical simulation. 2 Related Work Our framework combines three ideas: (i) complex-valued latent states, (i) Hamiltonian continuous-time dynamics, and (i) measurement-based decoding, into a single architecture for next-token prediction. Each of these ideas has antecedents in the literature, but no prior work combines all three, and several structural limitations of existing approaches motivate the specific design choices we make. We organize this section around the seven research programs most directly relevant to our proposal. 2.1 Unitary and Orthogonal Recurrent Networks The observation that recurrent neural networks suffer from vanishing and exploding gradients during backpropagation through time motivated a line of work that constrains the hidden-to-hidden transition matrix to be orthogonal or unitary. Arjovsky et al. arjovsky2016unitary introduced the Unitary Evolution RNN (uRNN), parameterizing the transition matrix as a product of diagonal phase matrices, Householder reflections, discrete Fourier transforms, and permutations. This composition is always unitary, which guarantees that the gradient norm neither grows nor shrinks across time steps. The authors demonstrated the benefit on synthetic benchmarks requiring long-term memory, such as the Copying Task and the Adding Problem, where the uRNN maintained stable gradients over hundreds of time steps where standard LSTMs degraded hochreiter1997long . However, the fixed factorization used by Arjovsky et al. covers only a strict subset of the unitary group U(N)U(N). Wisdom et al. wisdom2016full addressed this limitation by parameterizing the full unitary group via a gradient-based optimization on the Stiefel manifold absil2008optimization , using the Cayley transform to map skew-Hermitian matrices to unitary matrices. Their Full-Capacity Unitary RNN achieved lower error on the Adding Problem and the Copying Task at sequence lengths up to T=2000T=2000. Mhammedi et al. mhammedi2017efficient proposed parameterizing orthogonal matrices via a product of Householder reflections, obtaining a strictly expressive family with full coverage of O(N)O(N) at cost O(N2)O(N^2) per step while avoiding the manifold retraction required by gradient-based Stiefel methods. lezcano2019cheap subsequently offered a unified extrinsic parameterization of orthogonal and unitary matrices through the matrix exponential of skew-symmetric or skew-Hermitian matrices, equipped with an efficient Cayley-based retraction. Vorontsov et al. vorontsov2017orthogonality further studied orthogonality constraints in recurrent networks and their effect on gradient propagation. Our framework shares the use of unitary evolution and the Cayley transform with these models, but differs in three structural respects that are consequential for language modeling. First, both uRNNs and Full-Capacity Unitary RNNs typically define the recurrence as a discrete algebraic update ht+1=Wht+Vxth_t+1=Wh_t+Vx_t, where W is a unitary matrix and VxtVx_t is an additive input injection. In this formulation, the input token shifts the state but does not alter the dynamics that govern state evolution. In our framework, the input token at time t parameterizes the interaction Hamiltonian Hint(t)H_int(t), which is the generator of the evolution itself. This means the input controls the rotation axis and angular velocity of the state trajectory on the unit sphere, allowing the input to restructure the interference pattern among active interpretations rather than merely displacing the state vector. Second, prior unitary RNNs typically decode the hidden state through a standard real-valued affine projection followed by softmax. This discards the phase information in the complex hidden state, collapsing it to magnitudes before the output layer. Our Born-rule decoding computes p(k)=|⟨mk|ψ(t)⟩|2p(k)=| m_k|ψ(t) |^2, which is a quadratic function of the complex amplitudes. As we discuss in Section 3.5, this quadratic structure allows the output probabilities to depend on cross-term interference between components of |ψ(t)⟩ ψ(t), a mechanism that a linear projection followed by softmax cannot replicate without increasing the latent dimension. Third, while prior work focuses on gradient stability, it does not establish a formal representational separation between complex unitary and real orthogonal models. Our separation theorem (Section 5) constructs a family of disambiguation tasks and proves that a complex unitary model of dimension N computes these tasks exactly while any real orthogonal model with an affine-softmax readout requires dimension Ω(N2) (N^2). This identifies a specific readout-dependent bottleneck in standard real-valued architectures. 2.2 State-Space Models Structured State Space Models (SSMs) form the continuous time framework closest to our own. The S4 model introduced by Gu et al. gu2021efficiently defines the latent dynamics as a linear system x′(t)=Ax(t)+Bu(t)x (t)=Ax(t)+Bu(t), y(t)=Cx(t)+Du(t)y(t)=Cx(t)+Du(t), where A is initialized with the HiPPO matrix, a specific structure designed to optimally compress continuous signals. S4 discretizes this system and computes the resulting convolution via FFT, achieving strong results on the Long Range Arena benchmark. Mamba gu2023mamba extended SSMs by making the matrices B, C, and the discretization step size Δ functions of the input token, introducing input-dependent selection. This broke the time-invariance that enabled S4’s convolutional mode but allowed the model to selectively propagate or forget information, achieving language modeling perplexity competitive with Transformers. Smith et al. smith2023simplified subsequently simplified and extended the S4 lineage with the S5 model, showing further efficiency gains through parallel scans. Poli et al. poli2023hyena introduced the Hyena operator, a long-convolution approach that replaces attention with implicit parameterized convolutions, achieving competitive performance on language modeling while scaling sub-quadratically with sequence length; the analysis of Hyena underscores that the choice of state-update operator—linear convolution versus Hamiltonian flow—has direct consequences for the model’s expressive power and inductive bias. Our framework relates to SSMs in the following precise ways. When the interaction Hamiltonian is set to zero (Hint(t)=0H_int(t)=0) and the free Hamiltonian H0H_0 is diagonal, the Schrödinger equation idt|ψ⟩=H0|ψ⟩i ddt ψ=H_0 ψ reduces to N decoupled oscillators dtcj=−iωjcj ddtc_j=-i _jc_j. This is a diagonal linear system identical in structure to the diagonalized form of S4, except that our state variables are natively complex. In our framework, the complex structure is semantically load-bearing: phases encode relational information between interpretations, and the Born rule converts phase relationships into output probabilities through interference. Critically, the linearity of SSMs imposes a constraint: the state update at each step is a linear map applied to the current state. Our Hamiltonian framework introduces nonlinearity through the interaction term Hint(t)H_int(t): because Hint(t)H_int(t) is generated by a neural network applied to the token embedding and the current state (as detailed in Section 3.2), the map from input sequence to state trajectory is nonlinear, even though the instantaneous evolution under any fixed H(t)H(t) is linear. 2.3 Complex-Valued Neural Networks Trabelsi et al. trabelsi2018deep introduced a systematic framework for deep complex-valued networks, defining complex convolutions and batch normalization. They demonstrated improvements over real-valued baselines on tasks like MusicNet, attributing the gains to the ability of complex convolutions to learn phase-sensitive filters. Hirose and Yoshida hirose2012generalization provide an earlier theoretical treatment of complex-valued neural networks and their generalization properties. Virtue et al. virtue2017better showed that complex-valued networks can outperform their real-valued counterparts in MRI fingerprinting tasks, where the underlying physics produces complex-valued measurement signals and phase structure carries diagnostic information; this domain-specific result provides additional evidence that complex representations are most beneficial when the problem’s native structure is phase-bearing. However, the empirical literature is not uniformly positive; other works have found that real-valued models can match complex-valued performance if the parameter count is equalized by increasing the real network’s width. This suggests that complex values are not a panacea, but rather a specific inductive bias suitable for specific problems. Our separation theorem provides a formalization of when this bias is useful in sequence modeling. It identifies a specific mechanism—destructive interference between complex amplitudes accessed via a quadratic readout—and proves that it enables a dimension reduction by a quadratic factor relative to real-valued models with linear readouts. The proof constructs disambiguation functions in which the correct output depends on the relative phase between two context tokens. A complex unitary model encodes the interaction in the relative phases of the state vector, while a real orthogonal model must represent pairwise interactions as separate state dimensions. 2.4 Neural ODEs and Hamiltonian Networks Chen et al. chen2018neural showed that residual networks can be interpreted as discretizations of an ODE dhdt=f(h,θ,t) dhdt=f(h,θ,t). While powerful, unconstrained Neural ODEs face topological limitations; Dupont et al. dupont2019augmented showed that Neural ODEs in ℝnR^n cannot represent functions that require trajectories to cross. They resolved this by augmenting the state space. Kidger et al. kidger2021neural further developed neural controlled differential equations for sequence modeling, demonstrating that continuous-time recurrences can match or exceed discrete recurrent models on irregular time-series. Our framework resolves trajectory crossings through phase: two states with identical magnitudes |cj||c_j| but different phases occupy different points on the complex unit sphere. Phase provides an extra degree of freedom per complex coordinate, and the Hamiltonian dynamics naturally exploit this by rotating phases at rates determined by the learned frequencies. Greydanus et al. greydanus2019hamiltonian introduced Hamiltonian Neural Networks (HNNs) to enforce energy conservation in physical simulations. However, HNNs operate in real-valued phase spaces and conserve a scalar energy function. In contrast, our framework conserves the L2L^2 norm of a complex state vector. This conservation of total probability yields vector-valued probability currents (derived in Section 4) rather than a scalar energy value, serving as a diagnostic for information flow rather than physical energy. Stoudenmire and Schwab stoudenmire2016supervised demonstrated that tensor networks, which share the unitary-constraint philosophy of our approach, can be competitive classifiers, further motivating norm-preserving parameterizations in machine learning. 2.5 Liquid Neural Networks and Continuous-Time Recurrences Liquid Neural Networks (LNNs), introduced by Hasani et al. hasani2021liquid , define recurrent dynamics through a system of ordinary differential equations whose time constants are modulated jointly by the input signal and the current state. The core governing equation for the hidden state h(t)∈ℝdh(t) ^d is h˙(t)=−h(t)τ(h(t),x(t),θ)+fθ(h(t),x(t)), h(t)=- h(t)τ(h(t),x(t),θ)+f_θ(h(t),x(t)), (1) where τ>0τ>0 is a state- and input-dependent time constant determining how quickly each neuron decays toward the input-driven target, and fθf_θ is a learned nonlinearity. The resulting dynamics are “liquid” in the sense that each neuron’s effective temporal receptive field adapts to the current input, compressing or dilating its timescale in response to input statistics. Lechner et al. lechner2020neural demonstrated that wiring topologies inspired by the nematode C. elegans connectome produce remarkably compact and interpretable controllers for continuous control tasks, achieving state-of-the-art performance with far fewer parameters than generic Neural ODEs. The Closed-form Continuous-time (CfC) variant hasani2022closed derives an approximate closed-form solution to these ODEs, eliminating the need for an online numerical integrator during inference and recovering the computational efficiency of discrete-step recurrence while retaining the expressive inductive bias of continuous-time dynamics. LNNs share with our framework the fundamental architectural choice of grounding the state update in a differential equation rather than in a hand-designed recurrence rule, and both frameworks allow the input signal to modulate the transition, rather than contributing only an additive displacement. These structural similarities make LNNs the most direct continuous-time comparator to the architecture proposed here. Several differences distinguish the two frameworks, however, and they are consequential for the representational and stability properties relevant to language modeling. First, LNN dynamics operate in real-valued state spaces with unconstrained weight matrices; no architectural invariant controls the spectral properties of the Jacobian of fθf_θ, and stability is ensured through the sign of the decay term −h/τ-h/τ rather than through an exact structural constraint. Our framework uses a Hermitian Hamiltonian, which structurally guarantees that the continuous-time flow is norm-preserving (unitary) and that, via the Crank–Nicolson (Cayley) discretization (Section 3.4), this property holds exactly at every discrete step. Second, the CfC approximation replaces the ODE solution with an exponential-decay closed form; the quality of this approximation degrades when the nonlinear term in fθf_θ is large relative to the linear decay, creating a tradeoff between efficiency and representational fidelity that does not arise in our framework, where the Cayley update is exact with respect to norm preservation regardless of Hamiltonian magnitude. Third, LNNs decode through an affine-softmax layer, discarding any phase-like structure that the ODE trajectory may have produced; our Born-rule decoding is quadratic in the state, providing access to O(N2)O(N^2) pairwise features (Lemma 5.7) that the affine-softmax readout cannot access without a corresponding linear-dimension increase (Lemma 5.9). Taken together, the two frameworks represent complementary continuous-time strategies: LNNs prioritize adaptive timescales and biological plausibility in real-valued spaces, while our framework prioritizes algebraic structure preservation—norm conservation and phase-sensitive readout—in complex-valued spaces. Whether the additional structural constraints of our framework translate to improved performance on natural language tasks remains an empirical question addressed by the experimental protocols of Section 6.2. 2.6 Quantum Cognition and Quantum-Inspired NLP The hypothesis that human judgment exhibits interference effects has been developed by Busemeyer and Bruza busemeyer2012quantum . Empirical evidence, such as order effects and the conjunction fallacy, suggests that human probability judgments sometimes violate classical axioms in ways described by quantum probability. For example, the probability of “A and B” can be rated higher than the probability of “A” alone, which is naturally explained by constructive interference in a projector-based measurement model. We cite this literature not to claim that language models must simulate human cognitive processes, but to motivate the architecture as a potentially useful inductive bias. If human language production leaves traces of these interference-like structures, a model equipped with complex amplitudes and Born-rule decoding may represent them more efficiently than a classical probability model. The Born-rule measurement postulate has begun to attract direct empirical attention in NLP: recent work on text classification with Born’s rule guidotti2022text provides evidence that squared-amplitude measurement operators can serve as effective output mechanisms for natural language tasks, contextualizing our architectural choice within an emerging body of quantum-inspired NLP research. Similarly, quantum-inspired feature maps in reproducing kernel Hilbert spaces schuld2019quantum provide theoretical grounding for using complex inner products as similarity measures in high-dimensional feature spaces. Existing quantum cognition and quantum-inspired NLP models are largely either descriptive (fitting data with hand-selected operators) or restricted to small-scale demonstrations. Our framework is a trainable realization of this mathematical structure, scaling it to vocabulary sizes of tens of thousands and learning the operators from data. 2.7 Quantum Machine Learning and Quantum Neural Networks The term quantum machine learning (QML) encompasses two related but distinct paradigms that are important to distinguish from our approach. The first, and currently dominant, paradigm investigates how machine learning computations can be accelerated or qualitatively enriched by running on actual quantum processors, so-called near-term noisy intermediate-scale quantum (NISQ) devices preskill2018quantum . The second paradigm develops quantum-informed mathematical structures as inductive biases for models that run entirely on classical hardware, of which our framework is an instance. We discuss both, as each informs the design and interpretation of the present work. Parameterized quantum circuits and quantum neural networks. In the NISQ paradigm, quantum neural networks (QNNs) are typically implemented as parameterized quantum circuits (PQCs) cerezo2021variational ; mitarai2018quantum : a sequence of parameterized unitary gates drawn from a hardware-native gate set is applied to qubits initialized in a standard state, and a measurement in the computational basis extracts classical outputs. The measurement step is governed, as in our framework, by Born’s rule: a projective measurement on a state |ψ⟩ ψ produces outcome k with probability |⟨k|ψ⟩|2| k|ψ |^2. This shared measurement postulate makes the mathematics of PQC-based classifiers directly analogous to our Born-rule decoding layer. The key distinction is the substrate and the dimension: PQCs operate on quantum amplitude vectors that physically exist in a 2n2^n-dimensional Hilbert space for n qubits (with quantum interference arising from genuine quantum superposition), whereas our framework maintains a classically stored, explicitly computed state vector in ℂNC^N. The expressivity and trainability of QNNs are active research areas. biamonte2017quantum provided an influential early survey situating quantum machine learning within the broader landscape of quantum algorithms. Havlíček et al. havlicek2019supervised demonstrated that quantum feature maps can produce kernel methods with potential classification advantages on specific structured datasets. Abbas et al. abbas2021power characterized the effective dimension and trainability of QNNs, showing that they can exhibit high capacity per parameter in certain regimes. However, the “barren plateau” phenomenon mcclean2018barren poses a significant trainability obstacle analogous to (though more severe than) the vanishing-gradient problem in classical deep networks; recent work by Cerezo et al. cerezo2021cost has analyzed the conditions under which local cost functions resist barren plateaus in shallow circuits. Relation to the present framework. Our framework relates to the QML literature in three ways. First, we import the Born-rule measurement postulate and the Hermitian-Hamiltonian evolution law as algebraic constraints on a classical sequence model, not as descriptions of any physical quantum process. The representational advantages we establish (O(N2)O(N^2) effective features from a Born-rule readout (Lemma 5.7), exact norm-conservation from Hermitian dynamics) are purely algebraic consequences that apply on classical hardware. Second, the gradient stability analysis of Section 3.6 is structurally analogous to the barren-plateau analysis of PQCs: in both cases the central question is whether useful gradient signals propagate through a composition of unitary transformations. Our Crank–Nicolson discretization guarantees that the gradient of the loss with respect to the recurrent state propagates without change in norm along the state pathway, providing an architectural resolution to the state-gradient form of this problem, although gradients through the parameter pathway (inside gθg_θ) are not covered by this guarantee, as discussed in Section 3.6. Third, the data-encoding and measurement structure of our Born-rule layer resembles the data-reuploading classifiers of Pérez-Salinas et al. perezsalinas2020data , in which classical data is injected at multiple circuit layers to improve expressivity; in our model, the analogous role is played by the state-dependent Hamiltonian Hint(t)H_int(t), which re-encodes the current token at every step in a manner whose depth is effectively determined by the sequence length. Quantum reservoir computing and quantum-inspired classical algorithms. A third connection arises through quantum reservoir computing fujii2017harnessing ; ghosh2019quantum , which uses the dynamics of a (fixed, untrained) quantum system as a nonlinear feature map whose output is read by a simple trained linear layer. Quantum reservoir computers are the quantum analogue of echo state networks jaeger2004harnessing ; only the output weights are trained, making the approach computationally lightweight but limiting the model’s capacity to that of a kernel method applied to the reservoir’s fixed representational space. In contrast, our framework trains the full dynamics, the Hamiltonian, measurement, and initial state, end-to-end. The broader quantum-inspired classical algorithms literature tang2019quantum explores problems where the mathematical structure of quantum algorithms (amplitude amplification, phase estimation) can be simulated efficiently on classical hardware to yield formal algorithmic speedups; our work belongs to this tradition in motivation but targets representational rather than computational advantages. We emphasize that the present architecture is entirely classical in its implementation. The state |ψ(t)⟩ ψ(t) is stored as a vector of floating-point complex numbers; the Cayley update is a deterministic linear algebra operation; the Born-rule output is the squared modulus of an inner product, all computable on commodity hardware without access to quantum processors. The physical terminology (Hamiltonian, wave function, Born rule) names mathematical structures whose algebraic properties provide the guarantees we prove. This terminological choice is conventional in the quantum-inspired ML literature and should not be read as claiming a computational quantum advantage. 2.8 Transformers and Attention Mechanisms Transformers vaswani2017attention dominate large-scale language modeling. The attention mechanism provides direct access to context, bypassing the recurrent bottleneck. Our framework differs in its computational structure. First, the Transformer’s memory (KV cache) grows as O(TdL)O(TdL) with sequence length T, hidden dimension d, and number of layers L. More importantly, the standard scaled dot-product self-attention operation requires O(T2d)O(T^2d) time and O(T2+Td)O(T^2+Td) memory per layer due to the all-pairs attention matrix, a quadratic dependence on sequence length that becomes prohibitive for long documents. Efficient attention variants, including linear attention katharopoulos2020transformers , sparse attention such as BigBird zaheer2020big , and hardware-optimized exact attention dao2022flashattention , reduce this to O(Td)O(Td) or O(TlogTd)O(T T\,d), at the cost of approximating or restructuring the full attention pattern. Our recurrent state is fixed-size O(N)O(N), independent of sequence length. Second, the softmax attention mechanism computes a convex combination of value vectors. While the Multi-Layer Perceptron (MLP) layers following attention can implement subtraction or cancellation, the attention operation itself is an averaging process. In our framework, cancellation is an intrinsic property of the state accumulation: two latent components can cancel through destructive interference during the summation of state updates. This suggests a different mechanism for suppressing irrelevant information: rather than driving an attention weight to zero via a saturating softmax, the model can align phases to produce destructive interference. 3 The Quantum Sequence Model Single time step — internal mechanismToken xtx_tEmbed(xt)Embed(x_t)look-upNeural net gθg_θΦ(t)∈ℂN×r (t) ^N× rδ(t)∈ℝNδ(t) ^NΦΦ† diag(δ)diag(δ)Hermitianreal diagonalHint(t)=ΦΦ†+diag(δ)H_int(t)= +diag(δ)Hint=Hint†H_int=H_int H0=diag(λj)H_0=diag( _j)learnedH(t)=H0+Hint(t)H(t)=H_0+H_int(t)H=H†H=H Hint,I(t)=eiH0tHint(t)e−iH0tH_int,I(t)=e^iH_0t\,H_int(t)\,e^-iH_0tremove free oscillationsW(t)=(I+iΔt2Hint,I)−1(I−iΔt2Hint,I)W(t)= (I+i t2H_int,I )^\!-1 (I-i t2H_int,I )W†W=IW W=I (exact)|ψI(t+1)⟩=W(t)|ψI(t)⟩ _I(t+1)=W(t)\, _I(t)‖ψ‖≡1\|ψ\|≡ 1p(k|ψ)=|⟨mk|ψ⟩|2=tr(Mkρ)p(k|ψ)=|\! m_kψ\!|^2=tr(M_kρ)quadratic: N2N^2 features|mk⟩k=1V\|m_k \_k=1^Vlearned, orthonormalℒ=−∑tlogp(yt∣ψ(t))L=-\! _t p(y_t ψ(t))state feedback ψI(t) _I(t)Guarantee chain:ΦΦ†=(ΦΦ†)† =( ) (Hermitian) ⟶ H=H†H=H ⟶ W†W=IW W=I (Cayley) ⟶ ‖ψ‖≡1\|ψ\|≡1 ⟶ ∑kp(k)=1 _kp(k)=1 (valid distribution) Figure 1: Detailed architecture of a single time step in the quantum sequence model. The neural network gθg_θ (the sole unconstrained learned component, bold outline) receives the token embedding and the current interaction-picture state, and outputs the complex matrix Φ(t) (t) and real vector δ(t)δ(t). The outer product ΦΦ† is Hermitian by construction; adding the learned diagonal H0H_0 yields the full Hamiltonian H(t)H(t). The interaction picture removes the known free oscillations, and the Cayley transform discretizes the remaining evolution into an exactly unitary update W(t)W(t). Unitarity preserves ‖ψ‖=1\|ψ\|=1 at every step, which ensures the Born rule produces a valid probability distribution over the vocabulary. Left-side annotations mark learned components; right-side annotations trace the algebraic guarantees. The dashed line indicates recurrent state feedback, which makes the overall dynamics nonlinear despite each individual step being a linear (unitary) map. The Introduction argued that complex-valued representations support interference, a mechanism by which competing interpretations can suppress one another through phase opposition rather than through explicit gating. Section 2 surveyed four families of models that each capture part of this idea: unitary RNNs use complex-valued states but discard phase at the output layer; state-space models use continuous-time dynamics but operate in real-valued spaces; complex-valued networks employ complex parameters but impose no structural constraints on the dynamics; and Hamiltonian neural networks preserve conservation laws but operate in real phase spaces designed for physical simulation. No existing architecture combines complex-valued states, Hamiltonian evolution, and phase-sensitive decoding into a single trainable sequence model. This section constructs such an architecture in six steps. We first define the state space (Section 3.1), explaining what structure a complex unit-norm vector provides and why that structure is relevant to the interference mechanism. We then specify the evolution law (Section 3.2), decomposing it into a free component that establishes a spectrum of oscillation timescales and an input-driven component that couples latent dimensions in response to each observed token. A change of coordinates (Section 3.3) factors out the known free oscillations, isolating the input-driven dynamics in a form suited to numerical integration. A structure-preserving discretization (Section 3.4) converts the continuous evolution into a discrete update that maintains unitarity regardless of step size. A measurement rule (Section 3.5) extracts token probabilities from the complex state through a mechanism whose algebraic properties differ from the standard softmax readout. Finally, the initialization and training objective (Section 3.6) complete the model specification. Table 2: Summary of all mathematical notation used throughout the paper. All complex vectors are column vectors; ∥⋅∥\|·\| denotes the Euclidean (L2L^2) norm; (⋅)†(·) denotes the conjugate transpose. Symbol Space Definition Scalars and index sets N ℤ>0Z_>0 Latent (complex Hilbert-space) dimension V ℤ>0Z_>0 Vocabulary size T ℤ>0Z_>0 Sequence length d ℤ>0Z_>0 Token embedding dimension r ℤ>0Z_>0 Rank of Φ(t) (t); controls the interaction Hamiltonian’s expressivity-efficiency tradeoff (r≪Nr N) Δt t ℝ>0R_>0 Discretization step size (typically Δt=1 t=1: one step per token) λj _j ℝR j-th learnable oscillation frequency (eigenvalue of H0H_0) State vectors |ψ(t)⟩ ψ(t) ℂNC^N Schrödinger-picture latent state; ‖ψ(t)‖=1\|ψ(t)\|=1 |ψI(t)⟩ _I(t) ℂNC^N Interaction-picture latent state; |ψ(t)⟩=e−iH0t|ψI(t)⟩ ψ(t)=e^-iH_0t _I(t) cj(t)c_j(t) ℂC j-th amplitude: cj(t)=⟨j|ψ(t)⟩c_j(t)= jψ(t) pj(t)p_j(t) [0,1][0,1] Occupation probability of dimension j: pj(t)=|cj(t)|2p_j(t)=|c_j(t)|^2 Hamiltonians and evolution H(t)H(t) ℂN×NC^N× N Full Hamiltonian: H(t)=H0+Hint(t)H(t)=H_0+H_int(t); Hermitian H0H_0 ℂN×NC^N× N Free Hamiltonian: diag(λ0,…,λN−1)diag( _0,…, _N-1) Hint(t)H_int(t) ℂN×NC^N× N Interaction (input-driven) Hamiltonian; Hermitian Hint,I(t)H_int,I(t) ℂN×NC^N× N Interaction Hamiltonian in the interaction picture: eiH0tHint(t)e−iH0te^iH_0tH_int(t)e^-iH_0t (t2,t1)U(t_2,t_1) U(N)U(N) Unitary time-evolution operator mapping |ψ(t1)⟩ ψ(t_1) to |ψ(t2)⟩ ψ(t_2) W(t)W(t) U(N)U(N) Cayley (Crank–Nicolson) discrete update operator Interaction Hamiltonian parameterization gθg_θ — Neural network mapping (Embed(xt),[Re(cI(t)),Im(cI(t))])→(Φ(t),δ(t))(Embed(x_t),\,[Re(c_I(t)),Im(c_I(t))])→( (t),δ(t)) Φ(t) (t) ℂN×rC^N× r Low-rank factor of Hint(t)H_int(t); output of gθg_θ δ(t)δ(t) ℝNR^N Diagonal-shift vector; output of gθg_θ Φ~(t) (t) ℂN×rC^N× r Interaction-picture version: [Φ~]ja=eiλjtΦja(t)[ ]_ja=e^i _jt _ja(t) Embed(x)Embed(x) ℝdR^d Learned token embedding for vocabulary item x Measurement and output |mk⟩ m_k ℂNC^N Measurement vector for vocabulary item k M ℂN×VC^N× V Measurement matrix; columns are |mk⟩ m_k; satisfies MM†=INMM =I_N MkM_k ℂN×NC^N× N Rank-one projector |mk⟩⟨mk| m_k m_k ρ ℂN×NC^N× N Density matrix |ψ⟩⟨ψ| ψ ψ; rank-one Hermitian p(k∣ψ(t))p(k ψ(t)) [0,1][0,1] Born-rule token probability: |⟨mk|ψ(t)⟩|2| m_kψ(t)|^2 Probability currents Jj←k(t)J_j← k(t) ℝR Continuous-time probability current from dimension k to dimension j Jj←kmid(t)J_j← k^mid(t) ℝR Midpoint (exact discrete) probability current Separation theorem ND_N — Disambiguation task family parameterized by N CUSM — Complex Unitary Sequence Model (Definition 5.1) ROSM — Real Orthogonal Sequence Model (Definition 5.2) ρij _ij ℂN×NC^N× N Task density matrix Wj|ψi⟩⟨ψi|Wj†W_j _i _iW_j 3.1 State Space The interference mechanism described in the Introduction requires that the latent representation carry both magnitudes and phases: two contributions to a prediction cancel only if their phases oppose one another. To realize this, we define the latent state at discrete time step t as a unit-norm vector in a finite-dimensional complex Hilbert space: |ψ(t)⟩∈ℂN,⟨ψ(t)|ψ(t)⟩=1, ψ(t) ^N, ψ(t)|ψ(t) =1, (2) where N is a hyperparameter setting the model’s capacity and ⟨⋅|⋅⟩ ·|· is the standard Hermitian inner product. Expanding in the computational basis |0⟩,|1⟩,…,|N−1⟩\ 0, 1,…, N-1\, the state takes the form |ψ(t)⟩=∑j=0N−1cj(t)|j⟩, ψ(t)= _j=0^N-1c_j(t) j, (3) where each coefficient cj(t)∈ℂc_j(t) is called an amplitude. The unit-norm constraint requires ∑j|cj(t)|2=1 _j|c_j(t)|^2=1, so the squared magnitudes form a probability distribution over the N latent dimensions at every time step. Each amplitude carries two degrees of freedom. In polar form, cj(t)=rj(t)eiθj(t)c_j(t)=r_j(t)\,e^i _j(t): the magnitude rj(t)=|cj(t)|r_j(t)=|c_j(t)| determines how much weight the j-th dimension carries, while the phase θj(t)=argcj(t) _j(t)= c_j(t) encodes relational information that has no counterpart in real-valued representations under a linear readout. To see concretely why phases matter, consider two states that share the same magnitude distribution (e.g., say |c1|=|c2|=1/2|c_1|=|c_2|=1/ 2 with all other amplitudes zero) but differ in relative phase: in one state c2=c1c_2=c_1, and in the other c2=−c1c_2=-c_1. Under a linear readout w⊤hw h applied to the squared magnitudes, these two states produce identical outputs. But under the inner product with a measurement vector |m⟩ m, the overlap ⟨m|ψ⟩=m1∗c1+m2∗c2 m|ψ =m_1^*c_1+m_2^*c_2 differs between the two cases. In the first, the contributions from dimensions 1 and 2 reinforce; in the second, they cancel. This is interference, and it occurs because complex addition is phase-sensitive. Standard real-valued architectures achieve analogous suppression through dedicated learned mechanisms such as gating (in LSTMs hochreiter1997long and GRUs) or attention weight modulation (in Transformers vaswani2017attention ). These mechanisms require their own parameters and computational overhead. The complex state space provides suppression through the geometry of complex addition itself, without requiring additional parameters, though only when paired with a phase-sensitive readout such as the Born rule defined in Section 3.5. The connection to language modeling is as follows. When the model processes an ambiguous prefix, the state |ψ(t)⟩ ψ(t) distributes amplitude across latent dimensions associated (through training) with different interpretations. The relative phases between these amplitudes record how the interpretations relate to one another: whether they are compatible, conflicting, or independent. As subsequent tokens arrive, the evolution law (Section 3.2) rotates these phases, and the output rule (Section 3.5) converts the resulting interference pattern into token probabilities. The unit-norm constraint ensures that strengthening one interpretation necessarily weakens others: probability is conserved and redistributed, not created or destroyed. The conserved probability currents derived in Section 4 make this redistribution process visible. The basis vectors |j⟩ j do not correspond to pre-assigned linguistic categories. They are abstract directions in ℂNC^N that acquire semantic content through training, in the same way that hidden units in a feedforward network develop interpretable roles through optimization. What distinguishes this state space from a standard real-valued hidden vector is the geometry: the state lives on the complex unit sphere S2N−1⊂ℂNS^2N-1 ^N, which is a (2N−1)(2N-1)-dimensional manifold. Each complex coordinate contributes two real degrees of freedom (magnitude and phase), so two points on S2N−1S^2N-1 can share identical coordinate magnitudes while differing in their phases, and both the dynamics and the output rule are sensitive to these differences. The separation theorem in Section 5 formalizes this by showing that, for a specific family of disambiguation tasks, the phase degrees of freedom enable a complex state of dimension N with Born-rule output to encode relationships that require dimension Ω(N2) (N^2) in a real-valued model with affine-softmax readout. 3.2 Hamiltonian Decomposition Having defined the state space, we need a rule for how the state evolves as the model processes tokens. This rule must satisfy three properties identified in the preceding sections. First, it must preserve the unit norm of the state at every step, since both the probabilistic interpretation of the amplitudes and the Born-rule output (Section 3.5) depend on ‖ψ(t)‖=1\|ψ(t)\|=1. Second, it must allow the input token at each step to shape the dynamics, not merely shift the state by an additive bias, as in the input injection Wht+VxtWh_t+Vx_t used by unitary RNNs (Section 2), but control the axis and rate of the state’s rotation on the unit sphere. Third, it should accommodate multiple timescales of variation, so that different latent dimensions can be allocated to features that change at different rates. We govern the state evolution by the time-dependent Schrödinger equation, idt|ψ(t)⟩=H(t)|ψ(t)⟩,i ddt ψ(t)=H(t) ψ(t), (4) where H(t)H(t) is a Hermitian operator satisfying H(t)=H(t)†H(t)=H(t) . The Hermitian property is the single structural constraint on the dynamics, and it is sufficient to guarantee norm preservation. To verify this, differentiate the norm: dt⟨ψ|ψ⟩=⟨dψdt|ψ⟩+⟨ψ|dψdt⟩=⟨ψ|(iH)†|ψ⟩+⟨ψ|(−iH)|ψ⟩=i⟨ψ|H|ψ⟩−i⟨ψ|H|ψ⟩=0. ddt ψ|ψ = dψdt |ψ + ψ | dψdt = ψ|(iH) |ψ + ψ|(-iH)|ψ =i ψ|H|ψ -i ψ|H|ψ =0. (5) The time-evolution operator (t2,t1)U(t_2,t_1) mapping |ψ(t1)⟩ ψ(t_1) to |ψ(t2)⟩ ψ(t_2) is therefore unitary (†=IU U=I): the state rotates on the unit sphere but never drifts off it. This addresses a stability concern raised in Section 2: Haber and Ruthotto haber2017stable showed that norm-preserving dynamics provide the strongest stability guarantees for deep networks, but their antisymmetric-weight prescription only approximately preserves the norm under Euler discretization. Here, the Hermitian constraint guarantees norm preservation at the continuous level, and the Cayley discretization in Section 3.4 will maintain it at the discrete level as well. The practical consequence is that no gradient clipping, weight normalization, or spectral regularization is needed to keep the state norm stable. We decompose the Hamiltonian into two terms with distinct roles: H(t)=H0+Hint(t).H(t)=H_0+H_int(t). (6) The free Hamiltonian H0H_0. The first term is a diagonal matrix with learnable real entries: H0=diag(λ0,λ1,…,λN−1),λj∈ℝ.H_0=diag( _0, _1,…, _N-1), _j . (7) When the interaction is absent (Hint(t)=0H_int(t)=0), each amplitude evolves independently as cj(t)=cj(0)e−iλjtc_j(t)=c_j(0)\,e^-i _jt, which is a rotation in the complex plane at angular frequency λj _j. The parameters λj\ _j\ define a learnable spectrum of timescales: dimensions with large |λj|| _j| oscillate rapidly, while dimensions with small |λj|| _j| change slowly. Because H0H_0 is diagonal, it induces no coupling between dimensions; its role is to establish the baseline temporal structure that the interaction term will modulate. The motivation for separating out H0H_0 is both conceptual and computational. Conceptually, the spectrum λj\ _j\ provides the model with a bank of oscillators at different frequencies, analogous to a Fourier basis for temporal variation. When the model is trained on language data, the optimization can assign fast frequencies to dimensions that track rapidly varying features (such as local syntactic constraints that change with every token) and slow frequencies to dimensions that track gradually varying features (such as topic or discourse structure that persists over many sentences). Whether the trained model actually organizes its frequencies this way is an empirical question we address in Section 6.2; the architecture provides the capacity for this organization but does not enforce it. Computationally, factoring out H0H_0 enables the interaction-picture transformation (Section 3.3), which removes the known free oscillations from the evolution operator and improves the numerical accuracy of the Cayley discretization by reducing the effective operator norm that governs the integration error. This separation also connects to state-space models. As noted in Section 2, when Hint(t)=0H_int(t)=0 the Schrödinger equation reduces to N decoupled oscillators dtcj=−iλjcj ddtc_j=-i _jc_j, each evolving at a learned frequency. This is structurally identical to the diagonalized form of S4 gu2021efficiently , except that in our framework, the complex state variables are semantically load-bearing (i.e., their phases encode relational information that the output rule will convert into token probabilities) rather than an implementation artifact of a diagonalization. The interaction Hamiltonian Hint(t)H_int(t). The second term depends on the input token observed at time t and is responsible for all coupling between latent dimensions and all input-driven changes to the state trajectory. We construct Hint(t)H_int(t) to be Hermitian by design, regardless of the values produced by the underlying neural network. A network gθg_θ receives as input the concatenation of the token embedding Embed(xt)∈ℝdEmbed(x_t) ^d and the current interaction-picture state represented as real and imaginary parts [Re(cI(t)),Im(cI(t))]∈ℝ2N[Re(c_I(t)),\,Im(c_I(t))] ^2N (Section 3.3 defines this representation; the interaction-picture amplitudes cI(t)c_I(t) differ from the Schrödinger-picture amplitudes c(t)c(t) by the known phase factors eiλjte^i _jt, which gθg_θ can absorb into its learned parameters). The network gθg_θ is a standard feedforward neural network, distinct from, and not to be confused with, the attention weight matrices of Transformer models, that outputs two objects: a complex matrix Φ(t)∈ℂN×r (t) ^N× r and a real vector δ(t)∈ℝNδ(t) ^N. The interaction Hamiltonian is then Hint(t)=Φ(t)Φ(t)†+diag(δ(t)).H_int(t)= (t)\, (t) +diag(δ(t)). (8) This construction guarantees Hermiticity through two separate mechanisms. The matrix Φ(t)Φ(t)† (t) (t) is Hermitian because (ΦΦ†)†=ΦΦ†( ) = for any complex matrix Φ , and it is positive semidefinite with rank at most r, where r≪Nr N is a hyperparameter. The diagonal matrix diag(δ(t))diag(δ(t)) is Hermitian because it is real and diagonal. Their sum is therefore Hermitian for any output of gθg_θ, with no additional projection or correction step. The rank parameter r controls a tradeoff between expressivity and efficiency. An unrestricted N×N× N Hermitian matrix has N2N^2 real degrees of freedom, which would be expensive for gθg_θ to produce at each time step. The low-rank factorization reduces this to 2Nr2Nr real degrees of freedom (since Φ has NrNr complex entries, each with two real components). Each column ϕa(t)∈ℂN _a(t) ^N of Φ(t) (t) specifies a coupling pattern across all N dimensions: the rank-one matrix ϕaϕa† _a _a has entry [ϕa]j[ϕa]k∗[ _a]_j\,[ _a]_k^* at position (j,k)(j,k), coupling dimensions j and k with a strength and phase determined by the amplitudes and phases of ϕa _a at those positions. The full off-diagonal interaction ∑a=1rϕaϕa† _a=1^r _a _a is a superposition of r such coupling patterns. The diagonal term diag(δ(t))diag(δ(t)) supplements this by allowing the network to shift the effective oscillation frequency of individual dimensions in response to the current token, without introducing additional off-diagonal coupling. Together, the two terms give Hint(t)H_int(t) the structure needed to both redirect amplitude between dimensions (via off-diagonal coupling) and modulate their oscillation rates (via diagonal shifts). State-dependent Hamiltonians and nonlinearity. The fact that gθg_θ takes the current state as input, in addition to the token embedding, is a deliberate design choice that determines the model’s expressivity class. If gθg_θ depended only on the token embedding, the Hamiltonian at each step would be a fixed function of the input token, independent of the state. The evolution over a full sequence would then be a composition of token-determined unitary operators applied to the initial state, a sequence-to-state map that is linear in the initial state. This would place the model in the same expressivity class as linear state-space models such as S4 and Mamba gu2023mamba , whose fundamental limitation (discussed in Section 2) is that two tokens can only interact through accumulated products of linear transition matrices. By conditioning gθg_θ on the current state, the Hamiltonian at step t depends on the outcome of all previous steps, and the overall map from input sequence to state trajectory becomes nonlinear. The nonlinearity resides entirely in the map from (token, state) pairs to Hamiltonians; the evolution under any single fixed H(t)H(t) remains a linear unitary map. This separation is the mechanism by which the architecture achieves nonlinear input-output behavior while preserving the norm-conservation guarantee. The situation is analogous to a time-varying linear system in control theory: at each instant the dynamics are linear, but because the transition operator changes at each step as a function of both the input and the state, the composite behavior across steps is nonlinear. This arrangement also contrasts with the input mechanism used by unitary RNNs. As discussed in Section 2, both Arjovsky et al. arjovsky2016unitary and Wisdom et al. wisdom2016full inject the input as an additive bias: ht+1=Wht+Vxth_t+1=Wh_t+Vx_t, where W is a fixed unitary matrix. The unitary part of the dynamics is independent of the input, and the input only displaces the state after the rotation. In our framework, the input token parameterizes the Hamiltonian Hint(t)H_int(t), which is the generator of the rotation itself. The input thereby controls which dimensions couple and at what rate, rather than adding a post-rotation offset. For language modeling, this distinction matters: when a disambiguating token like “steep” follows the prefix “The bank was,” it should not merely shift the state toward one interpretation but should restructure the interference pattern among all active interpretations. An input-dependent Hamiltonian achieves this because the relative phases accumulated under Hint(t)H_int(t) determine which components interfere constructively and which cancel in the subsequent Born-rule readout. 3.3 The Interaction Picture The free Hamiltonian H0H_0 produces phase rotations at rates λj _j that persist regardless of the input. These oscillations are known analytically (i.e., the solution under H0H_0 alone is cj(t)=cj(0)e−iλjtc_j(t)=c_j(0)\,e^-i _jt), and they contribute to the operator norm of the full Hamiltonian H(t)=H0+Hint(t)H(t)=H_0+H_int(t) even when they carry no information about the input-driven dynamics. This matters for numerical accuracy: the Cayley discretization introduced in Section 3.4 is a second-order integrator whose local truncation error scales with the cube of the operator norm of the Hamiltonian being integrated hairer2006geometric . When some free frequencies λj _j are large, the norm ‖H(t)‖\|H(t)\| is dominated by H0H_0 rather than by the input-driven interaction Hint(t)H_int(t), and the integration error is correspondingly inflated. The interaction picture is a change of variables, standard in quantum mechanics sakurai1994modern , that factors out the free oscillations analytically and leaves an evolution equation driven only by the interaction. By removing H0H_0 from the Hamiltonian that the integrator must approximate, the interaction picture reduces the effective operator norm to ‖Hint(t)‖\|H_int(t)\|, which depends only on the input-driven coupling strength and is independent of the free frequencies. This improves the accuracy of the Cayley discretization without affecting its unitarity guarantee, which holds unconditionally. Define the interaction-picture state by |ψ(t)⟩=e−iH0t|ψI(t)⟩, ψ(t)=e^-iH_0t _I(t), (9) so that |ψI(t)⟩=eiH0t|ψ(t)⟩ _I(t)=e^iH_0t ψ(t) is the original state with the free rotations undone. Since H0H_0 is diagonal, this transformation acts componentwise: the j-th amplitude in the interaction picture is [cI(t)]j=cj(t)eiλjt[c_I(t)]_j=c_j(t)\,e^i _jt, which removes the free-oscillation factor from each component. The forward pass of the model operates on cI(t)c_I(t) throughout; the Schrödinger-picture state c(t)c(t) can be recovered at any time via (9) if needed, but it is not maintained explicitly during inference. As described in Section 3.2, the network gθg_θ receives the interaction-picture amplitudes cI(t)c_I(t) as input; since gθg_θ is a universal function approximator, it can learn to account for the known phase relationship cj(t)=e−iλjt[cI(t)]jc_j(t)=e^-i _jt\,[c_I(t)]_j between the two representations. To derive the evolution equation for |ψI(t)⟩ _I(t), substitute the ansatz (9) into the Schrödinger equation (4). The left-hand side becomes idt(e−iH0t|ψI⟩)=H0e−iH0t|ψI⟩+ie−iH0tdt|ψI⟩,i ddt (e^-iH_0t _I )=H_0e^-iH_0t _I+i\,e^-iH_0t ddt _I, (10) and the right-hand side is (H0+Hint(t))e−iH0t|ψI⟩(H_0+H_int(t))\,e^-iH_0t _I. The H0H_0 terms cancel on both sides. Left-multiplying the remainder by eiH0te^iH_0t yields idt|ψI(t)⟩=Hint,I(t)|ψI(t)⟩,i ddt _I(t)=H_int,I(t) _I(t), (11) where the interaction-picture Hamiltonian is Hint,I(t)=eiH0tHint(t)e−iH0t.H_int,I(t)=e^iH_0t\,H_int(t)\,e^-iH_0t. (12) Since H0H_0 is diagonal, this conjugation acts on individual matrix entries as [Hint,I(t)]jk=[Hint(t)]jkei(λj−λk)t.[H_int,I(t)]_jk=[H_int(t)]_jk\,e^i( _j- _k)t. (13) The diagonal entries (j=kj=k) are unchanged, while each off-diagonal entry acquires a time-dependent phase oscillating at the frequency difference λj−λk _j- _k between the two dimensions it couples. This frequency-dependent modulation has consequences for how latent dimensions interact over multiple time steps. Consider two dimensions j and k with similar natural frequencies (λj≈λk _j≈ _k). The factor ei(λj−λk)te^i( _j- _k)t oscillates slowly in discrete time, so the coupling [Hint(t)]jk[H_int(t)]_jk passes through to the interaction picture with little modification across consecutive tokens. Now consider two dimensions with very different frequencies (|λj−λk|| _j- _k| large relative to 1/Δt1/ t). The rapidly oscillating phase means that the effective coupling alternates in sign from one token to the next, making it difficult for the optimization to learn consistent amplitude transfer between these dimensions. The result is an inductive bias toward frequency-selective coupling: the optimization naturally favors amplitude exchange between dimensions on similar timescales and penalizes exchange between dimensions on very different timescales. This is analogous to the resonance phenomenon in coupled oscillators, where energy transfers most efficiently between oscillators whose natural frequencies are close. We emphasize that this is a property of the architecture’s inductive bias over training, not a per-step dynamical effect: at any single time step, the Cayley update applies the full interaction-picture Hamiltonian Hint,I(t)H_int,I(t) regardless of the magnitude of the frequency differences. Two properties of this transformation are important for the rest of the construction. First, it is an exact change of variables, not an approximation. The state |ψI(t)⟩ _I(t) contains the same information as |ψ(t)⟩ ψ(t) and can be recovered via (9) at any time. Second, |ψI(t)⟩ _I(t) is constant whenever Hint(t)=0H_int(t)=0: if no input arrives, the interaction-picture state does not change. This property is what makes the separation computationally useful: the Cayley integrator in Section 3.4 operates on equation (11), whose right-hand side has operator norm ‖Hint,I(t)‖=‖Hint(t)‖\|H_int,I(t)\|=\|H_int(t)\| (since unitary conjugation preserves the operator norm). The integration error therefore depends on the strength of the input-driven interaction, not on the potentially much larger free frequencies λj\ _j\. 3.4 Cayley Discretization The continuous evolution equation (11) must be converted to a discrete update for implementation on digital hardware. The choice of discretization is consequential because not all numerical integrators preserve unitarity, and any norm drift that occurs at a single step accumulates over the length of the sequence. To illustrate the problem, consider forward Euler: replace idt|ψI⟩=Hint,I|ψI⟩i ddt _I=H_int,I _I with |ψt+1⟩=(I−iΔtHint,I(t))|ψt⟩ _t+1=(I-i t\,H_int,I(t)) _t. The update matrix I−iΔtHint,I-i t\,H_int,I is not unitary for any nonzero Hint,IH_int,I and Δt t. To see this, note that for an eigenvector of Hint,IH_int,I with eigenvalue μ, the update multiplies its norm by |1−iΔtμ|=1+Δt2μ2>1|1-i t\,μ|= 1+ t^2μ^2>1, so the state norm grows at each step. Over a sequence of T tokens, this growth compounds, and for a language model processing sequences of thousands of tokens the accumulated drift is large enough to invalidate both the probabilistic interpretation of the state and the Born-rule output. Higher-order explicit methods (Runge–Kutta, Adams–Bashforth) suffer from the same structural deficiency: they approximate the unitary evolution operator with a polynomial in ΔtHint,I t\,H_int,I, and no polynomial of finite degree other than the trivial identity is unitary hairer2006geometric . We instead discretize using the Cayley transform, which arises from the implicit midpoint rule, also known as the Crank–Nicolson scheme crank1947practical ; hairer2006geometric ; iserles2000lie . This is the standard structure-preserving discretization for systems whose continuous-time flow is unitary (or more generally symplectic), and it is the canonical choice within the field of geometric numerical integration leimkuhler2004simulating . Evaluate the right-hand side of (11) at the midpoint 12(|ψt+1⟩+|ψt⟩) 12( _t+1+ _t) rather than at |ψt⟩ _t: |ψt+1⟩−|ψt⟩Δt=−iHint,I(t)|ψt+1⟩+|ψt⟩2. _t+1- _t t=-i\,H_int,I(t)\, _t+1+ _t2. (14) Rearranging to collect |ψt+1⟩ _t+1 on the left: (I+iΔt2Hint,I(t))|ψt+1⟩=(I−iΔt2Hint,I(t))|ψt⟩. (I+ i t2\,H_int,I(t) ) _t+1= (I- i t2\,H_int,I(t) ) _t. (15) The discrete update operator is therefore W(t)=(I+iΔt2Hint,I(t))−1(I−iΔt2Hint,I(t)),W(t)= (I+ i t2\,H_int,I(t) )^-1 (I- i t2\,H_int,I(t) ), (16) which is the Cayley transform of the skew-Hermitian matrix K=iΔt2Hint,I(t)K= i t2\,H_int,I(t). Proposition 3.1 (Unitarity of the Cayley update). If Hint,I(t)H_int,I(t) is Hermitian, then W(t)W(t) is unitary for every Δt>0 t>0. Proof. Define K=iΔt2Hint,I(t)K= i t2\,H_int,I(t). Because Hint,I(t)H_int,I(t) is Hermitian, K†=−iΔt2Hint,I(t)=−K =- i t2\,H_int,I(t)=-K, so K is skew-Hermitian. The update operator is W=(I+K)−1(I−K)W=(I+K)^-1(I-K). For the adjoint, W†=(I−K†)[(I+K†)]−1=(I+K)(I−K)−1W =(I-K )[(I+K )]^-1=(I+K)(I-K)^-1. Now observe that (I+K)(I−K)=I−K2=(I−K)(I+K)(I+K)(I-K)=I-K^2=(I-K)(I+K), so (I+K)(I+K) and (I−K)(I-K) commute and therefore so do their inverses. Computing the product: W†W W W =(I+K)(I−K)−1(I+K)−1(I−K) =(I+K)(I-K)^-1(I+K)^-1(I-K) =(I+K)(I+K)−1(I−K)−1(I−K)=I.∎ =(I+K)(I+K)^-1(I-K)^-1(I-K)=I. (17) The unitarity of W(t)W(t) holds for every step size Δt t and every spectrum of Hint,I(t)H_int,I(t). There is no stability condition requiring Δt t to be small relative to the eigenvalues of Hint,IH_int,I for the purpose of norm preservation: the state norm is exactly 11 after every step, regardless of the Hamiltonian’s magnitude. After T steps, ‖ψT‖=‖ψ0‖=1\| _T\|=\| _0\|=1 up to floating-point rounding errors, which we correct with a single renormalization as a numerical safeguard rather than an algorithmic necessity. Unitarity versus accuracy. It is important to distinguish the unconditional norm-preservation guarantee from the question of integration accuracy. The Cayley transform is a second-order integrator: its local truncation error (the difference between the Cayley update and the exact solution of (11) over one step) is O(Δt3‖Hint,I‖3)O( t^3\|H_int,I\|^3) hairer2006geometric . The norm of the error is controlled by ‖Hint,I‖\|H_int,I\|, which is why the interaction picture (Section 3.3) is beneficial: it removes the free-frequency contribution ‖H0‖\|H_0\| from the operator whose norm governs the error. However, the error is nonzero for any Δt>0 t>0 and any nonzero Hint,IH_int,I. Over T steps, phase errors accumulate: the discrete trajectory |ψI(t)⟩t=0T\ _I(t)\_t=0^T deviates from the exact continuous-time trajectory by an amount that grows with T. This phase drift does not violate any property used by the model. The Born-rule output at each step depends on the current discrete state, not on its fidelity to the continuous-time solution; the training objective (Section 3.6) optimizes the parameters of the discrete model directly. The continuous-time Schrödinger equation (4) serves as the architectural motivation for the update rule, not as a target that the discrete model must track. The guarantees that carry over exactly from continuous to discrete time are norm preservation (Proposition 3.1), the validity of the Born-rule output, and the conservation properties of the probability currents (Section 4). Trajectory accuracy is a separate property that depends on Δt t and ‖Hint,I‖\|H_int,I\|. It is worth noting that Wisdom et al. wisdom2016full also use the Cayley transform in their Full-Capacity Unitary RNN, but in a different role: they use it to parameterize a single fixed unitary transition matrix from a learnable skew-Hermitian matrix. In our framework, the Cayley transform is applied at each time step to a different interaction-picture Hamiltonian Hint,I(t)H_int,I(t), producing a time-varying sequence of unitary updates rather than a single fixed unitary matrix. Both uses exploit the same algebraic fact, namely that the Cayley transform maps skew-Hermitian matrices to unitary matrices, but ours applies it as a discretization of continuous dynamics (in the spirit of the Crank–Nicolson method crank1947practical ) rather than as a static parameterization. In implementation, we solve the linear system (15) at each step using a batched direct solver rather than explicitly forming the inverse (I+iΔt2Hint,I)−1(I+ i t2\,H_int,I)^-1. The low-rank-plus-diagonal structure of Hint,IH_int,I (inherited from the construction in equation (8)) can be exploited via the Woodbury matrix identity to reduce the cost from O(N3)O(N^3) to O(Nr2)O(Nr^2) per step. The Woodbury factorization involves inverting an r×r× r matrix G=Ir+iΔt2Φ~†D−1Φ~G=I_r+ i t2\, D^-1 , where D is the diagonal part of the coefficient matrix (Section 6.1). The conditioning of G depends on the spectrum of Φ~†D−1Φ~ D^-1 , which in turn depends on the current interaction Hamiltonian. For typical Hamiltonian magnitudes encountered during training, G is well conditioned because ‖Hint,I‖\|H_int,I\| is moderate and D is close to the identity. If training drives ‖Hint,I‖\|H_int,I\| to large values, the conditioning of G can degrade; monitoring the condition number of G during training provides a diagnostic for this scenario. 3.5 Born-Rule Decoding At each time step, the model must produce a probability distribution over a vocabulary of V tokens. The standard approach in neural language models is to compute a linear projection of the hidden state followed by softmax normalization. Section 2 discussed a limitation of this approach identified by Yang et al. yang2018breaking : when the hidden dimension N is smaller than the vocabulary size V, the matrix of log-probabilities across all contexts has rank at most N, creating a bottleneck that prevents the model from representing certain distributions regardless of its parameter values. We now define an alternative decoding mechanism whose algebraic structure differs from the softmax readout in a way that is directly relevant to this bottleneck and to the interference mechanism described throughout the paper. We associate with each vocabulary item k∈0,1,…,V−1k∈\0,1,…,V-1\ a learnable measurement vector |mk⟩∈ℂN m_k ^N and define the probability of token k as p(k∣ψ(t))=|⟨mk|ψ(t)⟩|2.p(k ψ(t))= | m_k|ψ(t) |^2. (18) This is a quadratic output rule: the probability depends on the squared modulus of a complex inner product between the state and a learned vector. The same mathematical form appears in quantum mechanics as the Born rule, and we adopt that name to make the connection to the quantum cognition literature (Section 2.6) explicit, while emphasizing that the justification for using this output rule in our framework is algebraic (quadratic access to phase information), not physical. The inner product ⟨mk|ψ(t)⟩=∑j=0N−1[mk]j∗cj(t) m_k|ψ(t) = _j=0^N-1[m_k]_j^*\,c_j(t) is a complex number whose squared modulus gives the output probability. For these probabilities to form a valid distribution over the vocabulary, they must sum to one. Collecting the measurement vectors as columns of a matrix M=[|m0⟩,…,|mV−1⟩]∈ℂN×VM=[ m_0,…, m_V-1] ^N× V, the condition ∑kp(k|ψ)=1 _kp(k|ψ)=1 for all unit-norm |ψ⟩ ψ is equivalent to ∑k=0V−1|mk⟩⟨mk|=MM†=IN. _k=0^V-1 m_k m_k=M =I_N. (19) This is the resolution-of-identity condition: the outer products of the measurement vectors must sum to the identity on ℂNC^N. When V≥NV≥ N (which is always the case in language modeling, where typical vocabularies contain tens of thousands of tokens while latent dimensions are in the hundreds), this requires M to have orthonormal rows. Enforcing the resolution of identity. In implementation, we store M as a learnable complex matrix and enforce the row-orthonormality constraint MM†=INMM =I_N via QR decomposition at each forward pass: we compute the thin QR factorization of M†∈ℂV×NM ^V× N and replace M†M with the Q factor, whose columns are orthonormal. This is a projection onto the Stiefel manifold St(N,V)=X∈ℂV×N:X†X=INSt(N,V)=\X ^V× N:X X=I_N\ absil2008optimization . Gradients of the loss with respect to the pre-projection parameter matrix flow through the QR decomposition via the chain rule; modern automatic differentiation frameworks support differentiation through QR townsend2016pymanopt . The resulting optimization trajectory is not identical to Riemannian gradient descent on the Stiefel manifold absil2008optimization (which would require a retraction and a metric correction), but it is a well-studied approximation that has been used successfully in prior work on orthogonal and unitary parameterizations wisdom2016full . An alternative approach would be to parameterize M directly via a product of Householder reflections or Givens rotations, which would eliminate the need for the QR projection; we use QR for simplicity and defer the comparison of manifold optimization strategies to future work. To verify that the Born rule produces valid probabilities under this constraint, compute ∑kp(k|ψ)=∑k⟨ψ|mk⟩⟨mk|ψ⟩=⟨ψ|(∑k|mk⟩⟨mk|)|ψ⟩=⟨ψ|IN|ψ⟩=‖ψ‖2=1, _kp(k|ψ)= _k ψ|m_k m_k|ψ = ψ| ( _k|m_k m_k| )|ψ = ψ|I_N|ψ =\|ψ\|^2=1, (20) where the last step uses the unit-norm property of the state. This derivation makes explicit the interdependence between the components of the architecture: the Born rule produces valid probabilities because the Cayley discretization preserves the unit norm that the Hermitian Hamiltonian conserves. If the discretization introduced even a small norm drift at each step, the output probabilities would no longer sum to one and a corrective normalization would be needed, defeating the purpose of the construction. Algebraic structure of the output. The algebraic structure of the Born-rule output differs from the softmax readout in a way that bears directly on the expressivity of the model. In a real-valued model with hidden state h∈ℝNh ^N and output weight matrix Wo∈ℝV×NW_o ^V× N, the log-probability of token k is logp(k)=wk⊤h−logZ p(k)=w_k h- Z, which depends on the state only through the linear projection wk⊤hw_k h. The matrix of these projections across all possible states and tokens has rank at most N yang2018breaking . Under the Born rule, the probability is p(k|ψ)=|∑j[mk]j∗cj|2=∑j∑j′[mk]j∗[mk]j′cj∗cj′.p(k|ψ)= | _j[m_k]_j^*\,c_j |^2= _j _j [m_k]_j^*\,[m_k]_j \,c_j^*\,c_j . (21) This is a quadratic form in the amplitudes, involving both the N diagonal terms |[mk]j|2|cj|2|[m_k]_j|^2|c_j|^2 (which depend only on magnitudes) and the (N2) N2 off-diagonal cross terms (which depend on relative phases). Each cross term has the form [mk]j∗[mk]j′cj∗cj′=|[mk]j||[mk]j′|rjrj′ei(αj′−αj+θj′−θj),[m_k]_j^*\,[m_k]_j \,c_j^*\,c_j =|[m_k]_j|\,|[m_k]_j |\,r_j\,r_j \,e^i( _j - _j+ _j - _j), (22) where αj=arg[mk]j _j= [m_k]_j and θj=argcj _j= c_j. The real part of this expression contributes to p(k|ψ)p(k|ψ), and its sign depends on the combined phase αj′−αj+θj′−θj _j - _j+ _j - _j. Two states with identical magnitudes but different relative phases therefore produce different token probabilities, an impossibility under a linear readout of a real-valued state. The quadratic structure means the Born-rule output depends on N2N^2 features of the state (the N squared magnitudes and N(N−1)N(N-1) real parameters in the off-diagonal cross terms), compared to the N features accessible to a linear readout. The formal version of this observation is Lemma 5.7 in Section 5, and the resulting dimensional separation is Theorem 5.10. We note that this comparison is specific to the affine-softmax readout used by standard architectures; a real-valued model equipped with a polynomial readout of degree two or higher could access a larger feature space, though such readouts are not standard in the architectures surveyed in Section 2. This quadratic structure is where the interference mechanism becomes concrete at the output level. The cross terms can contribute positively or negatively to the probability of token k, depending on the phase alignment. When the Hamiltonian evolution causes the phases θj _j to rotate such that previously constructive cross terms become destructive for a given token k, the probability p(k|ψ)p(k|ψ) decreases, not because a gate has been driven to saturation, but because the amplitudes are now interfering against that token. Conversely, when the phases align constructively for a different token k′k , its probability increases. This is the Born-rule realization of the disambiguation process described in the Introduction: the model suppresses one interpretation and reinforces another by rotating phases in its latent space, and the quadratic output converts these rotations into probability changes. The connection to the quantum cognition literature (Section 2.6) is also direct: the Born rule (18) is the same squared-amplitude measurement postulate that Busemeyer and Bruza busemeyer2012quantum use to explain interference effects in human judgment. In their framework, the measurement vectors correspond to possible responses in a psychological experiment and the state encodes the subject’s cognitive disposition; in ours, the measurement vectors correspond to vocabulary tokens and the state encodes the model’s belief about the next token given the preceding context. Whether the structural similarity between these two settings reflects a deeper connection between language statistics and cognitive interference patterns, or is merely a mathematical coincidence, is an empirical question that the present theoretical framework cannot answer. 3.6 Initialization and Training Objective The initial state |ψ(0)⟩ ψ(0) must be a unit-norm vector in ℂNC^N defined independently of the input sequence, since no tokens have been observed at t=0t=0. We parameterize it using two real-valued vectors ,∈ℝNa,b ^N as |ψ(0)⟩=+i‖+i‖. ψ(0)= a+ib\|a+ib\|. (23) The division by the norm ensures ‖ψ(0)‖=1\|ψ(0)\|=1 for any values of a and b, so standard gradient-based optimization kingma2015adam can update these parameters freely without risking an invalid initial state. The parameters are shared across all training sequences and optimized jointly with the rest of the model. After normalization, the magnitudes |aj+ibj|/‖+i‖|a_j+ib_j|/\|a+ib\| determine how much initial weight the model places on each latent dimension, and the phases arg(aj+ibj) (a_j+ib_j) set the initial phase relationships between dimensions before any tokens are processed. Together, these constitute a learned prior: the optimization discovers an initial distribution of amplitude and an initial pattern of phase relationships that provide a useful starting configuration for processing arbitrary input sequences. The training objective is the standard negative log-likelihood for autoregressive language modeling, computed using the Born-rule probabilities. Given a sequence x0,x1,…,xT−1x_0,x_1,…,x_T-1 with targets yt=xt+1y_t=x_t+1, the loss is ℒ=−∑t=0T−1logp(yt∣ψ(t))=−∑t=0T−1log|⟨myt|ψ(t)⟩|2.L=- _t=0^T-1 p(y_t ψ(t))=- _t=0^T-1 | m_y_t|ψ(t) |^2. (24) This objective drives the model to adjust its parameters so that, at each time step, the state |ψ(t)⟩ ψ(t) has large overlap (in the Born-rule sense) with the measurement vector |myt⟩ m_y_t corresponding to the correct next token. Minimizing (24) simultaneously shapes the initial state (determining the starting point of the trajectory), the free frequencies (determining the baseline timescale structure), the Hamiltonian generator (determining how inputs steer the evolution), and the measurement vectors (determining which state configurations correspond to which tokens). The full set of learnable parameters is: the initial-state vectors a and b, the free frequencies λjj=0N−1\ _j\_j=0^N-1, the measurement matrix M, and the weights θ of the neural network gθg_θ. Gradients of ℒL with respect to all parameters are computed by backpropagation through the sequence of Cayley integration steps. Each step involves solving the linear system (15), and the backward pass differentiates through this solve using implicit differentiation: given a system Ax=bAx=b, the gradient with respect to any parameter entering A or b requires solving one additional linear system with the coefficient matrix A†A . This avoids forming or storing the matrix inverse and is supported natively by modern automatic differentiation frameworks. The computational cost per time step is dominated by the linear solve, which scales as O(Nr2)O(Nr^2) when the Woodbury identity exploits the low-rank structure of Hint,IH_int,I. Gradient propagation through the state versus through the parameters. The unitarity of each Cayley step has a direct consequence for gradient propagation through the hidden state. The Jacobian of the Cayley update |ψt+1⟩=W(t)|ψt⟩ _t+1=W(t) _t with respect to the input state |ψt⟩ _t is the unitary matrix W(t)W(t) itself. The gradient of the loss with respect to |ψt⟩ _t is obtained by applying W(t)†=W(t)−1W(t) =W(t)^-1, which preserves the gradient’s norm. Over T−tT-t steps, the gradient passes through T−tT-t unitary matrices, and its norm remains unchanged. This prevents the vanishing and exploding gradient problem for the state gradient pathway vorontsov2017orthogonality , providing the same benefit as the unitary RNNs of Arjovsky et al. arjovsky2016unitary and Wisdom et al. wisdom2016full . However, the gradient of the loss with respect to the parameters of gθg_θ at step t involves an additional pathway: the gradient flows from |ψt+1⟩ _t+1 through W(t)W(t), into Hint,I(t)H_int,I(t), and then through the neural network gθg_θ. The network gθg_θ is an unconstrained feedforward network with no norm-preservation guarantee. Gradients passing through gθg_θ at each step can grow or shrink depending on the condition of gθg_θ’s Jacobian, and this growth or shrinkage is not mitigated by the unitarity of the state dynamics. The unitary guarantee therefore addresses one source of gradient instability (the recurrent state pathway) but not all sources (the parameter pathway through gθg_θ). Standard techniques for stabilizing feedforward networks, such as careful initialization, normalization within gθg_θ, and learning rate scheduling kingma2015adam , remain relevant for the parameter gradients. Therefore, the complete forward pass is as follows: at t=0t=0, the state is initialized via (23). At each subsequent step t, the network gθg_θ generates the interaction Hamiltonian Hint(t)H_int(t) from the current token embedding and current interaction-picture state; the interaction-picture Hamiltonian Hint,I(t)H_int,I(t) is computed via the conjugation (12); the Cayley update (15) advances the interaction-picture state by one step; and the Born rule (18) produces a distribution over the vocabulary from which the loss is computed. The loss (24) aggregates the log-probabilities across all time steps, and its gradients flow backward through the Cayley solves, through the Hamiltonian generator gθg_θ, and into all learnable parameters. The entire pipeline is differentiable end-to-end, with the unitarity of each Cayley step ensuring that gradients propagating along the state pathway neither vanish nor explode, while gradients through the parameter pathway are subject to the same considerations as in any deep feedforward network. (a) Unrolled sequence over T tokens|ψ(0)⟩ ψ(0)learned initx0x_0gθg_θH(0)H(0)W(0)W(0)|ψ(1)⟩ ψ(1)Bornx1x_1gθg_θH(1)H(1)W(1)W(1)Born⋯·s⋯·s|ψT−1⟩ _T\!-\!1xT−1x_T\!-\!1gθg_θH(T−1)H(T-1)W(T−1)W(T-1)|ψ(T)⟩ ψ(T)Bornℒ=−∑tlogp(yt|ψ(t))L=- _t p(y_t|ψ(t)) Key properties: ‖ψ(t)‖=1\|ψ(t)\|=1 at every t ∙ Each H(t)H(t) depends on xtx_t and ψ(t)ψ(t) ∙ Gradient preserved through W(t)†W(t) ∙ Parameters shared: gθg_θ, H0H_0, |mk⟩\|m_k \ (b) Interference — amplitude redistributionBefore “steep”financialriverother1other2W(t)W(t)After “steep”financialriverother1other2J1←0>0J_1← 0>0: amplitude flows financial → river ∙ : ∑j|cj|2=1 _j|c_j|^2=1(c) Born-rule readout|ψ⟩=∑jcj|j⟩ ψ= _jc_j jρ=|ψ⟩⟨ψ|ρ= ψ\! ψquadraticp(k)=tr(Mkρ)p(k)=tr(M_kρ)linear in ρ 2N^2 features: N magnitudes ++ N(N−1)N(N-1) cross-termsp(k|ψ)p(k|ψ) Figure 2: Multi-scale view of the quantum sequence model. (a) The model unrolled over T time steps. At each step the same network gθg_θ generates a token- and state-dependent Hamiltonian; the Cayley update advances the state on the complex unit sphere; and the Born rule reads out token probabilities. The loss aggregates log-probabilities across all steps. (b) Illustration of the interference mechanism. Processing the disambiguating token “steep” after the prefix “The bank was” causes probability to flow from the financial interpretation to the river interpretation via conserved, antisymmetric probability currents Jj←kJ_j← k. (c) Detail of the Born-rule readout. The quadratic map ψ↦ψψ†ψ ψ lifts the N-dimensional complex state into the N2N^2-dimensional space of Hermitian matrices, exposing both magnitude and phase cross-terms to the linear measurement tr(Mkρ)tr(M_kρ). 4 Conserved Probability Currents The derivation in this section mirrors the standard quantum-mechanical treatment of probability current in the Schrödinger equation sakurai1994modern ; griffiths2005introduction . Our contribution is threefold: we interpret the resulting currents as a diagnostic for latent probability flow in a trained language model, we show how the Hamiltonian decomposition H=H0+HintH=H_0+H_int isolates the input-driven component of each current, and we derive an exact discrete-time variant (Section 4.5) consistent with the Cayley update. Section 3 established that the Hermitian constraint on H(t)H(t) guarantees conservation of the total probability ∑j|cj(t)|2=1 _j|c_j(t)|^2=1 at every time step. This global conservation law ensures that the Born-rule output (Section 3.5) produces a valid probability distribution, but it does not describe how the model redistributes amplitude among its latent dimensions as it processes each token. Two models that both conserve total probability may differ entirely in the internal pattern of amplitude transfer they execute at each step. This section derives a finer-grained conservation law, a continuity equation with explicit pairwise flux terms, that resolves the internal flow of probability. The derivation proceeds in continuous time, using the Schrödinger equation (4) that motivates the architecture, and produces probability currents whose structural properties (antisymmetry, conservation) follow algebraically from the Hermiticity of H(t)H(t). Because the implemented model operates in discrete time via the Cayley update (Section 3.4), the continuous-time currents are approximations to the actual discrete probability changes; Section 4.5 quantifies this approximation and defines exact discrete currents that inherit the same structural properties. 4.1 Derivation of the Continuity Equation Define the occupation probability of latent dimension j at time t as pj(t)=|cj(t)|2,p_j(t)=|c_j(t)|^2, (25) where cj(t)=⟨j|ψ(t)⟩c_j(t)= j|ψ(t) is the j-th amplitude of the state vector. The normalization constraint ∑jpj(t)=1 _jp_j(t)=1 ensures that pj(t)j=0N−1\p_j(t)\_j=0^N-1 forms a probability distribution over the N latent dimensions at each time t. We now derive the equation governing how this distribution evolves. The Schrödinger equation (4) in component form reads ic˙j(t)=∑k=0N−1Hjk(t)ck(t),i\, c_j(t)= _k=0^N-1H_jk(t)\,c_k(t), (26) so that c˙j=−i∑kHjkck c_j=-i _kH_jk\,c_k. Taking the complex conjugate and applying the Hermiticity condition Hjk∗=HkjH_jk^*=H_kj yields c˙j∗(t)=i∑k=0N−1Hkj(t)ck∗(t). c_j^*(t)=i _k=0^N-1H_kj(t)\,c_k^*(t). (27) Differentiating the occupation probability (25) via the product rule: p˙j p_j =c˙j∗cj+cj∗c˙j = c_j^*\,c_j+c_j^*\, c_j =(i∑kHkjck∗)cj+cj∗(−i∑kHjkck) = (i _kH_kj\,c_k^* )c_j\;+\;c_j^* (-i _kH_jk\,c_k ) =∑k=0N−1(iHkjck∗cj−iHjkcj∗ck). = _k=0^N-1 (i\,H_kj\,c_k^*\,c_j\;-\;i\,H_jk\,c_j^*\,c_k ). (28) Each summand can be simplified by recognizing its algebraic structure. Using Hkj=Hjk∗H_kj=H_jk^* and the commutativity of scalar multiplication, the first term in each summand satisfies Hkjck∗cj=Hjk∗ck∗cj=(Hjkcj∗ck)∗H_kj\,c_k^*\,c_j=H_jk^*\,c_k^*\,c_j=(H_jk\,c_j^*\,c_k)^*. Setting w=Hjkcj∗ckw=H_jk\,c_j^*\,c_k, the summand becomes i(w∗−w)i(w^*-w). The identity i(w∗−w)=i(−2iIm(w))=2Im(w)i(w^*-w)=i(-2iIm(w))=2Im(w), valid for any w∈ℂw , then gives the continuity equation: dpjdt=∑k=0N−1Jj←k(t), \; dp_jdt= _k=0^N-1J_j← k(t),\; (29) where the probability current from dimension k to dimension j is Jj←k(t)= 2Im(Hjk(t)cj∗(t)ck(t)).J_j← k(t)\;=\;2\,Im\! (H_jk(t)\,c_j^*(t)\,c_k(t) ). (30) The current Jj←k(t)J_j← k(t) is a real number at each time t. When positive, probability is flowing from dimension k into dimension j; when negative, the flow is reversed. The continuity equation (29) states that the rate of change of probability at any dimension is entirely accounted for by the sum of pairwise currents into that dimension from every other dimension. 4.2 Properties of the Probability Current The currents inherit structural properties from the Hermiticity of H(t)H(t) that constrain them to be internally consistent: no current flows from a dimension to itself, and every unit of probability arriving at one dimension departs from exactly one other. Proposition 4.1 (Structural properties of the probability current). The probability currents defined in (30) satisfy the following three properties for all j, k, and t: 1. Antisymmetry: Jj←k(t)=−Jk←j(t)J_j← k(t)=-J_k← j(t). 2. Zero self-current: Jj←j(t)=0J_j← j(t)=0. 3. Global conservation: ∑j=0N−1dpjdt=0 _j=0^N-1 dp_jdt=0. Proof. For property (1), compute Jk←j=2Im(Hkjck∗cj)J_k← j=2Im(H_kj\,c_k^*\,c_j). Since Hkj=Hjk∗H_kj=H_jk^* and ck∗cj=(cj∗ck)∗c_k^*\,c_j=(c_j^*\,c_k)^*, the argument is Hjk∗(cj∗ck)∗=(Hjkcj∗ck)∗H_jk^*\,(c_j^*\,c_k)^*=(H_jk\,c_j^*\,c_k)^*. The identity Im(w∗)=−Im(w)Im(w^*)=-Im(w) then yields Jk←j=−2Im(Hjkcj∗ck)=−Jj←kJ_k← j=-2Im(H_jk\,c_j^*\,c_k)=-J_j← k. For property (2), set k=jk=j in (30): Jj←j=2Im(Hjj|cj|2)J_j← j=2Im(H_j\,|c_j|^2). The diagonal entry HjjH_j of a Hermitian matrix is real, and |cj|2∈ℝ|c_j|^2 , so their product is real and its imaginary part vanishes. For property (3), sum the continuity equation (29) over j: ∑jp˙j=∑j∑kJj←k _j p_j= _j _kJ_j← k. Every pair (j,k)(j,k) with j≠kj≠ k contributes Jj←k+Jk←j=0J_j← k+J_k← j=0 by antisymmetry, and the diagonal terms Jj←jJ_j← j vanish by property (2). The double sum is therefore zero. ∎ Since Jj←j=0J_j← j=0, the continuity equation (29) can equivalently be written as a sum restricted to k≠jk≠ j: dpjdt=∑k≠jJj←k(t). dp_jdt= _k≠ jJ_j← k(t). (31) This form makes explicit that the rate of change of occupation probability at dimension j equals the net inflow from all other dimensions. The (N2) N2 independent currents Jj←k(t)j<k\J_j← k(t)\_j<k (antisymmetry halves the count) form a complete, internally consistent ledger of probability transfers at each instant: every unit of probability that arrives at one dimension is drawn from a specific other dimension, with no residual or unaccounted-for flow. 4.3 Structure of the Currents Under the Hamiltonian Decomposition The decomposition H(t)=H0+Hint(t)H(t)=H_0+H_int(t) introduced in Section 3.2 determines which component of the Hamiltonian drives the probability currents. Because H0=diag(λ0,…,λN−1)H_0=diag( _0,…, _N-1) is diagonal, its off-diagonal entries vanish: [H0]jk=0[H_0]_jk=0 for j≠kj≠ k. Therefore, for any pair of distinct dimensions j≠kj≠ k, Hjk(t)=[Hint(t)]jk,H_jk(t)=[H_int(t)]_jk, (32) and the current between them depends exclusively on the interaction Hamiltonian: Jj←k(t)=2Im([Hint(t)]jkcj∗(t)ck(t))(j≠k).J_j← k(t)=2\,Im\! ([H_int(t)]_jk\,c_j^*(t)\,c_k(t) ) (j≠ k). (33) This has an immediate and important consequence: when no input is driving the interaction (Hint(t)=0H_int(t)=0), every off-diagonal current vanishes, and the occupation probabilities pj(t)\p_j(t)\ are individually conserved. Under H0H_0 alone, the amplitudes rotate in the complex plane at their respective frequencies λj _j, but their squared magnitudes do not change. The free Hamiltonian modulates phases without transferring probability; all probability redistribution is caused by the interaction term Hint(t)H_int(t), which is generated by the neural network gθg_θ from the current token embedding and state (Section 3.2). The currents therefore isolate the causal contribution of each token to the model’s internal dynamics, separating input-driven redistribution from the baseline phase evolution. The low-rank-plus-diagonal construction of Hint(t)H_int(t) from equation (8) further decomposes each current into contributions from individual coupling channels. Since [Hint(t)]jk=∑a=1rΦja(t)Φka∗(t)[H_int(t)]_jk= _a=1^r _ja(t)\, _ka^*(t) for j≠kj≠ k (the diagonal term diag(δ(t))diag(δ(t)) contributes only when j=kj=k), the current separates as Jj←k(t)=2∑a=1rIm(Φja(t)Φka∗(t)cj∗(t)ck(t)).J_j← k(t)=2 _a=1^rIm\! ( _ja(t)\, _ka^*(t)\,c_j^*(t)\,c_k(t) ). (34) Each coupling channel a∈1,…,ra∈\1,…,r\ contributes an independent term to the total current between dimensions j and k. The a-th term transfers probability at a rate determined jointly by two factors: the coupling strengths Φja(t) _ja(t) and Φka(t) _ka(t), which are set by the network gθg_θ in response to the current token and state, and the current amplitudes cj(t)c_j(t) and ck(t)c_k(t), which reflect the model’s accumulated response to all previous tokens. A current can flow between j and k only if both dimensions carry nonzero amplitude and the interaction Hamiltonian couples them; the direction of the flow (from k to j or from j to k) is determined by the sign of Im(ΦjaΦka∗cj∗ck)Im( _ja\, _ka^*\,c_j^*\,c_k), which depends on the relative phases of the coupling and the amplitudes. The state-dependence of the current formula merits emphasis. The same token can trigger different current patterns at different points in a sequence, because the state |ψ(t)⟩ ψ(t) encodes the cumulative effect of all preceding tokens. Two occurrences of the same word in different contexts will produce different currents, even though the interaction Hamiltonian Hint(t)H_int(t) depends on the same token embedding in both cases, because the amplitudes cj(t)c_j(t) and ck(t)c_k(t) entering the current formula differ. This context-sensitivity is not an additional feature that must be engineered; it is an algebraic consequence of the bilinear dependence of (30) on the Hamiltonian and the state. 4.4 Interpretation as a Diagnostic The probability currents provide a mechanistic account of how the model processes each token, distinct from post-hoc interpretability methods in both its derivation and its guarantees. Gradient-based attribution methods simonyan2014deep answer a counterfactual question (how would the output change if this input were perturbed?) by linearizing the model around its current operating point, an approximation whose fidelity depends on the local curvature of the loss landscape. Attention visualization vaswani2017attention reveals which context positions a Transformer attends to but not what computational role the attended content plays. Probing classifiers belinkov2022probing test whether specific information is decodable from a hidden state but do not reveal the process by which that information was incorporated. The probability currents answer a different question: at each time step, where is probability flowing, from which dimensions, to which dimensions, and at what rate? The structural properties established in Proposition 4.1 impose internal consistency constraints that post-hoc attributions are not guaranteed to satisfy: any probability gained by one dimension is explicitly sourced from one or more other dimensions, and the total is exactly conserved. These constraints follow algebraically from the Hermiticity of H(t)H(t) and cannot be violated regardless of the specific parameter values or training procedure. Two distinctions are necessary regarding the scope of these guarantees. First, the continuity equation (29) and the current formula (30) are derived from the continuous-time Schrödinger dynamics, while the implemented model advances in discrete Cayley steps. Section 4.5 below shows that the continuous-time currents approximate the actual discrete probability changes to first order in Δt t, and defines exact midpoint currents that reproduce the discrete changes without approximation error while retaining antisymmetry and conservation. For diagnostic use on the implemented model, the midpoint currents (equation (37)) are the appropriate tool; the continuous-time currents derived in this subsection provide the conceptual framework and the structural guarantees that the midpoint currents inherit. Second, the currents describe how probability redistributes across latent dimensions in the model’s internal state, not how probability redistributes across vocabulary tokens in the output distribution. The connection between internal redistribution and output change is mediated by the Born rule (Section 3.5): a current flowing from dimension k to dimension j increases |cj|2|c_j|^2 and decreases |ck|2|c_k|^2, and also modifies the cross terms cjcj′∗c_jc_j ^* in the Born-rule expansion, both of which affect the output distribution. The currents therefore provide a causal account of one specific aspect of the model’s computation (amplitude redistribution), from which output-level effects can be derived but are not directly read off. The channel decomposition (34) adds a further layer of resolution. It identifies which of the r coupling channels mediates each part of a given probability transfer, linking the redistribution to specific columns of the matrix Φ(t) (t) produced by gθg_θ. After training, the dominant currents triggered by a given token at a given context constitute a token-level summary of the model’s internal response: which dimensions gain amplitude, which lose it, and which coupling pathways carry the transfer. This diagnostic is available at every time step and for every input sequence. It does not require selecting a specific output class with respect to which attribution is computed (as gradient methods do), and it does not require training a separate model (as probing methods do). Its computational cost is O(N2)O(N^2) per time step for the full current matrix, or O(Nr)O(Nr) per time step if only the channel-decomposed currents for a preselected subset of dimension pairs are needed. For the typical regime in which N is in the hundreds and r≪Nr N, this cost is modest relative to the O(Nr2)O(Nr^2) cost of the Cayley update itself. 4.5 Relation to the Discrete Update The probability currents (30) are defined for the continuous-time Schrödinger dynamics. In practice, the model advances in discrete steps via the Cayley update (Section 3.4), and the relationship between the continuous-time currents and the actual discrete probability changes requires specification. This subsection establishes that relationship and defines an exact discrete alternative. After one Cayley step with step size Δt t, the discrete change in occupation probability is Δpj(t)=|cj(t+Δt)|2−|cj(t)|2. p_j(t)\;=\;|c_j(t+ t)|^2-|c_j(t)|^2. (35) The Cayley transform implements the implicit midpoint rule, which is a second-order integrator. The continuous-time currents therefore approximate the discrete changes as Δpj(t)=Δt∑k≠jJj←k(t)+O(Δt2). p_j(t)= t _k≠ jJ_j← k(t)\;+\;O( t^2). (36) The O(Δt2)O( t^2) error arises from the difference between the continuous-time derivative p˙j(t) p_j(t) and the finite-difference quotient Δpj(t)/Δt p_j(t)/ t. For step sizes Δt t of order unity (the natural choice when each discrete step corresponds to one token), the continuous-time currents are first-order accurate estimates of the actual discrete probability transfers. This means the continuous-time currents provide qualitatively correct but quantitatively approximate descriptions of the discrete model’s behavior: they correctly identify the direction and relative magnitude of probability flow between dimensions, but their numerical values may differ from the actual discrete changes by an O(Δt2)O( t^2) correction. Exact discrete probability changes Δpj(t) p_j(t) can always be computed from the Cayley update itself, since the update is explicitly unitary and the occupation probabilities before and after each step are available. However, expressing Δpj p_j in terms of the matrix entries of the Cayley update operator W(t)W(t) yields a sum over products of entries of W(t)W(t) and amplitudes that does not factor into a simple pairwise form with guaranteed antisymmetry. The continuous-time currents remain useful because they provide the pairwise, antisymmetric, conserving decomposition that makes the diagnostic interpretable. Midpoint currents: an exact discrete decomposition. The tension between the appealing structural properties of the continuous-time currents and the fact that the model operates in discrete time can be resolved by defining currents evaluated at the implicit midpoint that the Cayley scheme uses. Specifically, evaluate the continuous-time current formula (30) at the midpoint state 12(|ψ(t+Δt)⟩+|ψ(t)⟩) 12( ψ(t+ t)+ ψ(t)) and the Hamiltonian H(t)H(t), mirroring the implicit midpoint evaluation that defines the Cayley scheme (equation (14)). The resulting midpoint currents are Jj←kmid(t)=2Im(Hjk(t)c¯j∗(t)c¯k(t)),c¯(t)=12(c(t+Δt)+c(t)).J_j← k^mid(t)=2\,Im\! (H_jk(t)\, c_j^*(t)\, c_k(t) ), c(t)= 12 (c(t+ t)+c(t) ). (37) These midpoint currents satisfy antisymmetry by the same argument as Proposition 4.1 (the proof depends only on the Hermiticity of HjkH_jk and the algebraic structure of the imaginary part, both of which hold identically for the midpoint amplitudes c¯j c_j). They also reproduce the exact discrete probability changes without residual error: Δpj(t)=Δt∑k≠jJj←kmid(t). p_j(t)= t _k≠ jJ_j← k^mid(t). (38) To see why, note that the Cayley update is derived precisely by evaluating the right-hand side of the Schrödinger equation at the midpoint (equation (14)). The discrete change Δpj p_j can be written as Δpj=Δt(p˙j)mid p_j= t\,( p_j)_mid, where (p˙j)mid( p_j)_mid is the time derivative of pjp_j evaluated at the midpoint state. The derivation of the continuity equation in Section 4.1 applies identically to any state vector and any Hermitian matrix, so evaluating it at (c¯,H)( c,H) produces (p˙j)mid=∑kJj←kmid( p_j)_mid= _kJ_j← k^mid, and multiplying by Δt t yields (38). The midpoint currents therefore provide an exact, pairwise, antisymmetric decomposition of each discrete update, with no approximation error. The cost of computing them is the same as for the continuous-time currents (O(N2)O(N^2) for the full matrix or O(Nr)O(Nr) for channel-decomposed currents on selected pairs), plus the requirement of having the post-step state |ψ(t+Δt)⟩ ψ(t+ t) available, which it is after the forward pass. For diagnostic purposes on the implemented discrete model, the midpoint currents are the recommended tool: they inherit the structural properties of Proposition 4.1 (antisymmetry, zero self-current, global conservation) while exactly matching the discrete probability changes that the Cayley update produces. Summary of the two current variants. The continuous-time currents Jj←k(t)J_j← k(t) (equation (30)) are derived from the continuous Schrödinger equation and evaluated at the pre-step state c(t)c(t). They provide the conceptual foundation, the structural properties (Proposition 4.1), the Hamiltonian decomposition (Section 4.3), and the channel factorization (equation (34)). They approximate the discrete probability changes to first order, with O(Δt2)O( t^2) error. The midpoint currents Jj←kmid(t)J_j← k^mid(t) (equation (37)) are evaluated at the average of the pre-step and post-step states. They inherit all structural properties, support the same Hamiltonian and channel decompositions (by replacing cj,ckc_j,c_k with c¯j,c¯k c_j, c_k in all formulas), and reproduce the exact discrete probability changes. The choice between the two depends on whether qualitative interpretability (continuous-time currents, available without the post-step state) or quantitative fidelity to the discrete dynamics (midpoint currents, requiring the post-step state) is the priority. Having established the structural and diagnostic properties of the model’s dynamics, we turn in the next section to the question of representational capacity: for a concrete family of tasks, how does the state dimension required by a complex unitary model compare to the dimension required by a real orthogonal model? The formal answer, a quadratic lower bound of Ω(N2) (N^2) on the real state dimension required to match a complex model of dimension N, is proved in full as Theorem 5.10 and Lemma 5.9 of Section 5.4. 5 Expressivity Separation The preceding sections constructed a sequence model whose state evolves unitarily in ℂNC^N and whose output probabilities are extracted by the Born rule. Section 3.5 observed that the Born-rule output is a quadratic function of the complex amplitudes, involving (N2) N2 pairwise interference terms that depend on relative phases. Section 4 showed that the probability currents generated by input tokens redistribute amplitude through these same pairwise channels. These observations suggest that the complex-valued model encodes more information per dimension than a real-valued model with a linear readout, because the quadratic output structure converts N complex coordinates into O(N2)O(N^2) effective features. This section formalizes that suggestion as a separation theorem: we construct a concrete family of disambiguation tasks, prove that a complex unitary model of dimension N solves each task exactly, and prove that any real orthogonal model with an affine-softmax readout requires dimension Ω(N2) (N^2) for the same task. Born-Rule Readout (CUSM)|ψ⟩∈ℂN ψ ^Nρ=|ψ⟩⟨ψ|ρ= ψ\! ψVeronese-type lifting|cj|2⏟N+cj∗cj′⏟N(N−1) |c_j|^2_N+ c_j^*c_j _N(N-1)p(k)=tr(Mkρ)p(k)=tr(M_kρ)quadraticlinear in ρ N2N^2 featuresAffine-Softmax Readout (ROSM)h∈ℝdh ^dzk=wk⊤h+bkz_k=w_k h+b_kaffinelogp(k)=zk−logZ p(k)=z_k- Zp(k)=softmax(z)kp(k)=softmax(z)_klinearrank(logP)≤d+2rank( P)≤ d+2Dimensional comparisonNNCUSM:d≥N2−2d≥ N^2-2ROSM:quadratic gap: Θ(N) (N) Theorem 5.10 (conditional). A CUSM of dimension N solves ND_N exactly. Any ROSM with affine-softmax readout that matches p∗p^* requires d≥N2−2d≥ N^2-2, assuming rank(L∗)=N2rank(L^*)=N^2 (see Remark 5.11). Figure 3: Source of the expressivity separation. Left: The Born-rule readout applies a quadratic (Veronese-type) lifting ψ↦ψψ†ψ ψ , promoting the N-dimensional complex state to the N2N^2-dimensional space of rank-one Hermitian matrices. The measurement tr(Mkρ)tr(M_kρ) is then linear in the lifted space, accessing all N2N^2 features including the (N2) N2 pairwise phase cross-terms that encode interference. Right: The affine-softmax readout computes zk=wk⊤h+bkz_k=w_k h+b_k, which is linear in the d-dimensional real state. The rank constraint (Lemma 5.9) limits the log-probability matrix to rank d+2d+2. Bottom: Matching the N2N^2 features required by the Born-rule target forces d≥N2−2d≥ N^2-2, a quadratic gap over the complex model’s dimension N, under the full-rank condition on L∗L^* stated in Theorem 5.10. Three preliminary remarks are necessary to frame the result precisely. First, the separation established here is between two specific output mechanisms: the Born rule (quadratic in the state) and the affine-softmax readout (linear in the state). The lower bound on the real model’s dimension (Lemma 5.9) relies on the rank constraint imposed by the affine-softmax structure yang2018breaking . A real-valued model equipped with a different readout, for instance a polynomial readout of degree two or higher, would not be subject to the same rank bound and could potentially narrow the gap. The separation should therefore be understood as a comparison between the standard real-valued architecture (orthogonal dynamics with affine-softmax output, as used by all unitary and orthogonal RNNs surveyed in Section 2) and the proposed complex-valued architecture (unitary dynamics with Born-rule output), not as a fundamental limitation of real-valued representations independent of the output mechanism. Second, both model classes defined below (Definitions 5.1 and 5.2) use state-independent transitions: each token maps to a single fixed unitary or orthogonal matrix, with no dependence on the current state. This is strictly weaker than the full architecture of Section 3, which uses state-dependent Hamiltonians to achieve nonlinearity (Section 3.2). The separation theorem proves an advantage for the simplified CUSM, and since the full model subsumes the CUSM as a special case, the advantage carries over; but the theorem does not characterize the representational capacity of the full state-dependent model. Third, the formal task family ND_N defined below uses abstract context and query tokens, not natural-language words. Whether natural language contains subproblems whose structure is captured by ND_N is an empirical question that this section does not address; the task family serves to establish a provable dimensional gap, not to model a specific linguistic phenomenon. 5.1 Model Classes for the Separation For the separation theorem, we compare two simplified model classes that isolate the representational distinction between complex and real state spaces. Both assign one fixed transition matrix per alphabet token, with no state-dependent dynamics. Definition 5.1 (Complex Unitary Sequence Model (CUSM)). A CUSM of dimension N over an alphabet Σ and vocabulary =0,…,V−1V=\0,…,V-1\ consists of: • An initial state |ψ0⟩∈ℂN _0 ^N with ‖ψ0‖=1\| _0\|=1, • For each token x∈Σx∈ , a unitary matrix Wx∈U(N)W_x∈ U(N), • A set of measurement vectors |mk⟩k=0V−1\ m_k\_k=0^V-1 in ℂNC^N satisfying the resolution of identity ∑k|mk⟩⟨mk|=IN _k m_k m_k=I_N. Given an input sequence x1,x2,…,xTx_1,x_2,…,x_T, the model updates its state as |ψ(t)⟩=Wxt|ψ(t−1)⟩,t=1,…,T, ψ(t)=W_x_t ψ(t-1), t=1,…,T, (39) and produces output probabilities at each step via the Born rule: p(k∣ψ(t))=|⟨mk|ψ(t)⟩|2.p(k ψ(t))= | m_kψ(t) |^2. (40) The CUSM is a special case of the full architecture of Section 3, obtained by restricting gθg_θ to depend only on the token embedding (not on the state), setting H0=0H_0=0, and using Δt=1 t=1 in the Cayley discretization so that each token’s unitary is Wxt=Cayley(−iHint(xt)/2)W_x_t=Cayley(-iH_int(x_t)/2). The state-dependent nonlinearity discussed in Section 3.2 is absent from the CUSM; implementing it would strictly enlarge the class of representable functions. Definition 5.2 (Real Orthogonal Sequence Model (ROSM)). A ROSM of dimension d over an alphabet Σ and vocabulary =0,…,V−1V=\0,…,V-1\ consists of: • An initial state h0∈ℝdh_0 ^d with ‖h0‖=1\|h_0\|=1, • For each token x∈Σx∈ , an orthogonal matrix Qx∈O(d)Q_x∈ O(d), • Output weight vectors wkk=0V−1\w_k\_k=0^V-1 in ℝdR^d and a bias vector b∈ℝVb ^V. Given an input sequence x1,x2,…,xTx_1,x_2,…,x_T, the model updates its state as h(t)=Qxth(t−1),t=1,…,T,h(t)=Q_x_t\,h(t-1), t=1,…,T, (41) and produces output probabilities at each step via the affine-softmax readout: p(k∣h(t))=exp(wk⊤h(t)+bk)∑k′=0V−1exp(wk′⊤h(t)+bk′).p(k h(t))= \! (w_k h(t)+b_k ) _k =0^V-1 \! (w_k h(t)+b_k ). (42) Both model classes use norm-preserving transitions: the CUSM preserves ‖ψ‖=1\|ψ\|=1, and the ROSM preserves ‖h‖=1\|h\|=1. The distinction lies in the output mechanism. The CUSM uses the Born rule, which is a quadratic function of the complex state. The ROSM uses an affine projection followed by softmax, which depends on the state through the linear map h↦Wh+bh Wh+b. This difference in the algebraic degree of the output is the source of the separation. The affine-softmax readout in Definition 5.2 is the standard output mechanism used by recurrent neural networks, including the unitary and orthogonal RNNs of Arjovsky et al. arjovsky2016unitary and Wisdom et al. wisdom2016full . A real-valued model that replaced this readout with a quadratic function of the state (e.g., p(k)∝(wk⊤h)2p(k) (w_k h)^2) would itself perform an implicit polynomial lifting and would not be subject to the rank bound derived in Lemma 5.9. The separation theorem compares the two model classes as defined, not all possible combinations of dynamics and readouts. 5.2 The Disambiguation Task Family For each integer N≥2N≥ 2, we define a disambiguation task ND_N that requires the model to produce a specific probability distribution at the final position of a sequence, where the target distribution depends jointly on two tokens separated by an arbitrary number of filler tokens. Alphabet and sequence structure. The input alphabet is ΣN=a0,a1,…,aN−1∪b0,b1,…,bN−1∪σ _N=\a_0,a_1,…,a_N-1\∪\b_0,b_1,…,b_N-1\∪\σ\, containing N context tokens, N query tokens, and a single filler token. For a fixed sequence length T≥3T≥ 3, the input sequences have the form s=(ai,σ,…,σ⏟T−2,bj),i,j∈0,…,N−1.s=(a_i,\; σ,\;…,\;σ_T-2,\;b_j), i,j∈\0,…,N-1\. (43) There are N2N^2 such sequences, one for each pair (i,j)(i,j). Target distribution. The output vocabulary has V≥N2V≥ N^2 tokens. The target probability distribution at the final position (after processing bjb_j) is defined in terms of a specific Born-rule construction. Choose: 1. N unit vectors |ψ0⟩,…,|ψN−1⟩∈ℂN _0,…, _N-1 ^N (the context states), 2. N unitary matrices W0,…,WN−1∈U(N)W_0,…,W_N-1∈ U(N) (the query unitaries), 3. V measurement vectors |mk⟩k=0V−1\ m_k\_k=0^V-1 in ℂNC^N satisfying the resolution of identity ∑k|mk⟩⟨mk|=IN _k m_k m_k=I_N and the informationally complete condition (defined below). The target distribution for input pair (i,j)(i,j) is p∗(k∣i,j)=|⟨mk|Wjψi⟩|2=tr(MkWj|ψi⟩⟨ψi|Wj†),p^*(k i,j)= | m_kW_j _i |^2=tr\! (M_k\,W_j _i _iW_j ), (44) where Mk=|mk⟩⟨mk|M_k= m_k m_k. The context states, query unitaries, and measurement vectors are part of the task specification, not part of any model’s learned parameters. A model “solves” or “computes” the task ND_N if, for every input pair (i,j)(i,j) and every output token k, the model’s output probability at the final position equals p∗(k∣i,j)p^*(k i,j). Conditions on the task parameters. The context states and query unitaries must satisfy a spanning condition, and the measurement vectors must be informationally complete. We state these conditions precisely. Definition 5.3 (General-position unitaries). The context states |ψi⟩i=0N−1\ _i\_i=0^N-1 and query unitaries Wjj=0N−1\W_j\_j=0^N-1 are in general position if the N2N^2 density matrices ρij=Wj|ψi⟩⟨ψi|Wj†,i,j∈0,…,N−1, _ij=W_j _i _iW_j , i,j∈\0,…,N-1\, (45) are linearly independent as elements of the N2N^2-dimensional real vector space of N×N× N Hermitian matrices. Definition 5.4 (Informationally complete measurement). The measurement vectors |mk⟩k=0V−1\ m_k\_k=0^V-1 form an informationally complete measurement if the map ρ↦(tr(M0ρ),…,tr(MV−1ρ))ρ (tr(M_0ρ),\,…,\,tr(M_V-1ρ) ) is injective on the space of N×N× N Hermitian matrices. Equivalently, the V matrices Mk=|mk⟩⟨mk|\M_k= m_k m_k\ span the space of N×N× N Hermitian matrices as a real vector space. Informational completeness requires V≥N2V≥ N^2, since the space of Hermitian matrices has real dimension N2N^2. Such measurements exist for every N: for instance, the N2N^2 outer products of a set of N2N^2 vectors forming a weighted 2-design in ℂNC^N provide an informationally complete measurement renes2004symmetric . For the separation theorem, we fix any informationally complete measurement with V≥N2V≥ N^2, and choose context states and query unitaries in general position; the next subsection verifies that general-position configurations exist for every N. 5.3 Upper Bound: Complex Unitary Construction Proposition 5.5 (CUSM upper bound). For every N≥2N≥ 2 and every general-position configuration of context states and query unitaries, a CUSM of dimension N computes ND_N exactly. Proof. Set the initial state of the CUSM to |ψ0init⟩=|ψ0⟩ _0^init= _0 (the first context state). For each context token aia_i, define the transition unitary as Wai=ViW_a_i=V_i, where Vi∈U(N)V_i∈ U(N) is any unitary satisfying Vi|ψ0⟩=|ψi⟩V_i _0= _i. Such a unitary exists for every pair of unit vectors in ℂNC^N. For the filler token, set Wσ=INW_σ=I_N. For each query token bjb_j, set Wbj=WjW_b_j=W_j. Equip the model with the measurement vectors |mk⟩\ m_k\ from the task specification. After processing the input sequence (ai,σT−2,bj)(a_i,σ^T-2,b_j), the state is |ψ(T)⟩=Wj⋅INT−2⋅Vi⋅|ψ0⟩=Wj|ψi⟩. ψ(T)=W_j· I_N^T-2· V_i· _0=W_j _i. (46) The Born-rule output at position T is p(k∣ψ(T))=|⟨mk|Wjψi⟩|2=p∗(k∣i,j),p(k ψ(T))= | m_kW_j _i |^2=p^*(k i,j), (47) which matches the target distribution (44) for every pair (i,j)(i,j) and every output token k. ∎ Existence of general-position configurations. We verify that the general-position condition (Definition 5.3) is satisfiable for every N≥2N≥ 2. Lemma 5.6 (General position is generic). For every N≥2N≥ 2, the set of (|ψ0⟩,…,|ψN−1⟩,W0,…,WN−1)( _0,…, _N-1,W_0,…,W_N-1) that are not in general position has Lebesgue measure zero in the product space (S2N−1)N×U(N)N(S^2N-1)^N× U(N)^N. Proof. The N2N^2 density matrices ρij=Wj|ψi⟩⟨ψi|Wj† _ij=W_j _i _iW_j can be vectorized as elements of ℝN2R^N^2 using any fixed orthonormal basis E1,…,EN2\E_1,…,E_N^2\ for the real vector space of N×N× N Hermitian matrices (for instance, the generalized Gell-Mann matrices). Concretely, define the vectorization vec(ρ)∈ℝN2vec(ρ) ^N^2 by [vec(ρ)]α=tr(Eαρ)[vec(ρ)]_α=tr(E_α\,ρ). Stacking these vectors as rows of an N2×N2N^2× N^2 real matrix R, the general-position condition is equivalent to det(R)≠0 (R)≠ 0. The determinant det(R) (R) is a real-analytic function of the parameters (|ψi⟩,Wj)( _i,W_j). A real-analytic function that is not identically zero has a zero set of measure zero. It therefore suffices to exhibit a single configuration for which det(R)≠0 (R)≠ 0. For N=2N=2, we provide an explicit construction. Set W0=I2W_0=I_2 and W1=12(111−1),W_1= 1 2 pmatrix1&1\\ 1&-1 pmatrix, (48) and choose the context states |ψ0⟩=|0⟩=(10),|ψ1⟩=12(1i). _0= 0= pmatrix1\\ 0 pmatrix, _1= 1 2 pmatrix1\\ i pmatrix. (49) Note that |ψ1⟩ _1 is not orthogonal to |ψ0⟩ _0: ⟨ψ0|ψ1⟩=1/2≠0 _0 _1=1/ 2≠ 0. The four density matrices are ρ00 _00 =|0⟩⟨0|=(1000), = 0 0= pmatrix1&0\\ 0&0 pmatrix, (50) ρ10 _10 =|ψ1⟩⟨ψ1|=12(1−i1), = _1 _1= 12 pmatrix1&-i\\ i&1 pmatrix, (51) ρ01 _01 =W1|0⟩⟨0|W1†=12(1111), =W_1 0 0W_1 = 12 pmatrix1&1\\ 1&1 pmatrix, (52) ρ11 _11 =W1|ψ1⟩⟨ψ1|W1†=12(1i−i1). =W_1 _1 _1W_1 = 12 pmatrix1&i\\ -i&1 pmatrix. (53) Parametrize each 2×22× 2 Hermitian matrix (αu+ivu−ivγ) pmatrixα&u+iv\\ u-iv&γ pmatrix by the coordinate vector (α,u,v,γ)∈ℝ4(α,u,v,γ) ^4. In these coordinates, the four density matrices become the rows of R=(10001/20−1/21/21/21/201/21/201/21/2).R= pmatrix1&0&0&0\\ 1/2&0&-1/2&1/2\\ 1/2&1/2&0&1/2\\ 1/2&0&1/2&1/2 pmatrix. (54) A direct computation gives det(R)=−1/4≠0 (R)=-1/4≠ 0, confirming linear independence for this configuration. For general N>2N>2, the existence of a nonsingular configuration follows from the fact that the set of rank-one projectors |ψ⟩⟨ψ|:|ψ⟩∈ℂN,‖ψ‖=1\ ψ ψ: ψ ^N,\,\|ψ\|=1\ (the Veronese-type variety) spans the full N2N^2-dimensional space of Hermitian matrices. This spanning property can be verified directly: for any pair j<kj<k, the projectors |ej⟩⟨ej| e_j e_j generate the diagonal subspace; the projectors 12(|ej⟩+|ek⟩)(⟨ej|+⟨ek|) 12( e_j+ e_k)( e_j+ e_k) generate the real off-diagonal directions 12(|ej⟩⟨ek|+|ek⟩⟨ej|) 12( e_j e_k+ e_k e_j) (after subtracting the diagonal contribution); and the projectors 12(|ej⟩+i|ek⟩)(⟨ej|−i⟨ek|) 12( e_j+i e_k)( e_j-i e_k) generate the imaginary off-diagonal directions 12i(|ej⟩⟨ek|−|ek⟩⟨ej|) 12i( e_j e_k- e_k e_j). Together these produce N+2(N2)=N2N+2 N2=N^2 linearly independent Hermitian matrices, confirming that the Veronese-type variety spans ℝN2R^N^2. Since N2N^2 generic points on a spanning variety are linearly independent (the determinant of the corresponding N2×N2N^2× N^2 matrix is a nonzero polynomial in the parameters of the points), and the set of configurations for which this determinant vanishes has measure zero, the claim follows for all N. ∎ The role of non-orthogonal context states. The proof uses context states |ψi⟩ _i that are not mutually orthogonal. This is essential. If the context states form an orthonormal basis |i⟩i=0N−1\ i\_i=0^N-1, then for each fixed query unitary WjW_j, the N density matrices Wj|i⟩⟨i|Wj†i=0N−1\W_j i iW_j \_i=0^N-1 sum to the identity: ∑iWj|i⟩⟨i|Wj†=WjIWj†=I _iW_j i iW_j =W_jIW_j =I. This linear constraint reduces the number of independent density matrices per unitary from N to N−1N-1, and with N unitaries, the total number of independent density matrices is at most N(N−1)=N2−N(N-1)=N^2-N, which falls short of the N2N^2 needed for full span by N dimensions. Using non-orthogonal context states removes the constraint ∑iρij=I _i _ij=I, allowing each unitary to contribute N independent density matrices and enabling the full N2N^2-dimensional span. 5.4 Lower Bound for Real Orthogonal Models We now prove that any ROSM computing ND_N requires real dimension d≥N2−2d≥ N^2-2, under the assumption that the target log-probability matrix L∗L^* has full row rank N2N^2. The proof proceeds through three lemmas: (i) the algebraic structure of the Born-rule output (Lemma 5.7), a standard quadratic-lifting observation from quantum information theory nielsen2000quantum ; renes2004symmetric ; (i) the rank of the probability matrix P∗P^* (Lemma 5.8); and (i) the rank constraint imposed by the affine-softmax readout yang2018breaking (Lemma 5.9). The lower bound in Theorem 5.10 additionally assumes full row rank of the target log-probability matrix L∗L^*; we do not rely on rank preservation under entrywise logarithms, which does not hold in general. The Remark following Theorem 5.10 discusses this assumption and identifies a log-odds reformulation as a path toward an unconditional result. The new contribution is the combination of these ingredients to yield a formal conditional dimensional separation for sequence models (Theorem 5.10). Lemma 5.7 (Quadratic lifting). The Born-rule probability (44) depends on the state |ψ⟩ ψ only through the rank-one density matrix ρ=|ψ⟩⟨ψ|ρ= ψ ψ, and the dependence is linear. Specifically, fix an orthonormal basis E1,…,EN2\E_1,…,E_N^2\ for the N2N^2-dimensional real vector space of N×N× N Hermitian matrices (with respect to the inner product ⟨A,B⟩=tr(AB) A,B =tr(AB)), and define the vectorization vec(ρ)∈ℝN2vec(ρ) ^N^2 by [vec(ρ)]α=tr(Eαρ)[vec(ρ)]_α=tr(E_α\,ρ). Then p∗(k∣ψ)=tr(Mkρ)=vec(Mk)⊤vec(ρ),p^*(k ψ)=tr(M_k\,ρ)=vec(M_k) \,vec(ρ), (55) where vec(Mk)∈ℝN2vec(M_k) ^N^2 is the vectorization of Mk=|mk⟩⟨mk|M_k= m_k m_k in the same basis. Proof. The Born-rule probability is p∗(k∣ψ)=|⟨mk|ψ⟩|2=tr(Mk|ψ⟩⟨ψ|)=tr(Mkρ)p^*(k ψ)=| m_kψ|^2=tr(M_k ψ ψ)=tr(M_kρ). Since both MkM_k and ρ are Hermitian, tr(Mkρ)=∑αtr(EαMk)tr(Eαρ)=vec(Mk)⊤vec(ρ)tr(M_kρ)= _αtr(E_αM_k)\,tr(E_αρ)=vec(M_k) vec(ρ), where the second equality uses the orthonormality of Eα\E_α\. ∎ This lemma states that the Born-rule output is a quadratic function of the complex amplitudes cjc_j (since ρjk=cjck∗ _jk=c_jc_k^*), but a linear function of the vectorized density matrix vec(ρ)vec(ρ). The density matrix has N2N^2 real coordinates: N diagonal entries |cj|2|c_j|^2 and N(N−1)N(N-1) real parameters in the off-diagonal entries cjck∗c_jc_k^* (each complex off-diagonal entry contributes two real numbers, but Hermiticity halves the independent count). The Born rule, by computing a linear function of vec(ρ)vec(ρ), accesses all N2N^2 of these features, even though the state |ψ⟩ ψ itself has only 2N−12N-1 real degrees of freedom. The extra features arise from the cross terms cjck∗c_jc_k^*, which encode pairwise phase relationships between amplitudes. This quadratic lifting nielsen2000quantum is the mechanism through which the Born rule achieves O(N2)O(N^2) effective features from an N-dimensional complex state; the formal separation in Theorem 5.10 is a consequence of this lifting combined with the rank constraint on the softmax readout. Lemma 5.8 (Rank of the probability matrix). Let ND_N use context states and query unitaries in general position and an informationally complete measurement with V≥N2V≥ N^2. Define the N2×VN^2× V probability matrix P∗P^* by P(i,j),k∗=p∗(k∣i,j)=tr(Mkρij),P^*_(i,j),\,k=p^*(k i,j)=tr(M_k\, _ij), (56) where rows are indexed by pairs (i,j)∈0,…,N−12(i,j)∈\0,…,N-1\^2 and columns by output tokens k∈0,…,V−1k∈\0,…,V-1\. Then rank(P∗)=N2rank(P^*)=N^2 for generic choices of the task parameters. We do not derive any consequence for the entrywise logarithm L∗=logP∗L^*= P^*; entrywise logarithms do not preserve matrix rank in general, and the lower bound (Theorem 5.10) is stated conditionally on a full-rank assumption on L∗L^* rather than on any rank-preservation property of the logarithm. Proof. The proof proceeds in two steps. Step 1: The probability matrix has rank N2N^2. By Lemma 5.7, each row of P∗P^* is vec(ρij)⊤M^vec( _ij) M, where M^∈ℝN2×V M ^N^2× V has columns vec(Mk)vec(M_k). In matrix form, P∗=R^M^⊤P^*= R\, M , where R^∈ℝN2×N2 R ^N^2× N^2 has rows vec(ρij)⊤vec( _ij) . The general-position condition guarantees rank(R^)=N2rank( R)=N^2 (Definition 5.3), and informational completeness guarantees rank(M^)=N2rank( M)=N^2 (Definition 5.4). Therefore rank(P∗)=N2rank(P^*)=N^2. Step 2: All entries of P∗P^* are positive. Each entry P(i,j),k∗=|⟨mk|Wjψi⟩|2P^*_(i,j),\,k=| m_kW_j _i|^2. A squared modulus of an inner product of two nonzero vectors vanishes only when the vectors are orthogonal. For generic choices of |ψi⟩ _i, WjW_j, and |mk⟩ m_k, no state Wj|ψi⟩W_j _i is exactly orthogonal to any measurement vector |mk⟩ m_k, because the set of orthogonal pairs has measure zero in the product of unit spheres. Therefore P∗>0P^*>0 entrywise for generic parameters. ∎ Lemma 5.9 (Softmax rank bound). Let a ROSM of dimension d produce output probabilities p¯(k∣c) p(k c) for C input contexts via the affine-softmax readout (42), with hidden states hc∈ℝdh_c ^d, weight vectors wk∈ℝdw_k ^d, and bias b∈ℝVb ^V. Define the C×VC× V log-probability matrix L¯ L by L¯c,k=logp¯(k∣c) L_c,k= p(k c). Then rank(L¯)≤d+2.rank( L)≤ d+2. (57) Proof. The log-probability under the softmax readout is logp¯(k∣c)=wk⊤hc+bk−log∑k′exp(wk′⊤hc+bk′). p(k c)=w_k h_c+b_k- _k (w_k h_c+b_k ). (58) Define the logit matrix Z∈ℝC×VZ ^C× V by Zc,k=wk⊤hc+bkZ_c,k=w_k h_c+b_k, and the log-partition vector λ∈ℝCλ ^C by λc=log∑k′exp(Zc,k′) _c= _k (Z_c,k ). Then L¯=Z−λ 1V⊤ L=Z-λ\,1_V , where V∈ℝV1_V ^V is the all-ones vector. The logit matrix factors as Z=HWout⊤+Cb⊤Z=H\,W_out +1_C\,b , where H∈ℝC×dH ^C× d has rows hc⊤h_c and Wout∈ℝV×dW_out ^V× d has rows wk⊤w_k . By subadditivity of matrix rank, rank(Z)≤rank(HWout⊤)+rank(Cb⊤)≤d+1rank(Z) (HW_out )+rank(1_Cb )≤ d+1. Since λV⊤ 1_V is a rank-one matrix, rank(L¯)=rank(Z−λV⊤)≤rank(Z)+rank(λV⊤)≤(d+1)+1=d+2.∎rank( L)=rank(Z- 1_V ) (Z)+rank( 1_V )≤(d+1)+1=d+2. (59) This rank bound is a consequence of the same structural property identified by Yang et al. yang2018breaking in their analysis of the softmax bottleneck: because the logits depend on the hidden state through a linear projection, the resulting log-probability matrix is constrained to a low-dimensional affine subspace of ℝC×VR^C× V. The bound is specific to the affine-softmax readout. A readout that computes a nonlinear function of h before applying softmax (such as a two-layer MLP, or the mixture-of-softmax readout proposed in yang2018breaking ) would produce a log-probability matrix whose rank is not bounded by d+2d+2, and the lower bound in Theorem 5.10 would not apply to such a model. The separation theorem therefore compares two specific model families, not all possible architectures. Theorem 5.10 (Expressivity separation (conditional)). For every N≥2N≥ 2, there exists a disambiguation task ND_N with vocabulary size V=N2V=N^2 satisfying the general-position and informational-completeness conditions such that: 1. A CUSM of complex dimension N computes ND_N exactly. 2. Suppose additionally that the target log-probability matrix L∗L^*, with entries L(i,j),k∗=logp∗(k∣i,j)L^*_(i,j),k= p^*(k i,j), has full row rank N2N^2. Then any ROSM of real dimension d with affine-softmax readout that computes ND_N satisfies d≥N2−2d≥ N^2-2. Proof. Part (1) is Proposition 5.5. For part (2), suppose a ROSM of dimension d computes ND_N, meaning its output probabilities p¯(k∣i,j)=p∗(k∣i,j) p(k i,j)=p^*(k i,j) for all i,j,ki,j,k. Then the log-probability matrices coincide: L¯=L∗ L=L^*. By assumption, rank(L∗)=N2rank(L^*)=N^2. By Lemma 5.9, rank(L¯)≤d+2rank( L)≤ d+2. Combining these: N2=rank(L∗)=rank(L¯)≤d+2,N^2=rank(L^*)=rank( L)≤ d+2, (60) which gives d≥N2−2d≥ N^2-2. ∎ Remark 5.11. The full-rank condition on L∗L^* in part (2) of Theorem 5.10 is the only point in the lower-bound argument where a property beyond the general-position and informational-completeness conditions on the task parameters is required. We do not rely on the claim that entrywise logarithms preserve matrix rank: this does not hold in general, and rank(logP)≠rank(P)rank( P) (P) is possible even for 2×22× 2 matrices. The condition rank(L∗)=N2rank(L^*)=N^2 is a genericity-type assumption on the task-construction parameters; a complete analytic characterization of when it holds is left for future work. One promising route toward removing this assumption is to reformulate the rank argument in terms of relative log-odds. For a fixed reference token k0k_0, define ℓ(i,j),k∗=logp∗(k∣i,j)p∗(k0∣i,j) ^*_(i,j),k= p^*(k i,j)p^*(k_0 i,j) for k≠k0k≠ k_0. Under the affine-softmax readout (Definition 5.2), the relative logits of the ROSM are exactly affine in the hidden state: ℓ¯c,k=(wk−wk0)⊤hc+(bk−bk0), _c,k=(w_k-w_k_0) h_c+(b_k-b_k_0), (61) so the N2×(V−1)N^2×(V-1) matrix of relative logits has rank at most d, providing a cleaner rank constraint without any assertion about entrywise logarithms. Whether a rank-N2N^2 condition on ℓ∗ ^* follows from the task structure alone remains an open question that we leave for future work. The separation in Theorem 5.10 is between the complex dimension of the CUSM and the real dimension of the ROSM. A CUSM of complex dimension N uses 2N2N real numbers to represent its state (the real and imaginary parts of N complex amplitudes), while the ROSM requires at least N2−2N^2-2 real numbers. The ratio (N2−2)/(2N)(N^2-2)/(2N) grows as N/2N/2 for large N, so the real orthogonal model needs a state space that is a factor of Θ(N) (N) larger than the real representation of the complex state. In terms of the complex dimension N alone, the gap is quadratic: a complex model of dimension N achieves what a real model cannot achieve below dimension Ω(N2) (N^2). 5.5 The Source of the Gap The separation proved in Theorem 5.10 has a precise algebraic origin: it arises from the interaction between the quadratic Born-rule readout and the linear softmax readout, not from any intrinsic limitation of real-valued dynamics. Tracing the proof identifies two specific mechanisms. The quadratic lifting. A complex amplitude cj=rjeiθjc_j=r_je^i _j carries two real numbers per coordinate, but the quadratic gap does not arise from this factor-of-two storage advantage (which accounts for only a constant factor, not a quadratic one). The gap arises from the Born rule’s quadratic dependence on the state. The Born-rule probability for output k is p(k∣ψ)=|∑j=0N−1[mk]j∗cj|2=∑j=0N−1|[mk]j|2|cj|2+2∑j<j′Re([mk]j∗[mk]j′cj∗cj′).p(k ψ)= | _j=0^N-1[m_k]_j^*c_j |^2= _j=0^N-1|[m_k]_j|^2|c_j|^2+2 _j<j Re\! ([m_k]_j^*[m_k]_j \,c_j^*c_j ). (62) The first sum involves N terms that depend only on magnitudes |cj||c_j|. The second sum involves (N2) N2 cross terms, each depending on a relative phase θj′−θj _j - _j through cj∗cj′=rjrj′ei(θj′−θj)c_j^*c_j =r_jr_j e^i( _j - _j). The density matrix ρ=|ψ⟩⟨ψ|ρ= ψ ψ collects all N2N^2 of these features (diagonal magnitudes and off-diagonal phase products), and the Born rule reads them off via the linear map ρ↦tr(Mkρ)ρ (M_kρ) (Lemma 5.7). This is a Veronese-type lifting: the map |ψ⟩↦|ψ⟩⟨ψ| ψ ψ ψ sends the N-dimensional complex unit sphere into the N2N^2-dimensional space of Hermitian matrices, and the Born rule then computes a linear function in the lifted space nielsen2000quantum . An analogous lifting is exploited by Born machines in generative modeling cheng2018information ; liu2018differentiable , where this quadratic structure enables expressive probability estimation without explicit normalization. The linear bottleneck of softmax. The ROSM’s softmax readout computes logp(k∣h)=wk⊤h+bk−logZ p(k h)=w_k h+b_k- Z, which depends on the d-dimensional state h through the linear projection wk⊤hw_k h. The full log-probability vector lies in a (d+2)(d+2)-dimensional subspace (Lemma 5.9). When the Born-rule target requires N2N^2 independent features, this subspace must have dimension at least N2N^2, forcing d≥N2−2d≥ N^2-2. The ROSM cannot synthesize pairwise features from individual state components, because the softmax readout is linear in h yang2018breaking . A ROSM with a nonlinear readout (e.g., a polynomial or MLP-based readout) would not be subject to this linear bottleneck, and the lower bound of Theorem 5.10 would not apply to it. What the theorem does not claim. The separation applies specifically to the comparison between the Born-rule readout (quadratic, complex) and the affine-softmax readout (linear, real). It does not establish a fundamental limitation of real-valued dynamics independent of the readout mechanism. A real-valued model with a quadratic readout p(k)∝(wk⊤h)2p(k) (w_k h)^2 with appropriate normalization would perform its own Veronese-type lifting and could potentially match the Born rule’s feature count. Such readouts are not standard in the architectures surveyed in Section 2, but they are not precluded by any fundamental constraint. The theorem does not require the task ND_N to correspond to any specific linguistic phenomenon. The task family uses abstract context and query tokens, not natural-language words. Whether the algebraic structure exploited by ND_N (pairwise token interactions encoded in relative phases) is relevant to natural-language processing is an empirical question. The task serves to demonstrate a provable gap between two specific architectures, not to model a particular linguistic phenomenon. The theorem is also purely representational: it characterizes the minimum state dimension required to represent a given family of input-output mappings, not whether gradient-based optimization can find the parameters that realize this mapping. A model with sufficient representational capacity may fail to learn the target function if the loss landscape is unfavorable. The experimental protocols in Section 6.2 are designed to test whether the representational advantage identified here manifests in gradient-trained models. 6 Discussion The preceding sections developed three theoretical contributions: an architecture whose state evolves unitarily in ℂNC^N under a learned Hamiltonian and whose output follows the Born rule (Section 3), a continuity equation that decomposes each state update into antisymmetric pairwise probability currents (Section 4), and a separation theorem establishing that a complex unitary model of dimension N represents a family of disambiguation tasks that requires dimension Ω(N2) (N^2) in any real orthogonal model with affine-softmax readout (Section 5). This section examines four aspects of the framework that the theoretical development leaves open: the computational cost of the architecture, the concrete predictions that the theory makes and the experimental protocols that would test them, the boundaries of what the present analysis does and does not establish, and the theoretical directions that extend beyond the results proved here. Because the present paper is entirely theoretical, the experimental protocols described in Section 6.2 are intended as a precise specification of the empirical program needed to validate or refute the theory’s claims, not as a report of completed experiments. 6.1 Computational Complexity The architecture defined in Section 3 involves four computational stages at each time step: generating the interaction Hamiltonian from the current token and state, conjugating into the interaction picture, solving the Cayley linear system, and computing Born-rule output probabilities. We analyze the cost of each stage and the resulting per-step and per-sequence complexity. Per-step forward cost. At step t, the network gθg_θ receives the concatenation of the token embedding Embed(xt)∈ℝdEmbed(x_t) ^d and the real-imaginary representation of the current state [Re(c(t)),Im(c(t))]∈ℝ2N[Re(c(t)),\,Im(c(t))] ^2N, and produces the matrix Φ(t)∈ℂN×r (t) ^N× r and the vector δ(t)∈ℝNδ(t) ^N. The output has 2Nr+N2Nr+N real components. If gθg_θ is a feedforward network with LgL_g layers and maximum hidden width h, its evaluation cost is Cgθ=O(Lg⋅h⋅max(d+2N,h, 2Nr+N)),C_g_θ=O(L_g· h· (d+2N,\,h,\,2Nr+N)), which for d, N, and h of comparable magnitude simplifies to O(Lgh2)O(L_gh^2). The interaction-picture conjugation (equation (13)) replaces each entry Φja(t) _ja(t) with Φ~ja(t)=eiλjtΦja(t) _ja(t)=e^i _jt\, _ja(t), which costs O(Nr)O(Nr) complex multiplications. The resulting interaction-picture Hamiltonian has the form Hint,I(t)=Φ~(t)Φ~(t)†+diag(δ(t)),H_int,I(t)= (t)\, (t) +diag(δ(t)), (63) which inherits the low-rank-plus-diagonal structure of Hint(t)H_int(t) because the conjugation by the diagonal matrix eiH0te^iH_0t preserves diagonal structure and maps ΦΦ† to Φ~Φ~† . The Cayley update (15) requires solving the linear system (D+iΔt2Φ~Φ~†)⏟A|ψt+1⟩=(D′−iΔt2Φ~Φ~†)|ψt⟩⏟b, (D+ i t2\, \, )_A\, _t+1= (D - i t2\, \, ) _t_b, (64) where D=I+iΔt2diag(δ)D=I+ i t2diag(δ) and D′=I−iΔt2diag(δ)D =I- i t2diag(δ) are both diagonal. Each diagonal entry of D has the form 1+iΔt2δj1+ i t2 _j, with modulus 1+(Δtδj/2)2>0 1+( t\, _j/2)^2>0, so D is invertible for all real δj _j and all Δt>0 t>0. The right-hand side b is computed by first forming v=Φ~†|ψt⟩∈ℂrv= _t ^r at cost O(Nr)O(Nr), then Φ~v∈ℂN \,v ^N at cost O(Nr)O(Nr), and finally b=D′|ψt⟩−iΔt2Φ~vb=D _t- i t2\, \,v at cost O(N)O(N). The total cost for the right-hand side is O(Nr)O(Nr). The system A|ψt+1⟩=bA _t+1=b is solved via the Woodbury matrix identity. Because A=D+iΔt2Φ~Φ~†A=D+ i t2 is a rank-r update of the invertible diagonal matrix D, the solution is |ψt+1⟩=D−1b−iΔt2D−1Φ~(Ir+iΔt2Φ~†D−1Φ~)−1Φ~†D−1b. _t+1=D^-1b- i t2\,D^-1 (I_r+ i t2\, D^-1 )^-1 D^-1b. (65) The computation proceeds as follows: (i) compute y=D−1by=D^-1b at cost O(N)O(N); (i) compute z=Φ~†y∈ℂrz= y ^r at cost O(Nr)O(Nr); (i) form P=D−1Φ~∈ℂN×rP=D^-1 ^N× r at cost O(Nr)O(Nr); (iv) form the r×r× r Gram matrix G=Ir+iΔt2Φ~†PG=I_r+ i t2\, P at cost O(Nr2)O(Nr^2); (v) solve Gw=zGw=z for w∈ℂrw ^r at cost O(r3)O(r^3); (vi) compute |ψt+1⟩=y−iΔt2Pw _t+1=y- i t2\,Pw at cost O(Nr)O(Nr). The dominant cost is step (iv), giving a total Cayley solve cost of O(Nr2+r3)O(Nr^2+r^3). Since the design intent is r≪Nr N (Section 3.2), this simplifies to O(Nr2)O(Nr^2). A caveat on numerical stability is warranted. The Woodbury solve requires inverting the r×r× r matrix G=Ir+iΔt2Φ~†D−1Φ~G=I_r+ i t2\, D^-1 . The conditioning of G depends on the eigenvalues of Φ~†D−1Φ~ D^-1 , which in turn depend on the interaction Hamiltonian output by gθg_θ. For large step sizes Δt t or large Hamiltonian eigenvalues, G may become poorly conditioned, degrading the accuracy of the solve. No general bound on the condition number of G is available without additional constraints on the output range of gθg_θ. In practice, constraining the spectral norm of Φ(t) (t) through regularization of gθg_θ, or monitoring the condition number of G during training, would be necessary to ensure reliable numerics hairer2006geometric ; golub2013matrix . We flag this as a point requiring empirical investigation in any implementation. The Born-rule output (18) requires computing M†|ψ(t)⟩∈ℂVM ψ(t) ^V (where M∈ℂN×VM ^N× V is the measurement matrix) and then taking entry-wise squared moduli. The matrix-vector product costs O(NV)O(NV), and the subsequent squaring costs O(V)O(V). The total Born-rule cost is therefore O(NV)O(NV). Combining all stages, the per-step forward cost is O(Cgθ+Nr2+NV).O\! (C_g_θ+Nr^2+NV ). (66) For a vocabulary of size V in the tens of thousands and a latent dimension N in the hundreds with r in the low tens, the O(NV)O(NV) Born-rule cost dominates. This cost is the same order as the output-projection cost O(dV)O(dV) incurred by any architecture that computes a dense projection from a hidden state of dimension d to a vocabulary of size V, so the Born-rule decoding does not introduce additional asymptotic overhead relative to the standard softmax readout. The row-orthonormality constraint on the measurement matrix M, which is required for the Born rule to produce valid probabilities (equation (19)), is enforced via QR decomposition of M†∈ℂV×NM ^V× N at each forward pass. This projects M onto the Stiefel manifold absil2008optimization ; townsend2016pymanopt before each evaluation, but the stored parameter M is not itself constrained to the manifold; only its QR-decomposed version is used in the forward pass. The gradient of the loss with respect to the raw parameter flows back through the QR decomposition, which is a non-trivial operation on the Stiefel manifold whose correct handling requires that the automatic differentiation framework propagate gradients through the QR factorization accurately wisdom2016full ; townsend2016pymanopt . An alternative that avoids this issue entirely is to parameterize M directly on the Stiefel manifold using Riemannian gradient descent or a retraction-based optimizer absil2008optimization ; we do not prescribe a specific approach here but note that the choice of optimization strategy for M is consequential and requires careful treatment in any implementation. Backward pass. Backpropagation through the Cayley solve at each step requires solving one additional linear system with coefficient matrix A†=(I+iΔt2Hint,I)†=I−iΔt2Hint,IA =(I+ i t2\,H_int,I) =I- i t2\,H_int,I. Because Hint,IH_int,I is Hermitian, this adjoint matrix has the form D∗−iΔt2Φ~Φ~†D^*- i t2\, , where D∗=I−iΔt2diag(δ)D^*=I- i t2diag(δ) is again an invertible diagonal matrix. The Woodbury identity applies with the same structure and cost O(Nr2)O(Nr^2). The backward pass through the Born-rule output and through gθg_θ adds costs of the same order as the forward pass through each respective stage. The total per-step backward cost is therefore also O(Cgθ+Nr2+NV)O(C_g_θ+Nr^2+NV), and the per-step cost scales identically in both directions. It is important to distinguish two separate stability properties of the backward pass. The unitary structure of each Cayley step guarantees that the gradient of the loss with respect to the hidden state at time t propagates through T−tT-t steps without change in norm: since each step applies a unitary map to the state, the adjoint of that map is also unitary, and the chain rule preserves the gradient norm exactly vorontsov2017orthogonality . However, this guarantee does not extend to the gradients with respect to the parameters of gθg_θ. At each step t, the parameter gradient passes through gθg_θ, which is an unconstrained neural network without norm-preservation guarantees. Gradient vanishing and exploding can therefore occur inside gθg_θ at each time step, independently of the unitary structure of the state update. Standard mitigation strategies for deep networks, such as careful initialization, batch normalization within gθg_θ, or gradient clipping applied to the parameter gradients (but not to the state gradients), remain necessary. Memory. The recurrent state |ψ(t)⟩∈ℂN ψ(t) ^N requires O(N)O(N) storage, independent of the sequence length T. The measurement matrix M∈ℂN×VM ^N× V requires O(NV)O(NV) storage. The free frequencies λj\ _j\ require O(N)O(N) storage. The parameters of gθg_θ require O(|θ|)O(|θ|) storage, where |θ||θ| is the parameter count of the network. The QR decomposition of M†M costs O(VN2)O(VN^2) and depends only on the parameter matrix M, not on the input sequence; it is therefore performed once per parameter update and its cost is amortized over the T steps of the sequence to O(VN2/T)O(VN^2/T) per step. Comparison with existing architectures. The per-step state-update cost of O(Nr2)O(Nr^2) compares with O(N2)O(N^2) for the full-capacity unitary RNN of Wisdom et al. wisdom2016full (which applies a dense N×N× N unitary matrix at each step) and O(NlogN)O(N N) for the factored uRNN of Arjovsky et al. arjovsky2016unitary . The reduction from O(N2)O(N^2) to O(Nr2)O(Nr^2) is the cost of restricting the interaction Hamiltonian to rank r rather than using a full unitary parameterization. This restriction limits the number of independent coupling patterns per step to r (Section 3.2), which is a representational tradeoff: fewer coupling channels per step in exchange for lower per-step cost. Whether this tradeoff is favorable depends on whether the tasks encountered in practice require more than r simultaneous coupling patterns at a single time step, a question that the theory developed here does not answer and that requires empirical investigation. The O(N)O(N) memory cost for the recurrent state is shared by all recurrent architectures, including standard RNNs, LSTMs hochreiter1997long , unitary RNNs, and state-space models gu2023mamba , and contrasts with the O(TdL)O(TdL) key-value cache of a Transformer with L layers, hidden dimension d, and sequence length T vaswani2017attention . The distinction is quantitatively significant for long sequences: at T=104T=10^4, d=512d=512, and L=12L=12, the Transformer cache requires approximately 1.2×1081.2× 10^8 floating-point numbers, while the recurrent state requires 2N2N (accounting for real and imaginary parts). 6.2 Testable Predictions and Experimental Protocols The theoretical analysis yields five predictions, each traceable to a specific result or construction in the preceding sections. For each prediction, we state the claim, identify its theoretical basis, and describe an experimental protocol that would confirm or refute it. The predictions are ordered from those most directly implied by the formal results to those that are motivated by the architecture’s structure but not guaranteed by any theorem. Prediction 1: Quadratic dimensional scaling on synthetic disambiguation tasks. Theorem 5.10 asserts that for the task family ND_N (Section 5.2), a CUSM of complex dimension N achieves the optimal cross-entropy loss while any ROSM requires real dimension at least N2−2N^2-2 (under the full-rank condition on L∗L^*). This is a statement about representational capacity, not about learnability, so the prediction is that gradient-trained models exhibit a dimensional threshold consistent with the theorem’s bounds. The protocol is as follows. For each N∈4,8,16,32N∈\4,8,16,32\, construct an instance of ND_N by selecting context states and query unitaries in general position (Lemma 5.6 guarantees that random choices satisfy this condition with probability one) and an informationally complete measurement with V=N2V=N^2. Generate a training set by enumerating all N2N^2 input pairs (i,j)(i,j) and computing the target distributions from equation (44). The constraint enforcement strategies for both model classes must be specified: for the CUSM, the measurement matrix M is projected onto the Stiefel manifold via QR decomposition at each forward pass, and the unitary matrices WxW_x are parameterized via the Cayley transform of learnable skew-Hermitian matrices; for the ROSM, the orthogonal matrices QxQ_x are parameterized analogously via the matrix exponential of learnable skew-symmetric matrices. Both models are trained with the Adam optimizer kingma2015adam using a fixed learning rate of 10−310^-3, a cosine decay schedule, and at least five independent random seeds; the reported losses are means and standard deviations across seeds. Train CUSMs of dimensions N′∈N/2,N, 2NN ∈\N/2,\,N,\,2N\ and ROSMs of dimensions d∈N, 2N,N2/2,N2d∈\N,\,2N,\,N^2/2,\,N^2\ to minimize cross-entropy against the target distributions. For each trained model, report the gap between the achieved loss and the theoretical minimum ℒ∗=−1N2∑i,j∑kp∗(k∣i,j)logp∗(k∣i,j)L^*=- 1N^2 _i,j _kp^*(k i,j) p^*(k i,j). A loss gap is declared zero if it falls below 10−310^-3 nats. The theorem predicts that the CUSM loss gap reaches zero at N′=N =N, while the ROSM loss gap remains positive for all d<N2−2d<N^2-2 and reaches zero only at d≥N2−2d≥ N^2-2. If the ROSM achieves a zero loss gap at a dimension substantially below N2−2N^2-2, this would indicate either a failure of the general-position condition for the chosen task instance or a gap in the theorem’s lower bound argument. Prediction 2: Born-rule readout outperforms softmax readout on a shared complex state. The separation in Theorem 5.10 arises from the readout mechanism rather than from the dynamics alone. Lemma 5.9 shows that the affine-softmax readout constrains the log-probability matrix to rank at most d+2d+2 yang2018breaking , while the Born rule accesses N2N^2 features of the state through the quadratic lifting of Lemma 5.7. This predicts that replacing the Born-rule readout with a softmax readout applied to the same unitary state should degrade prediction quality, and that the degradation should grow with N. The protocol uses paired model variants. For each N∈64,128,256,512N∈\64,128,256,512\, train two models on the same language corpus with identical unitary dynamics (same H0H_0, same gθg_θ, same initial state). Model A uses Born-rule decoding as defined in Section 3.5. Model B decodes via p(k∣t)=softmax(Wo[Re(c(t));Im(c(t))]+b),p(k t)=softmax(W_o\,[Re(c(t));\,Im(c(t))]+b), where Wo∈ℝV×2NW_o ^V× 2N and b∈ℝVb ^V are learnable parameters. Both models have the same recurrent state dimension and differ only in the output layer. Compare validation perplexity as a function of N. The Born-rule model should achieve lower perplexity, and the perplexity ratio should increase with N, because the number of features accessible to the Born rule grows as N2N^2 while the number accessible to the softmax readout grows as 2N2N. We note that this prediction tests the readout mechanism in isolation; it does not compare against Transformer architectures or other strong baselines with different dynamics, and such a comparison would require controlling for total parameter count and computational cost, which is left to future work. Prediction 3: Probability current magnitude peaks at disambiguating tokens. Section 4.3 established that all inter-dimensional probability currents are driven exclusively by the interaction Hamiltonian Hint(t)H_int(t) (equation (33)), which is generated from the current token. When a token does not alter the model’s interpretation of the sequence, the interaction Hamiltonian can remain close to zero and the currents will be small. When a token resolves an ambiguity, the interaction Hamiltonian must redirect amplitude between competing interpretations, producing large currents. This predicts that the total current magnitude ‖J(t)‖=∑j<k|Jj←k(t)|\|J(t)\|= _j<k|J_j← k(t)| (67) is systematically larger at tokens that resolve semantic ambiguity than at tokens that do not. The protocol requires a trained CUSM and a test corpus annotated for lexical ambiguity. We propose using SemCor 3.0 miller1993semantic , which provides sense annotations for running text from the Brown corpus, as the primary evaluation corpus. The annotation scheme distinguishes tokens whose sense is contextually determined (and therefore potentially disambiguating) from function words and unambiguous content words. At each token position t in the test corpus, compute ‖J(t)‖\|J(t)\| using the midpoint current formula (37), which provides an exact decomposition of the discrete probability change without the O(Δt2)O( t^2) approximation error of the continuous-time currents. Partition the token positions into two groups: those annotated as carrying disambiguating sense information and all remaining positions. Compare the distributions of ‖J(t)‖\|J(t)\| between the two groups using a Wilcoxon rank-sum test with a Bonferroni correction for multiple comparisons across layers or model sizes. The prediction is that the disambiguating group has significantly higher total current (p<0.01p<0.01 after correction). A negative result would indicate either that the trained model does not use the current mechanism for disambiguation or that the annotation scheme does not capture the ambiguities relevant to the model’s internal dynamics. Prediction 4: Learned frequency spectrum correlates with linguistic timescales. Section 3.2 argued that the free frequencies λj\ _j\ provide the model with a bank of oscillators at different rates, and Section 3.3 showed that the interaction picture creates a frequency-selective coupling mechanism that favors amplitude exchange between dimensions with similar natural frequencies. Together, these properties provide the architectural capacity for the model to allocate different latent dimensions to features that vary at different temporal rates. Whether the optimization discovers this allocation is not guaranteed by any theorem in this paper; it is a hypothesis about the inductive bias of the architecture. The protocol is as follows. Train a CUSM on a standard language corpus. For each latent dimension j, compute the temporal autocorrelation function of the occupation probability pj(t)=|cj(t)|2p_j(t)=|c_j(t)|^2 on held-out test sequences: Rj(τ)=1T−τ∑t=1T−τ(pj(t)−p¯j)(pj(t+τ)−p¯j),R_j(τ)= 1T-τ _t=1^T-τ (p_j(t)- p_j ) (p_j(t+τ)- p_j ), (68) where p¯j=1T∑tpj(t) p_j= 1T _tp_j(t) is the time-averaged occupation. Define the decorrelation time τj∗ _j^* as the smallest τ for which Rj(τ)<Rj(0)/eR_j(τ)<R_j(0)/e. Plot τj∗ _j^* against |λj|| _j|. The hypothesis predicts a negative correlation: dimensions with large |λj|| _j| (fast oscillators) should track rapidly changing features and therefore have short decorrelation times, while dimensions with small |λj|| _j| (slow oscillators) should track gradually varying features and have long decorrelation times. A null result (no correlation between τj∗ _j^* and |λj|| _j|) would indicate that the optimization does not exploit the timescale structure afforded by the free Hamiltonian, and would motivate investigating alternative initialization or regularization strategies for the frequency spectrum. Prediction 5: Phase interference contributes to output quality. Section 5.5 decomposed the Born-rule probability into N diagonal terms that depend only on magnitudes |cj|2|c_j|^2 and (N2) N2 cross terms that depend on relative phases θj′−θj _j - _j. If the trained model encodes useful information in the relative phases of its state (as the separation theorem requires for the disambiguation tasks), then discarding the cross terms should degrade prediction quality. Define the diagonal-only readout p~(k∣ψ)=∑j=0N−1|[mk]j|2|cj|2∑k′∑j=0N−1|[mk′]j|2|cj|2, p(k ψ)= _j=0^N-1|[m_k]_j|^2\,|c_j|^2 _k _j=0^N-1|[m_k ]_j|^2\,|c_j|^2, (69) which retains only the magnitude-dependent terms and renormalizes. The prediction is that the full Born rule (18) achieves lower perplexity than the diagonal-only readout (69) on held-out data, and that the gap quantifies the contribution of phase interference to the model’s predictions. The protocol uses a single trained CUSM evaluated under two readout rules. At test time, compute the per-token log-probability under both the Born rule and the diagonal-only readout, using the same state trajectory (no retraining). Report the perplexity difference. Because the diagonal-only readout accesses N features (the squared magnitudes) while the full Born rule accesses N2N^2 features (including all cross terms), the gap should be substantial when the model has learned to encode pairwise relational information in phase. 6.3 Limitations of the Present Analysis The theoretical results established in Sections 3–5 are subject to several scope restrictions that we state explicitly, including some that affect the interpretation of the central results. The separation theorem applies to a simplified model class. The most significant structural limitation of the present work is that the CUSM used in the separation theorem (Definition 5.1) assigns a fixed, state-independent unitary to each token. The full architecture of Section 3 uses state-dependent Hamiltonians (Section 3.2), in which the network gθg_θ receives the current state as input and produces a Hamiltonian that depends on the full history of previous tokens. The full model therefore belongs to a strictly richer expressivity class than the CUSM. This means the separation theorem, which was motivated throughout Sections 1–3 as establishing an advantage for the full architecture, actually establishes an advantage only for the simplified state-independent special case. Any task solvable by the CUSM is also solvable by the full model, so the Ω(N2) (N^2) lower bound on real orthogonal models carries over, but the theorem does not characterize the additional expressivity gained by state-dependence. Proving a separation for the full state-dependent architecture is an open problem that requires different proof techniques, since the state-dependence breaks the algebraic structure exploited in Lemmas 5.7 and 5.8. The separation is between output mechanisms, not between complex and real dynamics. Theorem 5.10 compares the Born-rule readout with the affine-softmax readout. The lower bound in Lemma 5.9 holds for the specific form of the softmax readout defined in Definition 5.2: a linear projection of the state followed by softmax normalization yang2018breaking . A real-valued model equipped with a quadratic readout of the form p(k)∝(wk⊤h)2p(k) (w_k h)^2, or with a mixture-of-softmax readout yang2018breaking , would itself perform a lifting from ℝdR^d to a higher-dimensional feature space and could substantially narrow the gap. Such readouts are not standard in the architectures surveyed in Section 2, but they are not precluded by any fundamental constraint. The separation should therefore be understood as a statement about the advantage of Born-rule decoding over affine-softmax decoding, not as a claim that complex-valued dynamics are inherently superior to real-valued dynamics for all output mechanisms. This qualification applies wherever the paper describes the separation as being between “complex” and “real” models without further specification. The lower bound is conditional on a full-rank assumption. The lower bound in part (2) of Theorem 5.10 is conditioned on the assumption that the target log-probability matrix L∗L^* has full row rank N2N^2. This condition is stated as an explicit hypothesis in the theorem rather than derived from the task-construction parameters, and a complete analytic characterization of when it holds is deferred to future work. We do not claim that the full rank of the probability matrix P∗P^* (established by Lemma 5.8) implies the full rank of L∗=logP∗L^*= P^*: entrywise logarithms do not preserve matrix rank in general, and this step is not valid as a general mathematical statement (rank can increase or decrease under entrywise logarithm even in 2×22× 2 examples). The Remark following Theorem 5.10 identifies a log-odds reformulation that would yield a cleaner, potentially unconditional lower bound; developing this into a complete proof is left for future work. Additionally, Lemma 5.6 provides a complete explicit construction only for N=2N=2; for N>2N>2, the proof invokes a spanning property of the Veronese-type variety that is stated but not fully executed. The interaction picture does not reduce computation in the discrete-step regime. Section 3.3 motivates the interaction picture by arguing that it allows the integrator to “resolve only the timescales introduced by the input-driven interaction, not the potentially much faster free oscillations.” This argument is valid for adaptive continuous-time integrators that adjust their step size to resolve fast oscillations. In the discrete-step regime used by the model, the Cayley transform is applied once per token with a fixed step size Δt t, regardless of the spectrum of H0H_0. The interaction-picture Hamiltonian has off-diagonal entries [Hint,I(t)]jk=[Hint(t)]jkei(λj−λk)t[H_int,I(t)]_jk=[H_int(t)]_jk\,e^i( _j- _k)t, which oscillate rapidly when |λj−λk|| _j- _k| is large. The Cayley transform does not adaptively resolve these oscillations; it applies a single midpoint step regardless. The claimed computational advantage of the interaction picture, in the sense of reduced integrator burden, therefore does not materialize in the discrete-step regime. The interaction picture remains a useful change of variables that separates the free and interaction dynamics conceptually, and the frequency-selective coupling argument of Section 3.3 applies in an averaged sense over many steps, but the computational motivation as stated does not hold for the model as implemented. This does not affect the correctness of the architecture, only the justification for one of its components. The quantum cognition motivation is a hypothesis, not an inference. Sections 1 and 2.6 cite Busemeyer and Bruza busemeyer2012quantum and Pothos and Busemeyer pothos2013quantum to motivate the use of complex amplitudes and interference in a language model. The cited works document violations of classical probability in human judgment experiments, and argue that these violations are naturally explained by quantum probability. However, the inference that a language model would benefit from the same mathematical structure does not follow directly. Quantum cognition models describe experimental data from specific psychological paradigms; they are not general models of linguistic computation. Whether the statistical structure of text corpora exhibits the kind of interference effects documented in those experiments is an empirical question that requires corpus-level evidence and has not been established. The connection to quantum cognition should be read as a motivating analogy and a hypothesis to be empirically tested, not as a logical justification for the architecture. The architecture’s formal properties—norm preservation, quadratic output, conserved probability currents—stand on their own algebraic foundations and do not require the quantum cognition literature for their justification. No optimization theory. The analysis is entirely representational: it characterizes the set of functions that each model class can represent, not whether gradient-based optimization can find the parameters that realize a given function. This is a significant practical limitation. The training of the proposed architecture involves several optimization challenges that are not addressed here: constrained optimization of the measurement matrix M on the Stiefel manifold absil2008optimization ; differentiation through complex-valued Cayley solves at each step; state-dependent Hamiltonian generation through a neural network gθg_θ that receives the current state as input, creating a recurrent dependency in the gradient computation; and complex-valued parameter gradients throughout. The combination of these factors may produce optimization landscapes with unfavorable local minima, poor conditioning, or gradient interference between parameter groups that negate the representational advantage identified by the separation theorem. The gradient norm preservation guaranteed by the unitary Cayley steps applies only to the state gradient, as noted above, not to the parameter gradients of gθg_θ. The absence of any optimization analysis, even heuristic, is a gap that significantly limits the practical relevance of the paper’s theoretical contributions. Addressing this gap, whether through theoretical analysis of the training dynamics or through empirical characterization of the loss landscape on the synthetic tasks of Section 6.2, is a necessary step before the architecture can be recommended for practical use. The comparison with real-valued models does not engage with the full literature on complex networks. Section 2 cites Trabelsi et al. trabelsi2018deep and Hirose and Yoshida hirose2012generalization as primary references for the advantages of complex-valued networks. The empirical literature is not uniformly positive on this point. Several studies have found that real-valued models with increased width, or with gating mechanisms such as LSTM hochreiter1997long gates or Mamba’s gu2023mamba selective state-space mechanism, match or outperform complex-valued models of the same parameter count on standard benchmarks. The representation-theoretic advantage identified by Theorem 5.10 applies to a specific comparison between Born-rule and affine-softmax readouts; it does not imply that complex-valued models outperform real-valued models with richer nonlinear readouts or gating mechanisms in all settings. A fair empirical evaluation of the proposed architecture would need to control for parameter count, training cost, and the choice of competing architectures, and would need to acknowledge the mixed existing evidence on complex-valued networks. 7 Conclusion This paper introduced a sequence modeling framework in which the latent state is a unit-norm vector in ℂNC^N evolving under a learned, time-dependent Hermitian Hamiltonian, and in which output probabilities are computed via the Born rule. The framework rests on a single structural constraint: the Hamiltonian H(t)=H0+Hint(t)H(t)=H_0+H_int(t) is Hermitian by construction (Section 3.2), which guarantees that the evolution is unitary, which guarantees that the state norm is exactly preserved, which guarantees that the Born rule produces a valid probability distribution. Each subsequent component of the architecture follows from this chain: the interaction picture (Section 3.3) factors out the analytically known free oscillations to isolate the input-driven dynamics; the Cayley (Crank–Nicolson) discretization (Section 3.4) transfers the continuous-time unitarity guarantee to discrete hardware without approximation hairer2006geometric ; crank1947practical ; and the Born-rule decoding (Section 3.5) converts the preserved complex state into token probabilities through a quadratic function that is sensitive to all (N2) N2 pairwise phase relationships among latent dimensions. It should be noted that norm preservation and integration accuracy are distinct properties: the Cayley transform guarantees that the state norm is exactly 11 at every step regardless of step size, but it does not eliminate phase errors that accumulate over long sequences. The model can process arbitrarily long sequences without norm drift, but the accuracy of the state trajectory as an approximation to the continuous-time solution degrades with sequence length at a rate determined by the step size and the spectrum of the interaction Hamiltonian HintH_int. The paper established three results from this construction. First, the Hermitian constraint produces a continuity equation (Section 4), mirroring the standard quantum-mechanical treatment of probability current sakurai1994modern , whose flux terms are pairwise probability currents Jj←k(t)=2Im(Hjkcj∗ck)J_j← k(t)=2Im(H_jk\,c_j^*c_k) between latent dimensions. These currents are antisymmetric, vanish on the diagonal, and sum to zero at each dimension (Proposition 4.1), providing an exact, conservation-respecting decomposition of each state update into directed pairwise transfers. The currents are driven entirely by the interaction Hamiltonian Hint(t)H_int(t) (Section 4.3), which links each token’s contribution to a measurable redistribution of amplitude across the model’s internal state. For diagnostic purposes, the midpoint current formula (37), evaluated at the average of the pre-step and post-step states, provides an exact pairwise decomposition of the discrete probability change without the O(Δt2)O( t^2) approximation error of the continuous-time formula. Second, the Born rule’s quadratic dependence on the complex amplitudes enables a representational advantage over the affine-softmax readout used by standard real-valued models. Theorem 5.10 constructs a family of disambiguation tasks ND_N for which a complex unitary sequence model of dimension N computes the correct output distribution exactly, while any real orthogonal sequence model with affine-softmax readout yang2018breaking requires dimension at least N2−2N^2-2. This separation is between two specific output mechanisms (Born-rule decoding and affine-softmax decoding) applied to norm-preserving dynamics of the respective types. It should not be read as a general claim that complex-valued dynamics are superior to real-valued dynamics for all architectures and readout mechanisms; a real-valued model with a quadratic or mixture-of-softmax readout could substantially narrow the gap. The proof traces the separation to a specific algebraic mechanism: the Born rule implicitly lifts the N-dimensional complex state to the N2N^2-dimensional space of Hermitian matrices via the map |ψ⟩↦|ψ⟩⟨ψ| ψ ψ ψ (Lemma 5.7; this is the same Veronese-type lifting exploited by Born machines in generative modeling cheng2018information ), accessing N2N^2 linearly independent features through the cross terms cjcj′∗c_jc_j ^* that encode relative phases. The affine-softmax readout, which depends linearly on the state, cannot synthesize these pairwise features from a state of dimension less than N2−2N^2-2 (Lemma 5.9). The separation theorem is established for the simplified state-independent CUSM; extending it to the full state-dependent architecture of Section 3 is an open problem. Additionally, Theorem 5.10 is stated conditionally on the assumption that the target log-probability matrix L∗L^* has full row rank N2N^2; this assumption is not derived from the task construction, and the Remark following the theorem identifies a log-odds reformulation as a path toward unconditionalizing the result. Lemma 5.6 provides a full explicit construction only for N=2N=2; the general-N argument invokes a spanning property of the Veronese-type variety that is stated but not fully executed. Third, the per-step computational cost of the architecture is O(Cgθ+Nr2+NV)O(C_g_θ+Nr^2+NV) (Section 6.1), where the Born-rule output cost O(NV)O(NV) matches the output-projection cost of standard architectures with the same vocabulary size, and the Cayley solve cost O(Nr2)O(Nr^2) is controlled by the rank parameter r of the interaction Hamiltonian. The numerical stability of the Woodbury solve at each step depends on the conditioning of an r×r× r matrix whose eigenvalues are determined by the current Hamiltonian; this conditioning is not guaranteed by the architecture and requires monitoring or regularization in practice golub2013matrix . These three results are interdependent in a way that reflects the coherence of the framework rather than being independent contributions. The probability currents exist because the dynamics are unitary, which is the same property that makes the Born rule produce valid probabilities. The separation theorem quantifies the advantage of Born-rule decoding, which depends on the phase information that the unitary dynamics preserve. The computational cost of the Cayley discretization is what makes the unitarity guarantee practical, and the low-rank structure of the interaction Hamiltonian that controls this cost is the same structure that decomposes the probability currents into interpretable coupling channels (equation (34)). The framework is theoretical in its present form: no experiments on language data have been conducted, and the separation theorem applies to a constructed task family rather than to natural language directly. The claim that language modeling may benefit from interference mechanisms is motivated by an analogy with the quantum cognition literature busemeyer2012quantum ; pothos2013quantum and by emerging empirical evidence of Born-rule effectiveness in NLP tasks guidotti2022text , which document interference effects in human judgment and quantum-inspired text classification respectively; these analogies are hypotheses to be empirically tested, not logical derivations from those results. Whether the interference mechanism that drives the separation is relevant to the structure of natural language is the central empirical question that this work leaves open. The experimental protocols specified in Section 6.2 are designed to resolve this question: they test whether the quadratic dimensional scaling predicted by Theorem 5.10 manifests in gradient-trained models (Prediction 1), whether Born-rule decoding outperforms softmax decoding on a shared complex state (Prediction 2), whether probability currents concentrate at semantically disambiguating tokens (Prediction 3), whether the learned frequency spectrum correlates with linguistic timescales (Prediction 4), and whether phase cross terms contribute measurably to output quality (Prediction 5). A positive outcome on Predictions 1 and 2 would confirm the representational advantage identified by the separation theorem in the context of trained models. A positive outcome on Predictions 3 through 5 would provide evidence that the interference mechanism is not merely a theoretical possibility but is exploited by the trained model on natural text. The paper’s contribution is a formal argument that the mathematical structure characterized by complex-valued unit-norm states, Hermitian generators, and squared-amplitude measurement provides specific and quantifiable algebraic properties that distinguish it from standard real-valued sequence models with linear readouts: norm-preservation by construction, a quadratic output mechanism that accesses O(N2)O(N^2) pairwise features from an N-dimensional state, and a conserved continuity equation that decomposes each state update into antisymmetric pairwise probability flows. Whether these algebraic properties translate to practical gains on language modeling benchmarks is a question that the theory motivates but cannot answer alone, and that the specified experimental program is designed to address. References [1] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances in Neural Information Processing Systems, volume 30, 2017. [2] Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural Computation, 9(8):1735–1780, 1997. [3] Albert Gu and Tri Dao. Mamba: Linear-time sequence modeling with selective state spaces. In The Twelfth International Conference on Learning Representations, 2024. [4] Jerome R. Busemeyer and Peter D. Bruza. Quantum Models of Cognition and Decision. Cambridge University Press, Cambridge, 2012. [5] Emmanuel M. Pothos and Jerome R. Busemeyer. Can quantum probability theory be a viable explanatory tool for human judgment and decision making? Psychological Review, 120(4):914–956, 2013. [6] Ricky T.Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. In Advances in Neural Information Processing Systems, volume 31, 2018. [7] Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. In Advances in Neural Information Processing Systems, volume 32, 2019. [8] Eldad Haber and Lars Ruthotto. Stable architectures for deep neural networks. Inverse Problems, 34(1):014004, 2017. [9] Albert Gu, Karan Goel, and Christopher Ré. Efficiently modeling long sequences with structured state spaces. In The Tenth International Conference on Learning Representations, 2022. [10] Ernst Hairer, Christian Lubich, and Gerhard Wanner. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, volume 31 of Springer Series in Computational Mathematics. Springer, Berlin, 2nd edition, 2006. [11] John Crank and Phyllis Nicolson. A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Mathematical Proceedings of the Cambridge Philosophical Society, 43(1):50–67, 1947. [12] Martin Arjovsky, Amar Shah, and Yoshua Bengio. Unitary evolution recurrent neural networks. In Proceedings of the 33rd International Conference on Machine Learning, pages 1120–1128, 2016. [13] Scott Wisdom, Thomas Powers, John R. Hershey, Jonathan Le Roux, and Les Atlas. Full-capacity unitary recurrent neural networks. In Advances in Neural Information Processing Systems, volume 29, 2016. [14] P.-A. Absil, Robert Mahony, and Rodolphe Sepulchre. Optimization Algorithms on Matrix Manifolds. Princeton University Press, Princeton, NJ, 2008. [15] Zakaria Mhammedi, Andrew Hellicar, Ashfaqur Rahman, and James Bailey. Efficient orthogonal parametrisation of recurrent neural networks using householder reflections. In Proceedings of the 34th International Conference on Machine Learning, pages 2401–2409, 2017. [16] Mario Lezcano-Casado and David Martínez-Rubio. Cheap orthogonal constraints in neural networks: A simple parametrization of the orthogonal and unitary group. In Proceedings of the 36th International Conference on Machine Learning, pages 3794–3803, 2019. [17] Eugene Vorontsov, Chiheb Trabelsi, Samuel Kadoury, and Chris Pal. On orthogonality and learning recurrent networks with long term dependencies. arXiv preprint arXiv:1702.00071, 2017. [18] Jimmy T.H. Smith, Andrew Warrington, and Scott W. Linderman. Simplified state space layers for sequence modeling. In The Eleventh International Conference on Learning Representations, 2023. [19] Michael Poli, Stefano Massaroli, Eric Nguyen, Daniel Y. Fu, Tri Dao, Stephen Baccus, Yoshua Bengio, Stefano Ermon, and Christopher Ré. Hyena hierarchy: Towards larger convolutional language models. In Proceedings of the 40th International Conference on Machine Learning, pages 28043–28078, 2023. [20] Chiheb Trabelsi, Olexa Bilaniuk, Ying Zhang, Dmitriy Serdyuk, Sandeep Subramanian, João Felipe Santos, Soroush Mehri, Negar Rostamzadeh, Yoshua Bengio, and Christopher J. Pal. Deep complex networks. In International Conference on Learning Representations, 2018. [21] Akira Hirose and Shotaro Yoshida. Generalization Characteristics of Complex-Valued Feedforward Neural Networks in Relation to Signal Coherence, volume 23. IEEE, 2012. [22] Patrick Virtue, X. Yu Stella, and Michael Lustig. Better than real: Complex-valued neural nets for MRI fingerprinting. In Proceedings of the IEEE International Conference on Image Processing, pages 3953–3957, 2017. [23] Emilien Dupont, Arnaud Doucet, and Yee Whye Teh. Augmented neural ODEs. In Advances in Neural Information Processing Systems, volume 32, 2019. [24] Patrick Kidger, James Morrill, James Foster, and Terry Lyons. Neural controlled differential equations for irregular time series. In Advances in Neural Information Processing Systems, volume 33, pages 6696–6707, 2020. [25] Edwin Stoudenmire and David J. Schwab. Supervised learning with quantum-inspired tensor networks. In Advances in Neural Information Processing Systems, volume 29, 2016. [26] Ramin Hasani, Mathias Lechner, Alexander Amini, Daniela Rus, and Radu Grosse. Liquid time-constant networks. In Proceedings of the 35th AAAI Conference on Artificial Intelligence, pages 7657–7666, 2021. [27] Mathias Lechner, Ramin Hasani, Alexander Amini, Thomas A. Henzinger, Daniela Rus, and Radu Grosse. Neural circuit policies enabling auditable autonomy. Nature Machine Intelligence, 2(10):642–652, 2020. [28] Ramin Hasani, Mathias Lechner, Alexander Amini, Philipp Cellier, Daniela Rus, and Radu Grosse. Closed-form continuous-time neural networks. Nature Machine Intelligence, 4(11):992–1003, 2022. [29] Emanuele Guidotti and Alfio Ferrara. Text classification with Born’s Rule. In Advances in Neural Information Processing Systems, volume 35, 2022. [30] Maria Schuld and Nathan Killoran. Quantum machine learning in feature Hilbert spaces. Physical Review Letters, 122(4):040504, 2019. [31] John Preskill. Quantum computing in the NISQ era and beyond. Quantum, 2:79, 2018. [32] Marco Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C. Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R. McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, and Patrick J. Coles. Variational quantum algorithms. Nature Reviews Physics, 3(9):625–644, 2021. [33] Kosuke Mitarai, Makoto Negoro, Masahiro Kitagawa, and Keisuke Fujii. Quantum circuit learning. Physical Review A, 98(3):032309, 2018. [34] Jacob Biamonte, Peter Wittek, Nicola Pancotti, Patrick Rebentrost, Nathan Wiebe, and Seth Lloyd. Quantum machine learning. Nature, 549(7671):195–202, 2017. [35] Vojtěch Havlíček, Antonio D. Córcoles, Kristan Temme, Aram W. Harrow, Abhinav Kandala, Jerry M. Chow, and Jay M. Gambetta. Supervised learning with quantum-enhanced feature spaces. Nature, 567(7747):209–212, 2019. [36] Amira Abbas, David Sutter, Christa Zoufal, Aurelien Lucchi, Alessio Figalli, and Stefan Woerner. The power of quantum neural networks. Nature Computational Science, 1(6):403–409, 2021. [37] Jarrod R. McClean, Sergio Boixo, Vadim N. Smelyanskiy, Ryan Babbush, and Hartmut Neven. Barren plateaus in quantum neural network training landscapes. Nature Communications, 9(1):4812, 2018. [38] Marco Cerezo, Akira Sone, Tyler Volkoff, Lukasz Cincio, and Patrick J. Coles. Cost function dependent barren plateaus in shallow parametrized quantum circuits. Nature Communications, 12(1):1791, 2021. [39] Adrián Pérez-Salinas, Alba Cervera-Lierta, Elies Gil-Fuster, and José I. Latorre. Data re-uploading for a universal quantum classifier. Quantum, 4:226, 2020. [40] Keisuke Fujii and Kohei Nakajima. Harnessing disordered-ensemble quantum dynamics for machine learning. Physical Review Applied, 8(2):024030, 2017. [41] Sanjib Ghosh, Andrzej Opala, Michał Matuszewski, Tomasz Paterek, and Timothy C. H. Liew. Quantum reservoir processing. npj Quantum Information, 5(1):35, 2019. [42] Herbert Jaeger and Harald Haas. Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication. Science, 304(5667):78–80, 2004. [43] Ewin Tang. A quantum-inspired classical algorithm for recommendation systems. In Proceedings of the 51st Annual ACM Symposium on Theory of Computing, pages 217–228, 2019. [44] Angelos Katharopoulos, Apoorv Vyas, Nikolaos Pappas, and François Fleuret. Transformers are RNNs: Fast autoregressive transformers with linear attention. In Proceedings of the 37th International Conference on Machine Learning, pages 5156–5165, 2020. [45] Manzil Zaheer, Guru Guruganesh, Kumar Avinava Dubey, Joshua Ainslie, Chris Alberti, Santiago Ontanon, Philip Pham, Anirudh Ravula, Qifan Wang, Li Yang, and Amr Ahmed. Big bird: Transformers for longer sequences. In Advances in Neural Information Processing Systems, volume 33, pages 17283–17297, 2020. [46] Tri Dao, Daniel Y. Fu, Stefano Ermon, Atri Rudra, and Christopher Ré. FlashAttention: Fast and memory-efficient exact attention with IO-awareness. In Advances in Neural Information Processing Systems, volume 35, pages 16344–16359, 2022. [47] Jun John Sakurai. Modern Quantum Mechanics. Addison-Wesley, Reading, MA, revised edition, 1994. [48] Arieh Iserles, Hans Z. Munthe-Kaas, Syvert P. Nørsett, and Antonella Zanna. Lie-group methods. Acta Numerica, 9:215–365, 2000. [49] Benedict Leimkuhler and Sebastian Reich. Simulating Hamiltonian Dynamics. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, Cambridge, 2004. [50] Zhilin Yang, Zihang Dai, Ruslan Salakhutdinov, and William W. Cohen. Breaking the softmax bottleneck: A high-rank RNN language model. In International Conference on Learning Representations, 2018. [51] James Townsend, Niklas Koep, and Sebastian Weichwald. Pymanopt: A Python toolbox for optimization on manifolds using automatic differentiation. In Journal of Machine Learning Research, volume 17, pages 1–5, 2016. [52] Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015. [53] David J. Griffiths. Introduction to Quantum Mechanics. Pearson Prentice Hall, Upper Saddle River, NJ, 2nd edition, 2005. [54] Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman. Deep inside convolutional networks: Visualising image classification models and saliency maps. arXiv preprint arXiv:1312.6034, 2014. [55] Yonatan Belinkov. Probing classifiers: Promises, shortcomings, and advances. Computational Linguistics, 48(1):207–219, 2022. [56] Joseph M. Renes, Robin Blume-Kohout, Andrew J. Scott, and Carlton M. Caves. Symmetric informationally complete quantum measurements. Journal of Mathematical Physics, 45(6):2171–2180, 2004. [57] Michael A. Nielsen and Isaac L. Chuang. Quantum Computation and Quantum Information. Cambridge University Press, Cambridge, 2000. [58] Song Cheng, Jing Chen, and Lei Wang. Information perspective to probabilistic modeling: Boltzmann machines versus born machines. Entropy, 20(8):583, 2018. [59] Jin-Guo Liu and Lei Wang. Differentiable learning of quantum circuit born machines. Physical Review A, 98(6):062324, 2018. [60] Gene H. Golub and Charles F. Van Loan. Matrix Computations. Johns Hopkins University Press, Baltimore, MD, 4th edition, 2013. [61] George A. Miller, Claudia Leacock, Randee Tengi, and Ross T. Bunker. A semantic concordance. In Proceedings of the ARPA Workshop on Human Language Technology, pages 303–308, 1993.