# Clustering in Deep Stochastic Transformers

Lev Fedorov, Michaël E. Sander, Romuald Elie, Pierre Marion, Mathieu Laurière

<figure><img src="image_1.png" /></figure>

## To cite this version:

Lev Fedorov, Michaël E. Sander, Romuald Elie, Pierre Marion, Mathieu Laurière. Clustering in Deep Stochastic Transformers. 2026. hal-05486973

HAL Id: hal-05486973

https://hal.science/hal-05486973v1

Preprint submitted on 31 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Clustering in Deep Stochastic Transformers

Lev Fedorov<fnref n="1" /> Michaël E. Sander² Romuald Elie² Pierre Marion³ Mathieu Laurière<fnref n="1" />

**Abstract**

Transformers have revolutionized deep learning across various domains but understanding the precise token dynamics remains a theoretical challenge. Existing theories of deep Transformers with layer normalization typically predict that tokens cluster to a single point; however, these results rely on deterministic weight assumptions, which fail to capture the standard initialization scheme in Transformers. In this work, we show that accounting for the intrinsic stochasticity of random initialization alters this picture. More precisely, we analyze deep Transformers where noise arises from the random initialization of value matrices. Under diffusion scaling and token-wise RMS normalization, we prove that, as the number of Transformer layers goes to infinity, the discrete token dynamics converge to an interacting-particle system on the sphere where tokens are driven by a *common* matrix-valued Brownian noise. In this limit, we show that initialization noise prevents the collapse to a single cluster predicted by deterministic models. For two tokens, we prove a phase transition governed by the interaction strength and the token dimension: unlike deterministic attention flows, antipodal configurations become attracting with positive probability. Numerical experiments confirm the predicted transition, reveal that antipodal formations persist for more than two tokens, and demonstrate that suppressing the intrinsic noise degrades accuracy.

**1. Introduction**

The Transformer architecture (Vaswani et al., 2017) has revolutionized modern machine learning, establishing itself as the backbone of state-of-the-art models across distinct modalities (Devlin et al., 2019; Dosovitskiy et al., 2021; Moussad et al., 2023). Despite this empirical dominance, a rigorous understanding of the signal propagation in deep,

randomly initialized Transformers remains limited. The attention mechanism (Bahdanau et al., 2014) introduces a global, nonlinear interaction whose asymptotic behavior in deep networks is understood only in fragments.

Among existing theoretical directions, mean-field and infinite depth limits of self-attention dynamics (Sander et al., 2022) have yielded a rich framework, explaining the emergence of clustering (Geshkovski et al., 2023; Criscitiello et al., 2024), metastability (Bruno et al., 2024), and geometric structures that mirror empirical behavior (Castin et al., 2025; Burger et al., 2025). However, these works typically rely on idealized structural assumptions: they consider weights that are tied across the layers, or equivalently assume a fully deterministic layer dynamics, thereby neglecting the layer-wise randomness inherent to initialization.

In parallel, several theoretical studies have incorporated noise into attention models (Shalova & Schlichting, 2024; Balasubramanian et al., 2025). While these works capture qualitative phenomena such as phase transitions and metastability of token configurations, the injected noise is typically external and *idiosyncratic*, acting independently at the token level. Consequently, these models lack a self-contained derivation from the forward pass of implementable Transformer architectures and do not capture the *common-noise* structure induced by layer-wise random initialization.

A third relevant direction is diffusion limits for residual networks (Peluchetti & Favaro, 2020; 2021; Cont et al., 2023). These models capture the intrinsic stochasticity emerging from depth and random initialization, but cannot express the global token mixing of attention or incorporate practical normalization schemes like RMSNorm, which plays a crucial role in preventing explosion and fundamentally alters token dynamics (Karagodin et al., 2024; 2025).

We therefore lack a theoretical framework for Transformers that incorporates *intrinsic* stochasticity, that is noise arising directly from standard random initialization. Understanding this intrinsic noise is essential, as it corresponds to initialization schemes used in real-world applications. In this paper, we analyze a simplified self-attention-only Transformer in which randomness enters exclusively through the value matrices at initialization. Under diffusion scaling with RMS normalization, we prove that the discrete token dynamics converge to a system of stochastic differential

<footnote n="1">¹New York University, ²Google DeepMind, ³INRIA. Correspondence to: <lev.fedorov@nyu.edu>.</footnote>
equations (SDEs) constrained to the unit sphere. Building on this continuous-time limit, we show that intrinsic noise qualitatively alters token clustering: the interplay between attention temperature and token dimension leads to phase transitions that cannot occur in deterministic attention flows. These transitions offer a new perspective on rank collapse phenomena: While noise still drives tokens toward low-rank configurations, it simultaneously diversifies the possible terminal states, enabling pairwise separation. More precisely, we make the following contributions:

* **Architecturally grounded stochastic limit (Thm. 1).** We prove convergence in distribution of deep randomly initialized Transformers—where noise emerges from random initialization of Value matrices—to a system of SDEs on the sphere. Importantly, the limiting system is invariant to the initialization distribution, as long as it is centered with fixed variance.

* **Explicit phase transition boundaries for 2 tokens (Thm. 2).** We show that tokens can either converge to a single point or separate into two antipodal points on the sphere—and these are the only possible outcomes. The system can be in one of two regimes, and we derive explicit boundaries for this phase transition in terms of token dimension and the inverse softmax temperature.

* **Noise-induced accessibility of deterministic zero-probability states (Prop. 1).** We introduce noise directly into the deterministic dynamics studied in prior works, establishing a sharp phase transition at an explicit threshold: Below it, deterministic clustering to a single point persists; above it, noise enables tokens to diverge into antipodal configurations.

* **Numerical validation and robustness beyond the analytically tractable regime (Section 6).** We show that random initialization of value matrices leads to an improvement in accuracy on CIFAR-10 compared to constant weight initializations. We present numerical illustrations of predicted phase transitions and clustering behavior. We show that antipodal configurations persist beyond the two-tokens setting and probe consistency across varying time horizons. Rigorous convergence rates and theoretical guarantees beyond the two-token case remain open directions for future work.

**Notations.** Let $\mathbf{X} = (X^1, \dots, X^N)$ collect tokens $X^i \in \mathbb{R}^d$ (column vectors). $P_X^\perp[Y] = Y - \langle X, Y \rangle X$ projects onto the tangent space of $\mathbb{S}^{d-1}$ at $\mathbf{X}$. We denote convergence in distribution by $\xrightarrow{d}$, the space of continuous functions by $C(\mathbf{X}; Y)$, and a standard matrix-valued Brownian motion by $W_t \in \mathbb{R}^{d \times d}$. Finally, $\|\cdot\|$ is the Euclidean norm, $\mathcal{N}(0, 1)$ the standard normal, $\mathbb{N}_0$ the non-negative integers, and $\mathrm{Id}(d)$ the identity matrix.

## 2. Background and related work

We briefly review the theoretical frameworks that treat deep networks as continuous-time dynamical systems. These scaling limits provide the mathematical foundation for analyzing signal propagation in deep Transformers.

### 2.1. Residual Neural Networks

Residual networks (ResNets) (He et al., 2016) can be conveniently interpreted as discretizations of differential equations (Chen et al., 2018). Let $X_n \in \mathbb{R}^d$ denote the hidden representation at layer $n$. Consider the layer-wise update

$$X_{n+1} = X_n + \alpha_L V_n f(X_n), \quad n = 0, \dots, L-1, \quad (1)$$

where $\alpha_L$ is a depth-dependent scaling factor and $V_n \in \mathbb{R}^{d \times d}$ denotes the layer weight matrix, typically initialized randomly before training (e.g., with independent standard normal entries). As the depth $L \to \infty$, the choice of $\alpha_L$ determines the nature of the limiting dynamics, as summarized in Table 1.

<table><caption>Table 1. Scaling regimes and their continuous limits. Here, $W_t \in \mathbb{R}^{d \times d}$ is a matrix whose entries are independent scalar Brownian motions.</caption><thead><tr><th>α<sub>L</sub></th><th>INITIALIZATION</th><th>CONTINUOUS LIMIT</th></tr></thead><tbody><tr><td>1/L</td><td>V<sub>n</sub> = Id(d)</td><td>ODE: dX<sub>t</sub> = f(X<sub>t</sub>)dt</td></tr><tr><td>1/√L</td><td>V<sub>n</sub><sup>ij</sup> <sub>I,I,D</sub> N(0, 1)</td><td>SDE: dX<sub>t</sub><sup>⊥</sup> = f(X<sub>t</sub>)<sup>⊥</sup> dW<sub>t</sub></td></tr></tbody></table>

While the ODE limit ($\alpha_L = 1/L$) describes the flow of infinite-depth networks with near-identity weights (a.k.a. Neural ODE (Chen et al., 2018)), the SDE limit ($\alpha_L = 1/\sqrt{L}$) captures the intrinsic stochasticity arising from random initialization (Yang & Schoenholz, 2017; Cont et al., 2023; Marion et al., 2025). This “diffusion regime” is critical for understanding trainability and gradient stability in deep networks (Peluchetti & Favaro, 2020; Marion et al., 2025). Many works have studied signal propagation in deep residual networks, which include in particular Transformers, e.g. Zhang et al. (2019); Bachlechner et al. (2021); Hayou et al. (2021); Yang et al. (2024); Chizat (2025); Dong et al. (2025). See also references in He & Hofmann (2024). In this work, we study more specifically this stochastic perspective for Transformers, where the analysis is complicated by the global coupling between tokens.

### 2.2. Transformers and Self-Attention

The Transformer architecture (Vaswani et al., 2017) consists of stacked layers combining self-attention mechanisms, feedforward blocks, normalization, and residual connections. The defining component is the attention layer, which updates each token (i.e., a vector-valued representation of an
element in the input sequence) by aggregating information from the entire sequence.

At an abstract level, the attention mechanism associates to each token $X^i$ a weighted linear combination of the token representations,

$$A_{\beta}(X^i, \mathbf{X}; \Theta) = \sum_{j=1}^{N} w_{ij}(X^i, \mathbf{X}; \Theta) X^j,$$

where the weights $w_{ij}$ depend on pairwise interactions between $X^i$ and $X^j$ and are parameterized by $\Theta$. The precise functional form of these weights is architecture-dependent and will be specified below.

The output of the attention block (see Fig. 5) is then obtained by applying a linear transformation through the Value matrix $V \in \mathbb{R}^{d \times d}$, so that the full update takes the form

$$X^i \mapsto V A_{\beta}(X^i, \mathbf{X}; \Theta).$$

**Softmax Self-Attention (SA).** The standard attention mechanism introduced by Vaswani et al. (2017) normalizes exponential interactions via a softmax:

$$A_{\beta}^{S}(X^i, \mathbf{X}) := \frac{\sum_{j=1}^{N} \exp(\beta \langle QX^i, KX^j \rangle) X^j}{\sum_{j=1}^{N} \exp(\beta \langle QX^i, KX^j \rangle)}, \quad (2)$$

where $Q$ and $K$ are respectively the query and key matrices and $\beta$ is an inverse temperature parameter. This formulation constitutes the backbone of modern Transformer architectures and is the primary object of interest from a modeling standpoint.

**Unnormalized Self-Attention (USA).** Alongside SA, we consider its unnormalized counterpart,

$$A_{\beta}^{U}(X^i, \mathbf{X}) := \frac{1}{N} \sum_{j=1}^{N} \exp(\beta \langle QX^i, KX^j \rangle) X^j, \quad (3)$$

which has been widely used as a theoretical proxy in the analysis of attention dynamics (Sander et al., 2022; Geshkovski et al., 2023). USA preserves the same exponential interaction geometry as SA, while removing the normalization that couples tokens through the denominator.

Throughout the paper, statements formulated for $A_{\beta} \in \{A_{\beta}^{S}, A_{\beta}^{U}\}$ apply to both $A_{\beta}^{S}$ and $A_{\beta}^{U}$.

## 2.3. Related theoretical perspectives

A growing body of work has investigated deep Transformer dynamics through simplified yet mathematically tractable models. Three mainstream theoretical approaches have emerged as particularly relevant to the present work.

**Mean-Field Transformers.** The first approach studies Transformers in a mean-field regime, taking both the network depth and the number of tokens to infinity (Sander et al., 2022; Vuckovic et al., 2021). In this setting, attention dynamics can be expressed in terms of an empirical measure over tokens, leading to deterministic limiting equations that describe smoothness (Castin et al., 2023), universality (Furuva et al., 2025), token sample complexity (Bohbot et al., 2025), clustering (Geshkovski et al., 2023), and metastability (Geshkovski et al., 2024; Bruno et al., 2024; Castin et al., 2025). A key feature of this framework is that attention parameters are typically assumed to be fixed, with the Value matrix treated deterministically. As a result, stochasticity arising from layer-wise random initialization is not captured.

**Stochastic Transformers.** A second line of work introduces stochasticity directly at the level of token dynamics, often by modeling attention-induced interactions through Langevin-type equations (Shalova & Schlichting, 2024; Balasubramanian et al., 2025). In these models, each token is driven by an independent Brownian motion, leading to McKean-Vlasov limits and associated Fokker-Planck equations. While this approach captures important phenomena such as phase transitions, the injected noise is typically idiosyncratic and externally imposed, with no direct derivation from the discrete Transformer architecture.

**Signal propagation.** Another perspective on Transformers comes from the signal propagation literature, which aims at finding parameterizations for which the forward and backward passes are well-conditioned at initialization. Noci et al. (2022) show that the $\alpha_L = 1/\sqrt{L}$ scaling factor in front of the residual connection for attention reduces the magnitude of clustering in the attention layers (which they refer to as *rank collapse*). However, they still empirically observe partial clustering even with the $1/\sqrt{L}$ scaling (see Fig. 4 in their paper). Follow-up works propose Shaped Attention (Noci et al., 2023; He & Hofmann, 2024), which consists in an affine transform of the attention matrix, and obtains competitive results with respect to standard attention. Our work provides theoretical grounding on the rank collapse/clustering phenomenon at initialization with $1/\sqrt{L}$ scaling and i.i.d. Value matrices.

**Position of our work.** In contrast, our paper studies theoretically and numerically stochasticity that arises *intrinsically* from standard random initialization of deep Transformers. Randomness enters exclusively through the Value matrices, resulting in a diffusion limit driven by *common noise* for all tokens. This modeling choice is closed to practical deep Transformers implementations (Touvron et al., 2023; Mistral AI, 2025), where stochasticity arises from layer-wise random initialization rather than from externally injected noise. Therefore, the global coupling induced by attention
is preserved, leading to stochastic dynamics that are fundamentally distinct from models with token-wise noise.

## 3. Deep Stochastic Transformers

In this section, we describe the precise Transformer model that underlies our analysis. We work with a minimal yet expressive architecture built from three core components: (i) self-attention, (ii) residual connections with depth-dependent scaling, and (iii) token-wise RMS normalization. Position-wise feedforward blocks are deliberately omitted, as they act independently on each token and therefore do not contribute to inter-token interactions; this allows us to isolate the core mechanism of attention.

### 3.1. Attention Mechanism

As detailed in Section 2.2, we focus on Softmax Self-Attention (SA), defined in (2), and its unnormalized proxy (USA), defined in (3). We recall that throughout the rest of the paper, $A_β$ refers to either formulation unless specified. Importantly, we treat the Key ($K$) and Query ($Q$) matrices as fixed parameters absorbed into the attention function. While standard initialization involves randomness in all weight matrices, freezing $K$ and $Q$ allows the tractable derivation of the continuous-time limit in Section 4. As our results demonstrate, the intrinsic noise arising solely from $V$ is sufficient to break convergence of tokens to a single cluster. We leave the analysis of the fully stochastic setting for future research.

### 3.2. Residual Connections and Scaling

We adopt the diffusion scaling $α_L \sim 1/\sqrt{L}$ in the residual update (1), in direct analogy with the scaling regimes introduced in Table 1 for stochastic deep ResNets, leading to nontrivial stochastic limits as $L \to \infty$.

### 3.3. Token-Wise RMS Normalization

This operation constrains each token to the unit sphere, preventing norm explosion and imposing a geometric structure on the dynamics. RMS normalization (post-layer normalization) is applied to every token, independently: RMSNorm($X^i$) = $X^i/\|X^i\|$. This normalization projects the attention field onto the tangent space of the sphere, thereby modifying the induced flow (see Fig. 1) while preserving inter-token interactions. An analysis of normalization-induced geometric effects, with comparisons between different normalization schemes, is provided by Karagodin et al. (2025). RMSNorm is also a standard component of modern LLMs, including LLaMA (Touvron et al., 2023) and Mistral (Mistral AI, 2025).

<figure><img src="image_3.png" /><figcaption>Figure 1. Effect of RMS normalization on attention field. Left: unnormalized $A_β^U(X, X)$. Right: tangent projection $P_X^⊥[A_β^U(X, X)]$ constrains flow to sphere. Here $X = (X, X^1, X^2)$ and $β = 1$, $Q = K = \text{Id}(2)$.</figcaption></figure>

### 3.4. Deep Stochastic Transformer Model

The complete discrete-time model studied in this work integrates the three architectural elements introduced above: Attention mechanism $A_β$ (Section 3.1), Diffusion-scaled residual connection (Section 3.2), Token-wise RMS normalization (Section 3.3). Arranged as shown in Fig. 5,

$$X_{n+1}^i = \frac{X_n^i + \frac{1}{\sqrt{L}} V_{n+1} A_β(X_n^i, X_n)}{\|X_n^i + \frac{1}{\sqrt{L}} V_{n+1} A_β(X_n^i, X_n)\|}, n \in \mathbb{N}_0, \quad (4)$$

where $X_n = (X_n^1, \dots, X_n^N) \in (\mathbb{R}^d)^N$, $A_β \in \{A_β^S, A_β^U\}$ and $\{V_n\}_{n \ge 0} \subset \mathbb{R}^{d \times d}$ is a sequence of Value matrices, such that:

**Assumption 1.** The matrices $\{V_n\}_{n \ge 1} \subset \mathbb{R}^{d \times d}$ are i.i.d. with independent entries $\{v_n^{kl}\}_{k,l=1}^d$ satisfying

$$\mathbb{E}[v_n^{kl}] = 0, \quad \mathbb{E}[(v_n^{kl})^2] = \sigma^2, \quad |v_n^{kl}| \le M \quad \text{a.s.}$$

for some constant $M > 0$ independent of $n, k, l$.

Note that this assumption is satisfied by standard uniform and truncated normal initializations (Glorot & Bengio, 2010; He et al., 2015).

## 4. Continuous Time Limit

In this section, we describe the continuous-time limit of Deep Stochastic Transformers. The main technical part of our derivation comes from the nonlinearity introduced by RMS normalization.

### 4.1. Continuous-time limit of stochastic transformers

To pass to continuous time, we use the standard linear interpolation for every token trajectory $X^i$ and step $n \in \mathbb{N}$:

$$\begin{cases} \Delta X_n^i := X_{n+1}^i - X_n^i, & i \in \{1, \dots, N\}, \\ X_L^i(t) := X_n^i + \frac{t-n/L}{L} \Delta X_n^i, & \text{if } \frac{n}{L} \le t < \frac{n+1}{L}. \end{cases} \quad (5)$$
This yields the continuous-time process

$$X_L(t) = (X_L^1(t), \dots, X_L^N(t)) \in C(\mathbb{R}_+; (\mathbb{R}^d)^N).$$

We can now state the following diffusion limit theorem for deep normalized Transformers.

**Theorem 1** (Continuous-time limit of normalized attention). Suppose that Assumption 1 holds. Let $\mathbf{Y} = (Y^1, \dots, Y^N)$ be the solution to the Transformer SDEs, that is, the system of stochastic differential equations driven by a common matrix-valued Brownian motion $W_t \in \mathbb{R}^{d \times d}$:

$$dY_t^i = \frac{(1-d)\sigma^2}{2} \|A_\beta(Y_t^i, \mathbf{Y}_t)\|^2 Y_t^i dt + \\ + \sigma P_{Y_t^i}^\perp [dW_t A_\beta(Y_t^i, \mathbf{Y}_t)], \quad i = 1, \dots, N, \quad (6)$$

with initial condition $Y_0^i = X_0^i$. Then the linearly interpolated discrete dynamics $X_L$ defined in (5) converge in distribution<fnref n="1" /> to $\mathbf{Y}$ as $L \to \infty$, namely, $X_L \xrightarrow{d} \mathbf{Y}$. The latter convergence holds for $\mathbf{Y}$ and $X_L$ viewed as random elements of $C(\mathbb{R}_+; (\mathbb{R}^d)^N)$.

The proof of this result relies on the limiting properties of the nonlinear difference equation (4) and their characterization via a martingale problem (Kushner & Huang, 1981; Watanabe, 1984a;b). We defer the proof to App. B.3.

Thm. 1 is reminiscent of the celebrated Donsker invariance principle (Donsker, 1951), which states that random walks, when properly scaled, converge weakly to Brownian motions. In this context, we treat our discrete scheme (4) as a modified random walk on $(\mathbb{S}^{d-1})^N$. It is important to note that, unlike strong convergence, we do not require Gaussian steps. In fact, regardless of the specific random initialization, all will converge to the same limit as $L \to \infty$, provided the second moment is fixed, entries are centered and the support is bounded. This ensures convergence for any centered initialization of the Value matrices, with no need for Gaussian initialization.

It is worth saying a few words about the additional drift term in equation (6). A naive approach to stochasticizing the deep Transformer would be to introduce a matrix-valued Brownian motion directly into the continuous-time attention dynamics on the sphere:

$$dX_t^i = P_{X_t^i}^\perp [dW_t A_\beta(X_t^i, \mathbf{X}_t)],$$

in direct analogy with ResNets (see Table 1). However, such an equation does not preserve the spherical constraint. Even if the noise is tangential, applying Itô's formula to the evolution of $\|X_t^i\|^2$ given by the formula above reveals a nonzero radial drift generated by second-order terms, causing trajectories to leave the sphere. This is a classical phenomenon

in the theory of spherical Brownian motion (Stroock, 1971). As a simple illustration, consider $X_t = (X_t^1, X_t^2) \in \mathbb{S}^1$ evolving according to $dX_t = \frac{1}{\|X_t\|}(-X_t^2, X_t^1) dW_t$. Although $\langle X_t, dX_t \rangle = 0$, Itô's formula yields $d\|X_t\|^2 = dt$, so $X_t$ immediately leaves the sphere. By contrast, we have the following lemma.

**Lemma 1.** For any sphere-supported initial distribution $Y_0 \in (\mathbb{S}^{d-1})^N$, any solution of equation (6) remains on the unit sphere with probability 1.

## 4.2. Bridging Deterministic and Stochastic Dynamics

Beyond initialization, training naturally leads to dynamics that combine a deterministic component with a stochastic one. The intuition that training adds a deterministic component with strong correlations across successive layers is supported by numerical experiments in Cohen et al. (2021); Marion et al. (2025) for i.i.d. initialization, and by both numerical experiments and proofs in Marion et al. (2023) in the weight-tied setting. The hybrid model considered in this section, where the deterministic component corresponds to an identity value matrix shared across all layers, should be viewed as a first step toward understanding signal propagation after initialization.

To precisely trace how intrinsic noise alters the clustering behavior, we study a parametrized perturbation of the deterministic dynamics. We therefore introduce a *hybrid model*, in which the update step combines a deterministic drift of order $1/L$ with a stochastic perturbation of order $\varepsilon/\sqrt{L}$:

$$X_{n+1}^i = \frac{X_n^i + \omega_\varepsilon(n, L) A_\beta(X_n^i, \mathbf{X}_n)}{\|X_n^i + \omega_\varepsilon(n, L) A_\beta(X_n^i, \mathbf{X}_n)\|} \quad (7)$$

Here $n \in \mathbb{N}_0$, $X_0^i \in \mathbb{S}^{d-1}$, and

$$\omega_\varepsilon(n, L) = \frac{1}{L} + \varepsilon \frac{v_n}{\sqrt{L}},$$

where $\{v_n\} \subset \mathbb{R}$ are i.i.d. centered random variables with unit variance and bounded support, and $\varepsilon \ge 0$ controls the noise amplitude; $\varepsilon = 0$ corresponds to the deterministic case. Its continuous-time limit makes explicit how noise enters the drift and qualitatively alters the long-time clustering behavior.

The following two lemmas characterize this limit. The first shows that for $\varepsilon = 0$, we recover the deterministic attention flow studied in Geshkovski et al. (2025).

**Lemma 2.** When $\varepsilon = 0$, the rescaled dynamics (7) converge uniformly to the solution of

$$dY_t^i = P_{Y_t^i}^\perp [A_\beta(Y_t^i, \mathbf{Y}_t)] dt, \quad Y_0^i = X_0^i.$$

That is, $\sup_{i,n} \|X_n^i - Y_{n/L}^i\| \to 0$ as $L \to \infty$, where the supremum is taken over $1 \le i \le N$ and $0 \le n \le L - 1$.

<footnote n="1">¹Convergence in distribution is defined in App. B.2, see (12).</footnote>
The second lemma gives the stochastic limit for $\varepsilon > 0$, analogous to Thm. 1.

**Lemma 3.** Let $\mathbf{X}_L^\varepsilon(t)$ be the linear interpolation of (7) and let $\mathbf{Y}_t^\varepsilon$ solve the SDE

$$d\mathbf{Y}_t^i = \left[ \varepsilon^2 \mathcal{F}(\mathbf{Y}_t^i, \mathbf{A}_\beta(\mathbf{Y}_t^i, \mathbf{Y}_t)) + P_{\mathbf{Y}_t^i}^\perp[\mathbf{A}_\beta(\mathbf{Y}_t^i, \mathbf{Y}_t)] \right] dt \\ + \varepsilon P_{\mathbf{Y}_t^i}^\perp[\mathbf{A}_\beta(\mathbf{Y}_t^i, \mathbf{Y}_t)] dB_t, \quad (8)$$

where $\mathcal{F}(x, y) := \frac{3}{2}\langle x, y \rangle^2 x - \frac{1}{2}\|y\|^2 x - \langle x, y \rangle y$, and $B_t$ is a scalar Brownian motion. Then

$$\mathbf{X}_L^\varepsilon \xrightarrow{d} \mathbf{Y}^\varepsilon \quad \text{as } L \to \infty.$$

The latter one holds for $\mathbf{Y}^\varepsilon$ and $\mathbf{X}_L^\varepsilon$ viewed as random elements of $C(\mathbb{R}_+; (\mathbb{R}^d)^N)$.

This limiting SDE explicitly shows how the noise amplitude $\varepsilon$ reweights the drift (through $\varepsilon^2 \mathcal{F}$) and adds a diffusion term. In Section 5.3, we analyze the long-time behavior of (8) and demonstrate that a sufficiently large $\varepsilon$ provably switches the clustering outcome from a single point to the antipodal configuration.

# 5. Clustering

In this section, we begin by defining the notion of clustering in the stochastic setting. We then state the main results concerning the long-time behavior of Transformer SDEs and highlight the relationship between $(\beta, d)$ and the clustering behavior in the case $N = 2$.

## 5.1. Notion of Clustering

We first define clustering configurations for two tokens.

**Definition 1.** Let $Y^i$ and $Y^j$ ($i \neq j$) be two tokens among $N$ governed by the dynamics (6). We define **clustering to a single token** as the following random event:

$$A_1^{ij}: \quad \|Y_t^i - Y_t^j\| \to 0 \quad \text{as } t \to \infty. \quad (9)$$

Similarly, we define **clustering to two antipodal tokens** as the following random event:

$$A_2^{ij}: \quad \|Y_t^i + Y_t^j\| \to 0 \quad \text{as } t \to \infty. \quad (10)$$

When $N = 2$, we denote $A_1 = A_1^{12}$ and $A_2 = A_2^{12}$.

## 5.2. On the Relation Between $\beta$ and $d$ with Clustering

When the number of tokens is restricted to $N = 2$, the long-time behavior of the Transformer SDEs can be analyzed explicitly. To obtain closed-form phase-transition thresholds in terms of $(\beta, d)$, we introduce the following simplifying assumption, which is used only for this explicit calculation.

**Assumption 2.** For the purpose of explicit phase-boundary computations, the key and query matrices satisfy $Q^\top K = \text{Id}(d)$.

This identity specialization simplifies the interaction geometry and allows for closed-form expressions. Importantly, it is not required for the existence of the continuous-time limit or for the general formulation of the Transformer SDEs.

Under this assumption, explicit phase transitions can be characterized in the two-token case.

**Theorem 2.** Let $N = 2$ and let the tokens follow (6). Suppose that Assumption 2 holds. Then, $\mathbb{P}(A_1) + \mathbb{P}(A_2) = 1$. Moreover, $\mathbb{P}(A_1) > 0$ and furthermore:

* $d - 2 \geq \cosh(2\beta)$ implies $\mathbb{P}(A_2) = 0$,

* $d - 2 < \cosh(2\beta)$ implies $\mathbb{P}(A_2) > 0$.

The proof is deferred to App. C.3.

**Contrast with the deterministic case.** This phase transition is a direct consequence of intrinsic noise. In the deterministic setting with $K = Q = \text{Id}(d)$, tokens always converge to a single cluster for any $\beta, d$ (Criscitiello et al. (2024), see also Thm 6.1 in (Geshkovski et al., 2025)). Here, stochasticity unlocks a new regime: for $\beta > \frac{1}{2} \cosh^{-1}(d - 2)$, the antipodal configuration emerges as a possibility.

Importantly, Thm. 2 provides a partial solution to Problem 2.15 posed by Geshkovski et al. (2025), asking whether the deterministic clustering conclusions could be generalized to Transformers with random matrices $(Q, K, V)$.

## 5.3. Noise-Induced Clustering

The continuous-time limit of the hybrid model, given by Lem. 3, provides a transparent setting to quantify how noise reshapes clustering. For the two-token case, we can derive an explicit threshold that separates the single-cluster and antipodal regimes.

**Proposition 1 (Noise-induced clustering transition).** Let $A_\beta = A_\beta^U$, $N = 2$ and Assumption 2 hold. The solution of the limiting SDE (8) satisfies:

1. If $\varepsilon^2 < 2e^{-\beta}$, $\lim_{t \to \infty} \|Y_t^1 - Y_t^2\| = 0$, with probability 1, i.e., the tokens converge to a single cluster.

2. If $\varepsilon^2 > 2e^{-\beta}$, $\lim_{t \to \infty} \|Y_t^1 + Y_t^2\| = 0$, with probability 1, i.e., the tokens converge to the antipodal configuration.

Prop. 1 delivers a sharp, interpretable condition: a sufficiently large noise amplitude $\varepsilon$ (relative to the inverse temperature $\beta$) forces the system into the antipodal phase. This result complements Thm. 2 obtained for the full stochastic transformer. Numerical simulations (Fig. 4) confirm the theoretical threshold $\varepsilon_c(\beta) = \sqrt{2e^{-\beta}}$ and illustrate how the noise-controlled transition manifests in the discrete dynamics (7). Together, these findings underscore that the stochasticity inherent in random initialization is not a tech-
nical detail, but qualitatively alters the clustering landscape
of deep transformers.

As shown in App. C.1, the overlap diffusion—namely, the
one-dimensional diffusion governing the pairwise inner
product of two tokens—reveals why antipodal configura-
tions are reachable. In addition, apart from coincidence
and antipodality, the diffusion is strictly positive, making
intermediate states unstable. Near the boundaries, both drift
and diffusion vanish linearly.

6. Experimental Results

<figure><img src="image_4.png" /><figcaption>Figure 2. Effect of value-matrix initialization on CIFAR-10 validation accuracy. **Constant Weights** (magenta) correspond to *identity initialization* of the matrices *V* shared across all layers, while **Random Weights** (violet) correspond to *independent Gaussian initialization* of *V* in each layer. Curves show the mean validation accuracy over five independent training runs; shaded regions indicate the min. and max. accuracy across runs at each epoch.</figcaption></figure>

**Key Role of Randomness.** To legitimize our focus on ran-
dom initialization, we empirically demonstrate that such
initialization in deep self-attention stacks leads to better ac-
curacy compared to weight-tied initialization, using a simpli-
fied Vision Transformer trained on CIFAR-10 (Krizhevsky
et al., 2009). Our implementation is based on a public Py-
Torch reimplementation of Vision Transformer for CIFAR
datasets,<fnref n="2" /> following the original ViT architecture (Dosovit-
skiy et al., 2021) and using AutoAugment data augmenta-
tion (Cubuk & Zoph, 2018). The model consists of 128
self-attention layers with embedding dimension *d* = 96 and
RMS normalization. All models are trained for 300 epochs
with batch size 128 and linear learning-rate warmup over 8
epochs. Results are averaged over 5 random seeds. The key
and query matrices *K*, *Q* are initialized once, shared identi-
cally across all layers, and remain fully trainable throughout
optimization. We consider a single-head attention architec-
ture without any feed-forward (MLP) layers, so that each
block consists solely of attention followed by normalization,

exactly as in (4). The model is trained end-to-end using the
standard cross-entropy loss.

We compare two initialization for the value matrices $V$:
(i) **constant identity** initialization, where all layers share
the same value matrix initialized as the identity, and (ii)
**independent Gaussian** initialization, where each layer's
value matrix is initialized independently with i.i.d. Gaussian
entries of standard deviation $\sigma = 1$. Other architectural
choices, hyperparameters, and training settings are fixed.
Despite identical architectures and fully trainable parame-
ters, identity-initialized $V$ lead to a degradation in validation
accuracy. This indicates that independence at initialization
has a significant impact on the training dynamics of deep
attention stacks. Note that a similar observation is made by
He & Hofmann (2024, Fig. 23) in a language task.

**Numerical Verification of the Phase Transition.** We numerically investigate the phase transition predicted by Thm. 2 by directly simulating the *discrete Transformer dynamics* (4). Importantly, no time-discretization of the limiting SDE (6) is used; all experiments are performed with the original discrete-time model, rather than discretizing the continuous-time stochastic differential equation using methods like Euler-Maruyama. This distinction is crucial, as the properties we establish for the continuous-time dynamics, particularly with respect to behavior at infinity, are not guaranteed to hold in the discrete case.

We consider $N = 2$ tokens and dimensions $d \in \{4, \dots, 10\}$. For each pair $(\beta, d)$, we simulate 40,000 independent trajectories, initialized uniformly on the unit sphere. The discrete dynamics (4) are simulated with $L = 50,000$ and $T = 500$. At each layer, the value matrix $V_n$ is independently re-initialized with i.i.d. Gaussian entries of variance $\sigma^2 = 1$, while the key and query matrices are fixed and equal to the identity, $K = Q = \text{Id}(d)$. Clustering is detected at the final layer using inner-product thresholds: a trajectory is classified as *antipodal* if $\langle X_T^1, X_T^2 \rangle \le -1 + \varepsilon$, with $\varepsilon = 10^{-3}$. Empirical probabilities are obtained by averaging over all trajectories.

Fig. 3 shows the empirical probability of antipodal con-
vergence as a function of $\beta$ for each dimension $d$. A
sharp transition is observed near the theoretical boundary
$\beta_c(d) = \frac{1}{2} \cosh^{-1}(d-2)$ predicted by the continuous-time
analysis. This demonstrates that the phase transition derived
for the limiting SDE accurately describes the long-time be-
havior of the *discrete Transformer dynamics*. We checked
that the results are not sensitive to varying the time hori-
zon $T$ (see App. D.1 for the details).

**Multiple Tokens Case.** We extend the numerical study of the antipodal clustering regime to systems with $N \ge 3$ tokens. All experiments use a fixed depth $L = 100$ and time horizon $T = 500$, and are based on 10,000 indepen-

<footnote n="2">²https://github.com/omihub777/ViT-CIFAR</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Noise-induced phase transition for two tokens. Empirical probability of antipodal convergence as a function of $\beta$ for dimensions $d \in \{4, \dots, 10\}$. Each point is estimated from 40,000 independent trajectories of the discrete dynamics (4) with $L = 100$, the time horizon $T = 500$, and $\sigma = 1$. The dashed curve shows the theoretical phase boundary $\beta_c(d) = \frac{1}{2} \cosh^{-1}(d-2)$.</figcaption></figure>

dent trajectories per parameter setting. For $d=4$, we run simulations on a uniform grid $\beta \in [0, 8]$ with step size 0.5 and varying $N$, and record the empirical probability that the terminal configuration contains at least one antipodal pair. As shown in Fig. 6, this probability is getting strictly positive across different values of $N$, indicating that antipodal structure is not specific to the two-token case. Beyond mere occurrence, we emphasize persistence. For large $\beta$, we observe that clustering typically occurs well before $T=500$, after which the fraction of trajectories exhibiting antipodal structure remains stable over long times (Fig. 7), indicating long-lived antipodal states of the noisy dynamics. The stability of the antipodal configuration for more than 2 tokens remains an open question, though we conjecture that it is indeed stable, rather than merely meta-stable.

**Numerical verification of noise-induced threshold.** We illustrate the noise-induced clustering transition predicted by Prop. 1 using simulations of the discrete hybrid dynamics (7). Fig. 4 shows a phase diagram in the $(\beta, \varepsilon)$ plane obtained from $10^3$ independent trajectories with dimension $d=3$, depth $L=100$, and final time $T=50$. For each pair $(\beta, \varepsilon)$, we compute the empirical fractions of trajectories converging to a single cluster and to an antipodal configuration. The heatmap displays the maximum of these two fractions, with color indicating the dominant outcome (magenta: single cluster, violet: antipodal). The dashed curve $\varepsilon = \sqrt{2e^{-\beta}}$ is the theoretical threshold derived in Prop. 1. A sharp transition is observed across this curve, confirming that the analytically predicted threshold accurately captures the clustering behavior of the discrete dynamics.

**Rate of convergence and metastable behavior.** While a detailed analysis of metastability and convergence rates

is beyond the scope of the present work, we briefly comment on the qualitative behavior observed numerically. For small values of the inverse temperature $\beta$, in regimes where convergence to a single-cluster configuration dominates, the observed convergence rate is comparable to that of the corresponding deterministic dynamics (i.e., the deterministic specialization of (7) with $\varepsilon=0$, $A_\beta = A_\beta^S$, and $Q^\top K = \text{Id}(4)$). In contrast, for larger values of $\beta$, the stochastic dynamics exhibits a markedly faster convergence. Overall, within the parameter ranges considered here, we do not observe evidence of metastable behavior in the stochastic setting, in the sense of Geshkovski et al. (2024). Convergence plots illustrating these regimes are provided in App. D.2.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Noise-induced clustering threshold. Dominant clustering outcome of the discrete hybrid dynamics in the $(\beta, \varepsilon)$ plane. Colors encode the maximum of the empirical fractions of single-cluster (magenta) and antipodal (violet) configurations. Dashed curve $\varepsilon = \sqrt{2e^{-\beta}}$: theoretical threshold from Prop. 1.</figcaption></figure>

## Conclusion

In this paper, we studied deep Transformer dynamics in a regime where stochasticity arises *intrinsically* from standard random initialization of value matrices. By deriving a diffusion limit under RMS normalization, we showed that this initialization-induced noise is not a negligible perturbation but a structural component of signal propagation. In particular, intrinsic noise qualitatively changes clustering behavior, unlocking antipodal configuration regimes that are provably inaccessible in deterministic attention flows, and leads to sharp phase transitions governed by the attention temperature and dimension.

A primary direction for future work is to extend the proof of antipodal clustering to the multi-token setting ($N \ge 3$), as suggested by our numerical experiments. In addition, our analysis relies on a simplified architecture, with randomness isolated in the value matrices, serving as a first step toward a complete stochastic theory. Future research directions include extending this framework to capture the effects of random key and query matrices, incorporating feedforward blocks, and analyzing the implications of these initialization dynamics during training.
## Impact statement

This paper analyzes theoretical properties of deep Transformers where stochasticity arises from random initialization of value matrices. We do not foresee any specific ethical or societal implications arising directly from this work.

## References

Bachlechner, T., Majumder, B., Mao, H., Cottrell, G., and Auley, J. ReZero is all you need: Fast convergence at large depth. In de Campos, C. and Maathuis, M. (eds.), *Proceedings of the Thirty-Seventh Conference on Uncertainty in Artificial Intelligence*, volume 161, pp. 1352–1361. PMLR, 2021.

Bahdanau, D., Cho, K., and Bengio, Y. Neural machine translation by jointly learning to align and translate. *arXiv preprint arXiv:1409.0473*, 2014.

Balasubramanian, K., Banerjee, S., and Rigollet, P. On the structure of stationary solutions to mckean-vlasov equations with applications to noisy transformers. *arXiv preprint arXiv:2510.20094*, 2025.

Billingsley, P. *Convergence of probability measures*. John Wiley & Sons, 2013.

Bogachev, V. I. *Weak convergence of measures*, volume 234. American Mathematical Society Providence, 2018.

Bohbot, L., Letrouit, C., Peyré, G., and Vialard, F.-X. Token sample complexity of attention, 2025. URL https://arxiv.org/abs/2512.10656.

Bruno, G., Pasqualotto, F., and Agazzi, A. Emergence of meta-stable clustering in mean-field transformer models. *arXiv preprint arXiv:2410.23228*, 2024.

Burger, M., Kabri, S., Korolev, Y., Roith, T., and Weigand, L. Analysis of mean-field models arising from self-attention dynamics in transformer architectures with layer normalization. *Philosophical Transactions A*, 383(2298): 20240233, 2025.

Castin, V., Ablin, P., and Peyré, G. How smooth is attention? *arXiv preprint arXiv:2312.14820*, 2023.

Castin, V., Ablin, P., Carrillo, J. A., and Peyré, G. A unified perspective on the dynamics of deep transformers. *arXiv preprint arXiv:2501.18322*, 2025.

Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K. Neural ordinary differential equations. *Advances in neural information processing systems*, 31, 2018.

Chizat, L. The hidden width of deep resnets: Tight error bounds and phase diagrams. *arXiv preprint arXiv:2509.10167*, 2025.

Cohen, A.-S., Cont, R., Rossier, A., and Xu, R. Scaling properties of deep residual networks. In Meila, M. and Zhang, T. (eds.), *Proceedings of the 38th International Conference on Machine Learning*, volume 139 of *Proceedings of Machine Learning Research*, pp. 2039–2048. PMLR, 2021. URL https://proceedings.mlr.org/press/v139/cohen21b.html.

Cont, R., Rossier, A., and Xu, R. Asymptotic analysis of deep residual networks, 2023. URL https://arxiv.org/abs/2212.08199.

Criscitiello, C., Rebjock, Q., McRae, A. D., and Boumal, N. Synchronization on circles and spheres with nonlinear interactions. *arXiv preprint arXiv:2405.18273*, 2024.

Cubuk, E. D. and Zoph, B. Improving deep learning performance with autoaugment. Google Research Blog, June 2018. Accessed: 2026-01-05.

Devlin, J., Chang, M.-W., Lee, K., and Toutanova, K. Bert: Pre-training of deep bidirectional transformers for language understanding, 2019. URL https://arxiv.org/abs/1810.04805.

Dong, Y., Noci, L., Khodak, M., and Li, M. Attention retrieves, mlp memorizes: Disentangling trainable components in the transformer. *arXiv preprint arXiv:2506.01115*, 2025.

Donsker, M. D. An invariance principle for certain probability limit theorems. *Mem. Amer. Math. Soc.*, 6:12, 1951. ISSN 0065-9266,1947-6221.

Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., Uszkoreit, J., and Houlsby, N. An image is worth 16x16 words: Transformers for image recognition at scale, 2021. URL https://arxiv.org/abs/2010.11929.

Ferguson, T. S. *A course in large sample theory*. Routledge, 2017.

Furuya, T., de Hoop, M. V., and Peyré, G. Transformers are universal in-context learners. In *The Thirteenth International Conference on Learning Representations*, 2025.

Geshkovski, B., Letrouit, C., Polyanskiy, Y., and Rigollet, P. The emergence of clusters in self-attention dynamics. *Advances in Neural Information Processing Systems*, 36: 57026–57037, 2023.

Geshkovski, B., Koubbi, H., Polyanskiy, Y., and Rigollet, P. Dynamic metastability in the self-attention model. *arXiv preprint arXiv:2410.06833*, 2024.
Geshkovski, B., Letrouit, C., Polyanskiy, Y., and Rigollet, P. A mathematical perspective on transformers. *Bulletin of the American Mathematical Society*, 62(3):427–479, 2025.

Glorot, X. and Bengio, Y. Understanding the difficulty of training deep feedforward neural networks. In *Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics*, volume 9 of *Proceedings of Machine Learning Research*, pp. 249–256. PMLR, 2010.

Hayou, S., Clerico, E., He, B., Deligiannidis, G., Doucet, A., and Rousseau, J. Stable ResNet. In Banerjee, A. and Fukumizu, K. (eds.), *Proceedings of The 24th International Conference on Artificial Intelligence and Statistics*, volume 130, pp. 1324–1332. PMLR, 2021.

He, B. and Hofmann, T. Simplifying transformer blocks. In Kim, B., Yue, Y., Chaudhuri, S., Fragkiadaki, K., Khan, M., and Sun, Y. (eds.), *International Conference on Representation Learning*, volume 2024, pp. 8882–8910, 2024.

He, K., Zhang, X., Ren, S., and Sun, J. Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. In *Proceedings of the IEEE International Conference on Computer Vision (ICCV)*, pp. 1026–1034, 2015. doi: 10.1109/ICCV.2015.123.

He, K., Zhang, X., Ren, S., and Sun, J. Deep residual learning for image recognition. In *Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, June 2016.

Jacobsen, M. One-dimensional homogeneous diffusion. Unpublished lecture notes., 2008. URL https://web.math.ku.dk/~susanne/SummerSchool2008/coursematerial.html. Lecture notes from BioMedMath 2008: Stochastic Differential Equation Models with Applications to the Insulin-Glucose System and Neuronal Modelling, Middelfart, Denmark.

Karagodin, N., Polyanskiy, Y., and Rigollet, P. Clustering in causal attention masking. *Advances in Neural Information Processing Systems*, 37:115652–115681, 2024.

Karagodin, N., Ge, S., Polyanskiy, Y., and Rigollet, P. Normalization in attention dynamics. arXiv preprint arXiv:2510.22026, 2025.

Karlin, S. and Taylor, H. E. *A second course in stochastic processes*. Elsevier, 1981.

Krizhevsky, A., Hinton, G., et al. Learning multiple layers of features from tiny images. 2009.

Kurtz, T. G. Equivalence of stochastic equations and martingale problems. In *Stochastic analysis 2010*, pp. 113–130. Springer, 2010.

Kushner, H. and Huang, H. On the weak convergence of a sequence of general stochastic difference equations to a diffusion. *SIAM Journal on Applied Mathematics*, 40(3): 528–541, 1981.

Marion, P., Wu, Y.-H., Sander, M. E., and Biau, G. Implicit regularization of deep residual networks towards neural codes. *arXiv preprint arXiv:2309.01213*, 2023.

Marion, P., Fermanian, A., Biau, G., and Vert, J.-P. Scaling resnets in the large-depth regime. *Journal of Machine Learning Research*, 26(56):1–48, 2025. URL http://jmlr.org/papers/v26/22-0664.html.

Mistral AI. Mistral large. https://mistral.ai/technology/#models, 2025. Accessed: 2025-09-14.

Moussad, B., Roche, R., and Bhattacharya, D. The transformative power of transformers in protein structure prediction. *Proceedings of the National Academy of Sciences*, 120(32):e2303499120, 2023.

Noci, L., Anagnostidis, S., Biggio, L., Orvieto, A., Singh, S. P., and Lucchi, A. Signal propagation in transformers: Theoretical perspectives and the role of rank collapse. *Advances in Neural Information Processing Systems*, 35: 27198–27211, 2022.

Noci, L., Li, C., Li, M., He, B., Hofmann, T., Maddison, C. J., and Roy, D. The shaped transformer: Attention models in the infinite depth-and-width limit. *Advances in Neural Information Processing Systems*, 36:54250–54281, 2023.

Peluchetti, S. and Favaro, S. Infinitely deep neural networks as diffusion processes. In Chiappa, S. and Calandra, R. (eds.), *Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics*, volume 108 of *Proceedings of Machine Learning Research*, pp. 1126–1136. PMLR, 26–28 Aug 2020. URL https://proceedings.mlr.press/v108/peluchetti20a.html.

Peluchetti, S. and Favaro, S. Doubly infinite residual neural networks: a diffusion process approach. *Journal of Machine Learning Research*, 22(175):1–48, 2021.

Sander, M. E., Ablin, P., Blondel, M., and Peyré, G. Sinkformers: Transformers with doubly stochastic attention. In *International Conference on Artificial Intelligence and Statistics*, pp. 3515–3530. PMLR, 2022.
Shalova, A. and Schlichting, A. Solutions of stationary mckean-vlasov equation on a high-dimensional sphere and other riemannian manifolds. *arXiv preprint arXiv:2412.14813*, 2024.

Stroock, D. W. On the growth of stochastic integrals. *Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete*, 18(4):340–344, 1971.

Stroock, D. W. and Varadhan, S. S. *Multidimensional diffusion processes*. Springer, 2007.

Touvron, H., Lavril, T., Izacard, G., Martinet, X., Lachaux, M.-A., Lacroix, T., Rozière, B., Goyal, N., Hambro, E., Azhar, F., et al. Llama: Open and efficient foundation language models. *arXiv preprint arXiv:2302.13971*, 2023.

Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., and Polosukhin, I. Attention is all you need. *Advances in neural information processing systems*, 30, 2017.

Vuckovic, J., Baratin, A., and Combes, R. T. d. On the regularity of attention. *arXiv preprint arXiv:2102.05628*, 2021.

Watanabe, H. Diffusion approximations of some stochastic difference equations II. *Hiroshima Mathematical Journal*, 14:15–34, 1984a.

Watanabe, H. Diffusion approximations of some stochastic difference equations. ii. *Hiroshima mathematical journal*, 14(1):15–34, 1984b.

Yang, G. and Schoenholz, S. Mean field residual networks: On the edge of chaos. *Advances in neural information processing systems*, 30, 2017.

Yang, G., Yu, D., Zhu, C., and Hayou, S. Tensor programs VI: Feature learning in infinite depth neural networks. In *The Twelfth International Conference on Learning Representations*, 2024.

Zhang, H., Dauphin, Y., and Ma, T. Fixup initialization: Residual learning without normalization. In *International Conference on Learning Representations*, 2019.
A. Additional Figures

<figure><img src="image_7.png" /><figcaption>Figure 5. One layer of the Deep Transformer with added normalization and diffusion scaling (~1/√L). The diagram highlights that self-attention acts globally across all tokens. The residual connection is shown with a dashed outline to indicate that it represents a mode of connectivity rather than an operation.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Empirical probability of observing **at least one** antipodal pair in the terminal configuration as a function of $\beta$, for multiple values of $N$ at dimension $d = 4$. Simulations are performed on a uniform grid $\beta \in [0, 8]$ with step size 0.5, using 10,000 independent trajectories per parameter setting and a fixed time horizon $T = 500$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7. Time evolution of clustering fractions for a multi-token system ($d = 4, N = 50$) at inverse temperature $\beta = 5$. All curves are obtained by averaging over 2,000 independent trajectories and are shown up to time horizon $T = 50$.</figcaption></figure>
## B. Continuous Time Limit

First, we fix the basic notion of weak convergence used throughout this appendix (Appendix B.2). We then establish the continuous-time limit of the discrete Transformer dynamics under diffusion scaling, proving weak convergence to the limiting stochastic differential equation stated in Theorem 1 (Appendix B.3). Finally, we show that the limiting dynamics preserve the spherical constraint almost surely, so that trajectories remain on the unit sphere for all times (Lemma 1). At the end, we provide the proofs of Lemmas 2 and 3.

### B.1. Notation Convention

The space of infinitely smooth functions from $X$ to $Y$ is denoted by $C^{\infty}(X, Y)$. Here we deal with $(\mathbb{R}^d)^N$-valued diffusion processes and to simplify the notation, we treat $\mathbf{X} \in (\mathbb{R}^d)^N$ as column vector $\mathbf{X} = \begin{bmatrix} X^1 \\ \vdots \\ X^N \end{bmatrix}$, where $X^i \in \mathbb{R}^d$ is a column vector as well. Vector calculus notation for the gradient $\nabla\phi(\mathbf{X})$ and Hessian $\nabla^2\phi(\mathbf{X})$ changes accordingly.

### B.2. Background on Martingale Problem and Weak Convergence

We do not aim at a complete or self-contained exposition of the probabilistic framework used in this work. Instead, we follow the general setting developed by Stroock and Varadhan and only recall the basic definitions and notation that will be used throughout the paper. For a detailed and systematic treatment, we refer the reader to the classical references on weak convergence and stochastic processes (Stroock & Varadhan, 2007; Billingsley, 2013; Bogachev, 2018).

Consider the canonical probability space $\Omega = C(\mathbb{R}_+; \mathbb{R}^D)$ of continuous functions taking values in $\mathbb{R}^D$, equipped with the metric

$$d(\omega, \omega') := \sum_{n=0}^{+\infty} \frac{1}{2^n} \frac{\sup_{0 \le t \le n} \|x(t, \omega) - x(t, \omega')\|}{1 + \sup_{0 \le t \le n} \|x(t, \omega) - x(t, \omega')\|}, \quad (11)$$

where $x(t, \omega) := \omega(t)$ denotes the position of the trajectory $\omega$ at time $t$. With this metric, $(\Omega, d)$ is a Polish$^3$ space. Endowed with its Borel $\sigma$-algebra $\mathcal{B}(\Omega)$, this is the underlying measurable space on which we work.

Let $\{\mu_n\}_{n=1}^{\infty} \subset \mathcal{P}(\Omega)$ be a sequence of probability measures. We say that $\mu_n$ converges weakly to $\mu \in \mathcal{P}(\Omega)$, and write

$$\mu_n \xrightarrow{w} \mu,$$

if for every bounded continuous function $f: \Omega \to \mathbb{R}$,

$$\int f(\omega) \mu_n(d\omega) \to \int f(\omega) \mu(d\omega), \quad \text{as } n \to \infty.$$

If $X_n$ and $X$ are $\Omega$-valued random variables such that $\mu_n = \text{Law}(X_n)$, $\mu = \text{Law}(X)$, then the above convergence is called convergence in distribution. In this case, we write

$$X_n \xrightarrow{d} X. \quad (12)$$

**Martingale problem.** We recall the notion of the martingale problem in the sense of Stroock and Varadhan, adapted to the present setting. We emphasize that, although one may formulate the martingale problem together with an explicit filtration, in this work we do not specify a filtration a priori. Instead, throughout the paper we implicitly work with the natural filtration generated by the coordinate process.

Let $E = \mathbb{R}^D$ and let $\Omega$ be the canonical path space. Let $L$ be a (second-order) diffusion operator on $E$ with domain $\text{dom}(L) \subset C_b(E)$, and let $\nu \in \mathcal{P}(E)$ be a prescribed initial distribution.

**Definition 2.** A probability measure $P \in \mathcal{P}(\Omega)$ is said to be a solution of the martingale problem for $(L, \nu)$ if

1. the initial marginal of $P$ is $\nu$, that is

$$P \circ X_0^{-1} = \nu;$$

<footnote>$^3$That is, complete and separable; see Bogachev (2018).</footnote>
2. for every $f \in \text{dom}(L)$, the process

$$M_f(t) := f(X_t) - f(X_0) - \int_0^t Lf(X_s) ds, \quad t \ge 0,$$

is a P-martingale with respect to the natural filtration.

In this case, we write $P \in \text{MP}(L, \nu)$.

### B.3. Continuous-Time Limit of Deep Stochastic Transformers

*Proof of Theorem 1.* First, we carry out some preliminary steps. Note that, to establish weak convergence on $\Omega$, it suffices to prove it on an arbitrary but fixed finite time horizon $T \in \mathbb{R}_+$. This follows directly from the form of the metric (11). We therefore fix a time horizon $T$. We decompose the proof into the following steps:

1. Linearize the dynamics (4) and show that the linearized dynamics converges to the nonlinear one;

2. Show that the limit of the linearized dynamics is given by (6);

3. Transfer the convergence result from the linearized dynamics to the original nonlinear dynamics.

Now we implement the steps outlined above. For this purpose, we introduce an auxiliary function $\tilde{A}_\beta: \mathbb{R}^d \times (\mathbb{R}^d)^N \to \mathbb{R}^d$ such that $\tilde{A}_\beta \in C^\infty(\mathbb{R}^d \times (\mathbb{R}^d)^N, \mathbb{R}^d)$ and

$$\tilde{A}_{\beta}(X, \mathbf{X}) := \begin{cases} A_{\beta}, & \text{if } \|X\|, \|X^i\| \le 2 \text{ for } 1 \le i \le N, \\ 0, & \text{if } \|X\|, \|X^i\| \ge 4. \end{cases}$$

Clearly, this function is not unique, coincides with $A_\beta$ on $\mathbb{S}^{d-1} \times (\mathbb{S}^{d-1})^N$, and is globally bounded. We note that any result proved for the recurrent relation below

$$X_{n+1}^i = \frac{X_n^i + \frac{1}{\sqrt{L}} V_{n+1} \tilde{A}_\beta(X_n^i, \mathbf{X}_n)}{\|X_n^i + \frac{1}{\sqrt{L}} V_{n+1} \tilde{A}_\beta(X_n^i, \mathbf{X}_n)\|}, n \in \mathbb{N}_0, \mathbf{X}_0 \in (\mathbb{S}^{d-1})^N, \quad (13)$$

also holds for the recurrent relation (4), since the corresponding solutions coincide. At this stage, the introduction of the modified function merely allows us to state results on the entire space without additional comments, while they automatically apply to the compact subset of interest.

**Step 1. [Linearization]** We note that the following expansion is valid for $X \in \mathbb{S}^{d-1}$ and $\mathbf{X} \in (\mathbb{S}^{d-1})^N$:

$$\frac{X + \delta V \tilde{A}_{\beta}(X, \mathbf{X})}{\|X + \delta V \tilde{A}_{\beta}(X, \mathbf{X})\|} = X + \delta P_{\mathbf{X}}^{\perp} [V \tilde{A}_{\beta}(X, \mathbf{X})] + \delta^2 G(X, V \tilde{A}_{\beta}(X, \mathbf{X})) + \delta^3 R(X, V \tilde{A}_{\beta}(X, \mathbf{X})),$$

where $G$ is given by

$$G(x, y) := \frac{3}{2} \langle x, y \rangle^2 x - \frac{1}{2} \langle y, y \rangle x - \langle x, y \rangle y,$$

and $R(X, V \tilde{A}_{\beta}(X, \mathbf{X}))$ is a uniformly bounded remainder term. Discarding the remainder, we obtain a truncated discrete-time dynamics $\{Y_n^i\}_{n \in \mathbb{N}_0} \subset \mathbb{R}^d$ with $\delta = \frac{1}{\sqrt{L}}$ in place of (13):

$$Y_{n+1}^i = Y_n^i + \frac{1}{\sqrt{L}} P_{Y_n^i}^{\perp} [V_{n+1} \tilde{A}_\beta(Y_n^i, \mathbf{Y}_n)] + \frac{1}{L} G(Y_n^i, V_{n+1} \tilde{A}_\beta(Y_n^i, \mathbf{Y}_n)), \quad i \in \{1, \dots, N\}, \quad (14)$$

and we assume that $\mathbf{Y}_0 = \mathbf{X}_0$. Note that, for matrices $\{V_n\}$ with uniformly bounded entries, the functional families $\{G(X, V_n \tilde{A}_\beta(X, \mathbf{X}))\}_{n=0}^\infty$ and $\{P_{\mathbf{X}}^\perp [V_n \tilde{A}_\beta(X, \mathbf{X})]\}_{n=0}^\infty$ are uniformly Lipschitz$^4$, with Lipschitz constants $M_G$ and $M_P$, respectively. Then, a standard Grönwall-type argument implies that, for any fixed time horizon $T$,

$$\max_{1 \le i \le N} \max_{1 \le k \le [TL]} \mathbb{E} \|X_k^i - Y_k^i\|^2 \to 0 \quad \text{as } L \to \infty. \quad (15)$$

<footnote>$^4$That is, all functions in the family are M-Lipschitz for a common constant $M \in \mathbb{R}_{\ge 0}$.</footnote>
We carry out this argument below. Indeed, consider the difference $\Delta_n^i := X_n^i - Y_n^i$:

$$ \Delta_{n+1}^i = \Delta_n^i + \frac{1}{\sqrt{L}} \left( P_{X_n^i}^\perp [V_{n+1} \tilde{\Lambda}_\beta(X_n^i, \mathbf{X}_n)] - P_{Y_n^i}^\perp [V_{n+1} \tilde{\Lambda}_\beta(Y_n^i, \mathbf{Y}_n)] \right) + \\ + \frac{1}{L} \left( G(X_n^i, V_{n+1} \tilde{\Lambda}_\beta(X_n^i, \mathbf{X}_n)) - G(Y_n^i, V_{n+1} \tilde{\Lambda}_\beta(Y_n^i, \mathbf{Y}_n)) \right) + \frac{1}{L^{3/2}} R_n^i, $$

where $R_n^i = R(X_n^i, V^{n+1} \tilde{\Lambda}_\beta(X_n^i, \mathbf{X}_n))$. Moreover, we have

$$ \Delta_{n+1}^i = \frac{1}{\sqrt{L}} \sum_{k=1}^{n} \left( P_{X_k^i}^\perp [V_{k+1} \tilde{\Lambda}_\beta(X_k^i, \mathbf{X}_k)] - P_{Y_k^i}^\perp [V_{k+1} \tilde{\Lambda}_\beta(Y_k^i, \mathbf{Y}_k)] \right) + \\ + \frac{1}{L} \sum_{k=1}^{n} \left( G(X_k^i, V_{k+1} \tilde{\Lambda}_\beta(X_k^i, \mathbf{X}_k)) - G(Y_k^i, V_{k+1} \tilde{\Lambda}_\beta(Y_k^i, \mathbf{Y}_k)) \right) + \frac{1}{L^{3/2}} \sum_{k=1}^{n} R_k^i. $$

We also define the difference between all tokens as $\Delta_n := \mathbf{X}_n - \mathbf{Y}_n$. Note that, for the Euclidean norms in $\mathbb{R}^{d \times N}$ and $\mathbb{R}^d$, we have the following equality

$$ \|\Delta_n\|^2 = \sum_{k=1}^{n} \|\Delta_n^k\|^2. \qquad (16) $$

Now, applying the Cauchy-Schwarz inequality and using independence of the entries of $\{V_n\}$ for different values of $n$ (discrete Itô isometry), we obtain

$$ \mathbb{E}\|\Delta_{n+1}^i\|^2 \le 3 \left( \frac{M_P^2}{L} + \frac{M_G^2}{L} \right) \sum_{k=1}^{n} \mathbb{E}\|\Delta_k\|^2 + \frac{3R^2}{L^2}, $$

where $R = \max_{i,n} R_n^i$. Summing the inequality above and applying the property of the norm (16), we obtain

$$ \mathbb{E}\|\Delta_{n+1}\|^2 \le 3N \left( \frac{M_P^2}{L} + \frac{M_G^2}{L} \right) \sum_{k=1}^{n} \mathbb{E}\|\Delta_k\|^2 + \frac{3R^2 N}{L^2}. $$

Thus, by discrete Grönwall inequality we have $\mathbb{E}\|\Delta_{n+1}\|^2 \le \frac{3R^2 N}{L^2} \exp(3N(M_P^2 + M_G^2)T)$ for $n+1 \le \lfloor LT \rfloor$. Since $N$, $R$, $T$, $M_P$, and $M_G$ are treated as fixed constants, this yields the convergence (15).

Clearly, the same uniform approximation result holds for the linear interpolations of the sequences $\{X_n^i\}_{n \in \mathbb{N}_0} \subset \mathbb{R}^d$ and $\{Y_n^i\}_{n \in \mathbb{N}_0} \subset \mathbb{R}^d$, denoted respectively by $\mathbf{X}_L(t)$ and $\mathbf{Y}_L(t)$:

$$ \max_{1 \le i \le N} \sup_{0 \le t \le T} \mathbb{E}\|X_L^i(t) - Y_L^i(t)\|^2 \to 0 \quad \text{as } L \to \infty. \qquad (17) $$

We now proceed to the second step.

**Step 2. [Weak Convergence]** This step crucially relies on the work of Watanabe (1984a), which establishes convergence in terms of the martingale problem together with the uniqueness of the solution. Watanabe's theorem reduces the proof to identifying the limiting process generator via the first two conditional moments of the one-step increment.

Indeed, according to Theorem 1 therein, we can explicitly compute the coefficients of the generator $R$ corresponding to the solution of the limiting martingale problem $MP(R, \text{Law}(\mathbf{X}_0))$.

Assumptions I)-VI) of Watanabe (1984a) are satisfied under our setting for the following reasons. First, all coefficient functions involved in the scheme are infinitely differentiable; in particular, they are more regular than required, since the theorem only assumes a finite number of bounded derivatives. Second, the driving random matrices $\{V_n\}_{n \ge 1}$ form an i.i.d. sequence, which is stronger than the stationarity assumptions imposed in the general framework. Third, all derivatives of the coefficients are locally bounded, a direct consequence of the smoothness of $\tilde{\Lambda}_\beta$ together with the compactness of its support. Finally, the mean drift condition holds trivially, since the noise has zero mean due to the independence and centering of $\{V_n\}$, see Assumption 1.
Expectations for the the Assumption VII) could be computed explicitly using the following identities:

$$ \mathbb{E}\langle x, Vy \rangle^2 = \sigma^2 \cdot \|x\|^2 \cdot \|y\|^2, \quad \mathbb{E}\langle Vy, Vy \rangle = \sigma^2 \cdot d \cdot \|y\|^2, \quad \mathbb{E}\langle x, Vy \rangle Vy = \sigma^2 \cdot \|y\|^2 \cdot x \in \mathbb{R}^d. \quad (18) $$

Treating entire sequence of the tokens as a vector in $(\mathbb{R}^d)^N$, we get that

$$ c(\mathbf{X}) = \sigma^2 \frac{1-d}{2} \begin{bmatrix} \|\tilde{A}_\beta(X^1, \mathbf{X})\|^2 X^1 \\ \vdots \\ \|\tilde{A}_\beta(X^N, \mathbf{X})\|^2 X^N \end{bmatrix} \in (\mathbb{R}^d)^N, \quad b_{k\ell}(\mathbf{X}) = 0, $$

Computing of $A(\mathbf{X}) = (a_{k\ell}(\mathbf{X}))_{1\le k, \ell \le d \cdot N}$ is more delicate, because we have a diffusion in $(\mathbb{R}^d)^N$:

$$ A(\mathbf{X}) = \mathbb{E} \underbrace{\begin{bmatrix} P_{\tilde{X}^1}^{\perp} [V \tilde{A}_{\beta}(X^1, \mathbf{X})] \\ \vdots \\ P_{\tilde{X}^N}^{\perp} [V \tilde{A}_{\beta}(X^N, \mathbf{X})] \end{bmatrix}}_{\in \mathbb{R}^{d \times N}} \begin{bmatrix} [P_{\tilde{X}^1}^{\perp} [V \tilde{A}_{\beta}(X^1, \mathbf{X})]]^T & \cdots & [P_{\tilde{X}^N}^{\perp} [V \tilde{A}_{\beta}(X^N, \mathbf{X})]]^T \end{bmatrix} $$

A direct identification shows that the generator $R$, defined by

$$ R\phi(\mathbf{X}) = \langle c(\mathbf{X}), \nabla\phi(\mathbf{X}) \rangle + \frac{1}{2} \mathrm{Tr}[A(\mathbf{X})\nabla^2\phi(\mathbf{X})], $$

coincides with the generator $R'$ associated with the following system:

$$ dY_t^i = \frac{(1-d)\sigma^2}{2} \|\tilde{A}_\beta(Y_t^i, Y_t)\|^2 Y_t^i dt + \sigma P_{Y_t^i}^{\perp} [dW_t \tilde{A}_\beta(Y_t^i, Y_t)], \quad Y_0 = X_0, \quad i = 1, \dots, N, \quad (19) $$

treated as an $(\mathbb{R}^d)^N$ valued diffusion process. Moreover, we enjoy the fact that $MP(R', \mathrm{Law}(Y_0))$ is a weak solution of (19), see for instance Theorem 1 in Kurtz (2010) (local boundedness and measurability requirements are satisfied). As a bonus, we get that this weak solution is distributionally unique, because such is solution of the martingale problem.

By a standard argument, which we reproduce in Lemma 1, any solution of (19) remains almost surely supported on the unit sphere. Consequently, along such trajectories one has $\tilde{A}_\beta = A_\beta$ in (19). This completes the identification of the limit and shows that the linearized dynamics converges to (6).

**Step 3. [Passage to the Limit]** We apply a Slutsky-type argument; see Theorem 6b) in Ferguson (2017), which extends without modification to Polish spaces. Indeed, note that

1. convergence (17) implies convergence in probability: $X_L - Y_L \xrightarrow{\mathbb{P}} 0$;

2. according to Step 2, we have $Y_L \xrightarrow{d} Y$.

It follows that $(X_L - Y_L) + Y_L \xrightarrow{d} Y$, which completes the proof of the theorem. $\square$

*Proof of Lemma 1.* The result can be easily deduced by applying Itô's formula to $\|Y_t\|^2 = \langle Y_t, Y_t \rangle$. We immediately obtain

$$ d\langle Y_t^i, Y_t^i \rangle = 2\langle Y_t^i, dY_t^i \rangle + \langle dY_t^i, dY_t^i \rangle. $$

Using the fact that $\langle P_x^{\perp}[y], x \rangle = 0$ for all $y \in \mathbb{R}^d$ and performing some simplifications, we obtain

$$ d\|Y_t^i\|^2 = \sigma^2 \|A_\beta(Y_t^i, Y_t)\|^2 (\|Y_t^i\|^2 - d)(\|Y_t^i\|^2 - 1) dt, $$

which implies that $Y_{t+s}^i \in \mathbb{S}^{d-1}$ for all $s > 0$ as long as $Y_t^i \in \mathbb{S}^{d-1}$. Since we have the initial condition $\|Y_0^i\| = 1$, we immediately obtain the desired statement. $\square$
Proof of Lemma 2. Similarly to the first step of the proof of Theorem 1, for $x \in \mathbb{S}^{d-1}$ and $y \in \mathbb{R}^d$ we consider the following function:

$$F_{\delta}(x, y) = (x + \delta y)/\|x + \delta y\| = \frac{(x + \delta y)}{\sqrt{\langle x, x \rangle + 2\delta\langle x, y \rangle + \delta^2\langle y, y \rangle}}.$$

Trivially, we obtain

$$F_{\delta}(x, y) = (x + \delta y) (1 + 2\delta A(x, y) + \delta^2 B(x, y))^{-1/2},$$

where we denote

$$A(x, y) := \langle x, y \rangle, \quad B(x, y) := \langle y, y \rangle.$$

By directly applying the Taylor expansion and assuming that $x \in \mathbb{S}^{d-1}$ and $y$ belongs to a compact subset $K \subset \mathbb{R}^d$, we obtain a more convenient form for sufficiently small $\delta > 0$:

$$F_{\delta}(x, y) = x + \delta P_{x}^{\perp}[y] + \delta^2 \left( \frac{3}{2} A^2 x - \frac{1}{2} B x - A y \right) + \delta^3 R(x, y),$$

where $\|R(x, y)\| \le C$. We note that in our case, $\delta = \frac{1}{L}$. This corresponds to a deterministic scaling regime, in contrast with the diffusive scaling $\delta = L^{-1/2}$ used in Theorem 1. Also, note that we are working on the product of the unit spheres. Thus, the remainder of order $1/L^2$ is bounded. This means that for every token number $i \le N$, we have

$$X_{n+1}^{i} - X_{n}^{i} = \frac{1}{L} P_{X_{n}^{i}}^{\perp} [A_{\beta}(X_{n}^{i}, \mathbf{X}_{n})] + \bar{o}(\frac{1}{L}),$$

where the remainder converges to zero uniformly for all $n \le L$. A standard Grönwall argument implies the uniform convergence to the solution of the following problem:

$$\begin{cases} \frac{dY_t^i}{dt} = P_{Y_t^i}^{\perp} [A_{\beta}(Y_t^i, \mathbf{Y}_t)], & \text{for } i \le N, \\ \mathbf{Y}_0 = \mathbf{X}_0. \end{cases}$$

Thus, we have obtained the required statement. $\square$

Proof of Lemma 3. We note that for the dynamics (7), we apply exactly the same technique as for Theorem 1. In this case, instead of the truncated dynamics (14), we obtain the following dynamics:

$$Y_{n+1}^{i} = Y_{n}^{i} + \frac{\epsilon v_{n+1}}{\sqrt{L}} P_{Y_{n}^{i}}^{\perp} [\tilde{A}_{\beta}(Y_{n}^{i}, \mathbf{Y}_{n})] + \frac{1}{L} [G(Y_{n}^{i}, \epsilon v_{n+1} \tilde{A}_{\beta}(Y_{n}^{i}, \mathbf{Y}_{n})) + P_{Y_{n}^{i}}^{\perp}[\tilde{A}_{\beta}(Y_{n}^{i}, \mathbf{Y}_{n})]], \quad i \in \{1, \dots, N\}.$$

Then, for the linearized dynamics above and the recurrence relation (7), we apply exactly the same argument. This leads to the result of the lemma. $\square$

## C. Long-Time Behavior and Clustering

**Outline of the Analysis.** This section analyzes the long-time behavior of the overlap dynamics underlying clustering. We first derive a closed one-dimensional description of the overlap process and classify its boundary points using the scale-speed framework of Karlin & Taylor (1981) (Lemmas 4, 5, and 6). Based on this classification, we then prove the clustering and anti-clustering results by exploiting the martingale properties of the associated scale function, making explicit a classical implication of one-dimensional diffusion theory in our setting.

### C.1. Overlap Dynamics and Reduction to One Dimension

We begin by observing that clustering events admit a simple geometric characterization in terms of pairwise inner products. Let $(Y_t^i)_{i=1}^N \subset \mathbb{S}^{d-1}$ be token trajectories solving the Transformer SDE (6). For any pair $(i, j)$, the squared distance and squared sum satisfy

$$\|Y_t^i - Y_t^j\|^2 = 2(1 - \langle Y_t^i, Y_t^j \rangle), \quad \|Y_t^i + Y_t^j\|^2 = 2(1 + \langle Y_t^i, Y_t^j \rangle),$$

since the tokens remain on the unit sphere. Consequently,

$$\|Y_t^i - Y_t^j\| \to 0 \iff \langle Y_t^i, Y_t^j \rangle \to 1,$$
which corresponds to clustering, while

$$ \|Y_t^i + Y_t^j\| \to 0 \iff \langle Y_t^i, Y_t^j \rangle \to -1, $$

corresponding to convergence toward an antipodal configuration. Thus, clustering and anti-clustering events are entirely characterized by the long-time behavior of the overlap process

$$ Z_t^{i,j} := \langle Y_t^i, Y_t^j \rangle \in [-1, 1]. $$

For general $N$, applying Itô's formula to $Z_t^{i,j}$ yields a one-dimensional stochastic differential equation whose coefficients depend on the remaining tokens through the attention weights. Hence the overlap dynamics are not closed in general. The situation becomes particularly tractable in the case $N = 2$, where the overlap forms a closed Markov diffusion.

**Lemma 4** (Closed overlap dynamics for two tokens: unnormalized vs. self-attention). Assume $N = 2$ and $K = Q = \text{Id}(d)$. Let $(Y_t^1, Y_t^2) \in (\mathbb{S}^{d-1})^2$ solve (6) and define the overlap

$$ Z_t := \langle Y_t^1, Y_t^2 \rangle \in [-1, 1]. $$

Then $(Z_t)_{t \ge 0}$ is a one-dimensional diffusion. More precisely:

(i) **Unnormalized self-attention** $A_\beta = A_\beta^U$. There exists a scalar Brownian motion $B_t$ such that

$$ dZ_t = \frac{\sigma^2}{4}(1 - Z_t^2) \left( 2(d-2)e^{\beta(1+Z_t)} - Z_t(e^{2\beta} + e^{2\beta Z_t}) \right) dt + \frac{\sigma}{\sqrt{2}}(1 - Z_t^2) \sqrt{e^{2\beta} + e^{2\beta Z_t}} dB_t. \quad (20) $$

(ii) **Normalized self-attention** $A_\beta = A_\beta^S$. There exists a scalar Brownian motion $\tilde{B}_t$ such that

$$ dZ_t = \sigma^2 \frac{(1 - Z_t^2)}{(e^\beta + e^{2\beta Z_t})^2} \left( 2(d-2)e^{\beta(1+Z_t)} - Z_t(e^{2\beta} + e^{2\beta Z_t}) \right) dt + \sigma \frac{(1 - Z_t^2)}{e^\beta + e^{2\beta Z_t}} \sqrt{e^{2\beta} + e^{2\beta Z_t}} d\tilde{B}_t. \quad (21) $$

*Proof.* We write the overlap as

$$ Z_t = \langle Y_t^1, Y_t^2 \rangle, $$

and apply Itô's formula to the bilinear map $(x, y) \mapsto \langle x, y \rangle$. We obtain the identity

$$ dZ_t = \langle dY_t^1, Y_t^2 \rangle + \langle Y_t^1, dY_t^2 \rangle + \langle dY_t^1, dY_t^2 \rangle. \quad (22) $$

We now substitute the dynamics of $(Y_t^1, Y_t^2)$ from (6). We recall that

$$ P_x^\perp [dW_t f(x)] = \sum_{k,l=1}^d c_{kl}(x) dW_t^{kl}, \quad c_{kl}(x) = f_l(x)(e_k - x_k x) \in \mathbb{R}^d, $$

where $e_k := [0 \dots \underbrace{1}_{k} \dots 0]^T$ and $f \in \mathbb{R}^d$. Expanding similarly expressions in (22) and using pair-wise independence of entries, we get:

$$ dZ_t = \sigma^2 \left( \frac{1}{2}(1-d)Z_t \left( \|A_\beta(Y_t^1, Y_t)\|^2 + \|A_\beta(Y_t^2, Y_t)\|^2 \right) + (Z_t^2 + d - 2) \langle A_\beta(Y_t^1, Y_t), A_\beta(Y_t^2, Y_t) \rangle \right) dt + \\ + \sigma \sum_{k,l=1}^d \left[ Y_t^{2,k} A_\beta^l(Y_t^1, Y_t) + Y_t^{1,k} A_\beta^l(Y_t^2, Y_t) - Z_t \left( Y_t^{1,k} A_\beta^l(Y_t^1, Y_t) + Y_t^{2,k} A_\beta^l(Y_t^2, Y_t) \right) \right] dW_t^{kl}, $$

where $A_\beta^k(Y, Y) := \langle A_\beta(Y, Y), e_k \rangle$ and $Y_t^{i,k} := \langle Y_t^i, e_k \rangle$. Normalizing the diffusion component, there exists a scalar Brownian motion $B_t \in \mathbb{R}$ (correlated with the original ones by the quadratic covariation property) such that:

$$ \sum_{k,l=1}^d \left[ Y_t^{2,k} A_\beta^l(Y_t^1, Y_t) + Y_t^{1,k} A_\beta^l(Y_t^2, Y_t) - Z_t \left( Y_t^{1,k} A_\beta^l(Y_t^1, Y_t) + Y_t^{2,k} A_\beta^l(Y_t^2, Y_t) \right) \right] dW_t^{kl} \\ = \sqrt{(1 - Z_t^2) \left( \|A_\beta(Y_t^1, Y_t)\|^2 + \|A_\beta(Y_t^2, Y_t)\|^2 - 2Z_t \langle A_\beta(Y_t^1, Y_t), A_\beta(Y_t^2, Y_t) \rangle \right)} dB_t. \quad (23) $$
Thus, the final expression for arbitrary two tokens is:

$$dZ_t = \sigma^2 \left( \frac{1}{2}(1-d)Z_t \left( \|\mathbf{A}_{\beta}(Y_t^1, \mathbf{Y}_t)\|^2 + \|\mathbf{A}_{\beta}(Y_t^2, \mathbf{Y}_t)\|^2 \right) + (Z_t^2 + d - 2) \langle \mathbf{A}_{\beta}(Y_t^1, \mathbf{Y}_t), \mathbf{A}_{\beta}(Y_t^2, \mathbf{Y}_t) \rangle \right) dt + \sigma \sqrt{(1-Z_t^2) \left( \|\mathbf{A}_{\beta}(Y_t^1, \mathbf{Y}_t)\|^2 + \|\mathbf{A}_{\beta}(Y_t^2, \mathbf{Y}_t)\|^2 - 2Z_t \langle \mathbf{A}_{\beta}(Y_t^1, \mathbf{Y}_t), \mathbf{A}_{\beta}(Y_t^2, \mathbf{Y}_t) \rangle \right)} dB_t.$$

For $N = 2$, we substitute the explicit formulas for unnormalized and normalized attention. For unnormalized attention $A_{\beta}^{U}$, we have

$$A_{\beta}^{U}(Y^{i}, \mathbf{Y}) = \frac{1}{2} (Y^{i} e^{\beta} + Y^{3-i} e^{\beta Z_t}), \quad i = 1, 2.$$

For normalized attention $A_{\beta}^{S}$, we have

$$A_{\beta}^{S}(Y^{i}, \mathbf{Y}) = \frac{Y^{i} e^{\beta} + Y^{3-i} e^{\beta Z_t}}{e^{\beta} + e^{\beta Z_t}}, \quad i = 1, 2.$$

Substituting these expressions and simplifying algebraically yields exactly (20) for case (i) and (21) for case (ii), as claimed. □

## C.2. Boundary Behavior of the Overlap Diffusion

By Lemma 4, the overlap process $(Z_t)_{t \ge 0}$ is a one-dimensional Itô diffusion on $[-1, 1]$ with smooth coefficients on $(-1, 1)$ that vanish only at the boundaries $\{-1, 1\}$; the diffusion coefficient is strictly positive in the interior. Hence, the long-time behavior is determined entirely by whether $Z_t$ can reach $\pm 1$ and, once there, remain there.

We employ Feller's boundary classification for one-dimensional diffusions (Karlin & Taylor, 1981, Chapter 15, §6, p.226). Both boundaries $\pm 1$ are inaccessible, but their attracting/non-attracting character depends on $(d, \beta)$. The analysis reveals that, for each $(d, \beta)$ and attention mechanism, the boundaries can become attracting. As a result, $Z_t$ converges to either $+1$ or $-1$ with positive probability, with no other possible limits.

**Lemma 5 (Boundary classification for two tokens, full noise).** Consider the overlap diffusion $Z_t$ for $N = 2$ tokens under the full noise model, for either unnormalized self-attention $A_{\beta}^{U}$ (SDE (20)) or normalized self-attention $A_{\beta}^{S}$ (SDE (21)). The boundary behavior is identical for both attention mechanisms and is given as follows:

1. For dimension $d = 2$, both boundaries $+1$ and $-1$ are attracting.

2. For dimensions $d \ge 3$:

(a) The boundary $+1$ is always attracting.

(b) The boundary $-1$ is attracting if and only if $d - 2 < \cosh(2\beta)$. Otherwise, it is non-attracting.

*Proof.* Both SDEs can be written in the standard form

$$dZ_t = \mu(Z_t)dt + \sigma(Z_t)dB_t,$$

with drift and diffusion coefficients that are smooth on $(-1, 1)$ and vanish at the endpoints. Crucially, the ratio $\frac{2\mu(z)}{\sigma^2(z)}$—which determines the scale function—is identical for the two attention mechanisms. Indeed,

$$\frac{2\mu(z)}{\sigma^2(z)} = \frac{2(d-2)e^{\beta(1+z)} - z(e^{2\beta} + e^{2\beta z})}{(1-z^2)(e^{2\beta} + e^{2\beta z})},$$

as can be verified directly from (20) and (21). Therefore, the derivative of the scale function (with base point $x_0 \in (-1, 1)$):

$$s(x) = \exp \left[ - \int_{x_0}^{x} \frac{2\mu(y)}{\sigma^2(y)} dy \right] = \exp \left[ - \int_{x_0}^{x} \frac{2(d-2)e^{\beta(1+y)} - y(e^{2\beta} + e^{2\beta y})}{(1-y^2)(e^{2\beta} + e^{2\beta y})} dy \right],$$

is the same for both dynamics.
Define $\psi(y) := \frac{2(d-2)e^{\beta(1+y)}}{e^{2\beta} + e^{2\beta y}}$. Then

$$ \mathbf{s}(x) = \begin{cases} \sqrt{\frac{1-x_0^2}{1-x^2}}, & d=2, \\ \sqrt{\frac{1-x_0^2}{1-x^2}} \exp\left[-\int_{x_0}^{x} \frac{\psi(y)}{1-y^2} dy\right], & d \ge 3. \end{cases} $$

For $d \ge 3$, near the right boundary we have $\psi(1) = d - 2$ and

$$ \int_{x_0}^{1-\varepsilon} \frac{\psi(y)}{1-y} dy = \frac{d-2}{2} |\log \varepsilon| + O(1), \quad \varepsilon \to 0^+. $$

Hence $s(1-\varepsilon) \asymp \varepsilon^{-1/2+\frac{d-2}{2}}$. The integral $\int_{x_0}^{1} s(x)dx$ diverges for all $d \ge 2$, so the boundary +1 is always attracting.

Near the left boundary, $\psi(-1) = \frac{d-2}{\cosh(2\beta)}$ and

$$ \int_{-1+\varepsilon}^{x_0} \frac{\psi(y)}{1+y} dy = \frac{d-2}{2 \cosh(2\beta)} |\log \varepsilon| + O(1), \quad \varepsilon \to 0^+. $$

Thus $s(-1+\varepsilon) \asymp \varepsilon^{-1/2+\frac{d-2}{2\cosh(2\beta)}}$. The boundary $-1$ is attracting precisely when $\int_{-1}^{x_0} s(x)dx < \infty$, which occurs iff $-\frac{1}{2} - \frac{d-2}{2\cosh(2\beta)} > -1$, i.e. $d-2 < \cosh(2\beta)$.

For $d=2$, the explicit scale function yields $s(1-\varepsilon) \asymp \varepsilon^{-1/2}$ and $s(-1+\varepsilon) \asymp \varepsilon^{-1/2}$, so both boundaries are attracting. $\square$

**Lemma 6 (Attainability).** For the diffusion processes defined in (20), (21), (26), with $Z_0 \in (-1, 1)$, we have the following properties:

1. $|Z_t| \le 1$ a.s. for all $t \ge 0$;

2. Both boundaries $\pm 1$ are inaccessible; that is, $|Z_t| < 1$ for any finite $t$.

*Proof.* The first item follows directly from the definition $Z_t = \langle X_t^1, X_t^2 \rangle$ for two unit vectors $X_t^1, X_t^2$ and the Cauchy–Schwarz inequality. For the second item, we could verify it directly by computing the corresponding integral $\Sigma$ in Karlin & Taylor (1981), but instead we consider the function $f(x) = \tanh^{-1}(x)$, defined as

$$ f(x) := \frac{1}{2} \ln \left( \frac{1+x}{1-x} \right). $$

For the generic dynamics

$$ dZ_t = (1 - Z_t^2)(c(Z_t) dt + \gamma(Z_t) dB_t), $$

which correspond to (20), (21), (26), where $c$ and $\gamma$ are smooth and bounded on $[-1, 1]$, we consider $Y_t = f(Z_t)$. Applying Itô's formula, we obtain:

$$ dY_t = (c \circ \tanh(Y_t) + Z_t \cdot (\gamma \circ \tanh(Y_t))^2) dt + \gamma \circ \tanh(Y_t) dB_t, \quad Y_0 = f(Z_0). $$

By Itô's Uniqueness and Existence Theorem, $Y_t$ is well-defined and does not explode to $\pm\infty$ in finite time. Thus, $Z_t = \tanh(Y_t)$ never reaches $\pm 1$. $\square$
### C.3. Proof of the Phase Transition Result (Theorem 2)

With the boundary classification of the overlap diffusion established in Lemmas 5,6, we are now equipped to prove the main result of this section. Before doing so, we introduce some additional notation. Based on the function $\mathbf{s}(x)$, we define a scale function:

$$ \mathbf{S}(x) = \int_{x_0}^{x} \mathbf{s}(y) dy. $$

The main property is that directly applying Itô's formula, we get that $\mathbf{S}(Z_t)$ is a local martingale. Thanks to the classification of the boundary points, we know precisely the range of parameters $(\beta, d)$ for which $\mathbf{S}(\pm 1)$, defined as the corresponding limit, is finite, and vice versa. Note that this behavior does not depend on the choice of the initial point $x_0 \in (-1, 1)$. We then define a convenient set of stopping times:

$$ \tau_a := \inf \{t \ge 0 \mid Z_t = a\} \text{ for } a \in (-1, 1). $$

Although Theorem 2 follows directly from the boundary classification above, the limiting behavior itself requires a proof. While this result is folklore and appears in various forms across different textbooks, we have struggled to find it stated in exactly the form we need. For the sake of clarity, we therefore include a proof based on unpublished lecture notes by Jacobsen (2008).

*Proof of Theorem 2.* For the case when two initial token positions do not coincide, we have $Z_0 = z_0 \in (-1, 1)$. We restrict ourselves to the case when $z_0$ is non-random, but the results remain the same for random initial distributions by convolution. We note that both boundaries $\pm 1$ are unattainable in finite time by Lemma 6.

The proof relies on the properties of the function $\mathbf{S}(x)$ defined above. If $\mathbf{S}(1) < \infty$, then

$$ \lim_{t \to \infty} Z_t = 1 \quad \mathbb{P}^{z_0}\text{-almost surely on } A_-, \quad \text{where } A_- = \bigcup_{a: a < z_0} (\tau_a = \infty), \quad \mathbb{P}^{z_0}(A_-) = \frac{\mathbf{S}(z_0) - \mathbf{S}(-1)}{\mathbf{S}(1) - \mathbf{S}(-1)}. \quad (24) $$

Moreover, if $\mathbf{S}(-1) > -\infty$, then

$$ \lim_{t \to \infty} Z_t = -1 \quad \mathbb{P}^{z_0}\text{-almost surely on } A_+, \quad \text{where } A_+ = \bigcup_{b: b > z_0} (\tau_b = \infty), \quad \mathbb{P}^{z_0}(A_+) = \frac{\mathbf{S}(1) - \mathbf{S}(z_0)}{\mathbf{S}(1) - \mathbf{S}(-1)}. $$

Indeed, let's show that for $\mathbf{S}(1) < \infty$, we have (24). Clearly, for $a < z_0 \{\tau_a = \infty\} \uparrow A_-$ as $a \to -1$, then

$$ \mathbb{P}^{z_0}(A_-) = \lim_{a \downarrow -1} \mathbb{P}^{z_0}(\tau_a = \infty). $$

Note that $\mathbf{S}(Z_{t \wedge \tau_a}) := \mathbf{S}(Z_{\min\{t, \tau_a\}})$ is a bounded local martingale, hence a true martingale. Boundedness follows from the fact that $\mathbf{S}(Z_{t \wedge \tau_a}) \in [a, \mathbf{S}(1)]$ and we assumed that $\mathbf{S}(1) < \infty$. Then, by the martingale convergence theorem, we can well define $\mathbf{S}(Z_{\tau_a})$ as follows:

$$ \mathbf{S}(Z_{\tau_a}) := \begin{cases} \mathbf{S}(a) & \text{on } (\tau_a < \infty), \\ \lim_{t \to \infty} \mathbf{S}(Z_t) & \text{on } \tau_a = \infty. \end{cases} $$

So, we have from one side

$$ \mathbb{E}^{z_0}[\mathbf{S}(Z_{\tau_a})] = \mathbf{S}(z_0), \quad (25) $$

from another, by inaccessibility of 1

$$ \mathbb{P}^{z_0}[\tau_a = \infty] = \lim_{z_0 < b} \mathbb{P}^{z_0}[\tau_a > \tau_b] = \lim_{b \to 1} \frac{\mathbf{S}(z_0) - \mathbf{S}(a)}{\mathbf{S}(b) - \mathbf{S}(a)} = \frac{\mathbf{S}(z_0) - \mathbf{S}(a)}{\mathbf{S}(1) - \mathbf{S}(a)}, $$

Substituting it to (25), we get

$$ \mathbf{S}(z_0) = \mathbb{E}^{z_0}[\mathbf{S}(Z_{\tau_a})] = \mathbf{S}(a) \cdot (1 - \mathbb{P}^{z_0}[\tau_a = \infty]) + \mathbb{E}^{z_0}[\mathbf{S}(Z_{\tau_a}) \mathbf{Ind}\{\tau_a = \infty\}], $$

and thus simplifying we get that

$$ \mathbb{E}^{z_0}[\mathbf{S}(Z_{\tau_a}) \mathbf{Ind}\{\tau_a = \infty\}] = \mathbf{S}(1)\mathbb{P}(\tau_a = \infty). $$

Since $\mathbf{S}(Z_{\tau_a}) \le \mathbf{S}(1)$, we have that $\mathbf{S}(Z_{\tau_a}) = \mathbf{S}(1) \mathbb{P}^{z_0}$ almost surely on $\{\tau_a = \infty\}$. Then, (24) is proved. From this and the boundary classification in Lemma 5, we get the statement of the theorem. $\square$
### C.4. Hybrid Model Boundary Behavior

We now turn to the hybrid clustering regime.

*Proof of Lemma 1.* The overlap equation similar to (21) and (20) is given as follows:

$$dZ_t = e^{\beta Z_t} (Z_t^2 - 1) \left[ \left( \frac{\varepsilon^2 e^{\beta}}{2} + Z_t \varepsilon^2 e^{\beta Z_t} - 1 \right) dt - \varepsilon dB_t \right]. \quad (26)$$

The reasoning about convergence is exactly the same as in the proof of Theorem 2, since Lemma 6 applies to both. Here we only perform the boundary classification.

In this case, for $s(x)$ we obtain the following expression:

$$s(x) = \exp \left( -2 \int_{x_0}^{x} \frac{\varepsilon^2 e^{\beta(1-y)}}{\varepsilon^2(y^2-1)} dy \right) = \exp \left( \int_{x_0}^{x} \frac{(e^{\beta} - \frac{2}{\varepsilon^2})e^{-\beta y} + y}{1-y^2} dy \right).$$

Simplifying the integral, we get:

$$s(x) = \frac{1-x_0^2}{1-x^2} \exp \left( \left( e^{\beta} - \frac{2}{\varepsilon^2} \right) \int_{x_0}^{x} \frac{e^{-\beta y}}{1-y^2} dy \right).$$

Hence,

$$s(1-\delta) \sim \frac{1}{\delta} \cdot \exp \left( -\frac{1}{2} \left( e^{\beta} - \frac{2}{\varepsilon^2} \right) e^{-\beta} \log(\delta) \right), \quad s(-1+\delta) \sim \frac{1}{\delta} \cdot \exp \left( \frac{1}{2} \left( e^{\beta} - \frac{2}{\varepsilon^2} \right) e^{\beta} \log(\delta) \right), \text{ as } \delta \to 0.$$

Simplifying, we get:

$$s(1-\delta) \sim \delta^{-\frac{3}{2}+\frac{e^{-\beta}}{\varepsilon^2}}, \quad s(-1+\delta) \sim \delta^{-1+\frac{e^{2\beta}}{2}-\frac{e^{\beta}}{\varepsilon^2}}.$$

From the standard integrability conditions, we have

1. $\frac{1}{\varepsilon^2} > \frac{e^{\beta}}{2}$ — the boundary $s = 1$ is attracting, while $s = -1$ is non-attracting.

2. $\frac{1}{\varepsilon^2} < \frac{e^{\beta}}{2}$ — the boundary $s = 1$ is non-attracting, while $s = -1$ is attracting.

This completes the proof. □

## D. Additional Experimental Results

This section provides additional numerical experiments demonstrating the convergence dynamics of token clustering. We validate that the chosen time horizons are sufficiently large for clustering behavior to stabilize, and demonstrate the robustness of our theoretical findings across different parameter regimes.

### D.1. Phase Transition Analysis Across Parameter Regimes

To verify that our choice of time horizon in the main paper (Figure 3) is conservative, we present convergence curves at shorter time horizons ($T = 50$ and $T = 100$ versus $T = 500$ on Figure 3). Specifically, we track two key quantities:

1. The proportion of tokens that fail to cluster into any configuration

2. The proportion of observed clustering configurations, which stabilizes well below the maximum time horizon
We examine three representative values of the inverse temperature parameter $\beta$: the regimes ($\beta \to 0$ and $\beta \to \infty$) where behavioral differences are most pronounced, and an intermediate regime. These regimes are chosen to demonstrate that qualitative features persist across the parameter space. All other experimental settings are identical to those used in the main phase diagram (Figure 3), with the sole difference being the reduced time horizon. Reported curves are obtained by averaging over $2 \cdot 10^3$ independent trajectories.

We focus on three representative regimes of the inverse temperature parameter $\beta$, corresponding to qualitatively distinct clustering behaviors. In the first regime, two-token clustering is possible for all dimensions. In the second regime, such clustering occurs only for a subset of dimensions. In the third regime, two-token clustering again becomes possible uniformly across all dimensions. The selected values of $\beta$ are chosen to illustrate these three regimes.

<figure><img src="image_10.png" /><figcaption>(a) Below-threshold regime ($\beta < 1$): convergence toward antipodal clustering does not happen for all dimensions.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) Near-threshold regime ($\beta = 1$). Convergence to the antipodal configuration happens only for $d < 6$.</figcaption></figure>

**Figure 8.** Convergence behavior across representative values of the inverse temperature $\beta$. Panels (a)–(c) correspond respectively to regimes in which all trajectories converge to a single-cluster configuration, in which convergence to an antipodal configuration is attainable for all dimensions considered, and in which the antipodal configuration is attainable only for certain dimensions while remaining inaccessible for others.

## D.2. Rate of Convergence and Metastable behavior

We provide additional numerical illustrations supporting the qualitative observations on convergence rates and metastable behavior discussed in the main text. In particular, we focus on two representative regimes of the inverse temperature
parameter $\beta$, corresponding to smaller and larger values.

**Small-$\beta$ regime (comparable convergence rates).** In this regime, both the stochastic and deterministic dynamics converge predominantly to a single-cluster configuration. Numerically, the observed convergence rates of the two dynamics are comparable. The deterministic dynamics is given by the specialization of (7) with $\varepsilon = 0$, $A_\beta = A_\beta^S$, and $Q^K K = \text{Id}(4)$. See Figure 9a.

**Large-$\beta$ regime (absence of metastability in the stochastic dynamics).** For larger values of $\beta$, both dynamics still converge to structured clustered configurations; however, their qualitative convergence behavior differs markedly. The deterministic dynamics exhibits metastable behavior, characterized by long-lived transient configurations before convergence. In contrast, no such metastability is observed numerically in the stochastic dynamics, which converges significantly faster. See Figure 9b.

Overall, across the regimes considered here, the stochastic dynamics does not display metastable behavior in the sense of Geshkovski et al. (2024), in contrast to its deterministic counterpart at large values of $\beta$.

<figure><img src="image_12.png" /><figcaption>Figure 9. Comparison of convergence behavior in two representative regimes. **Left:** $\beta = 0.5$, number of tokens $N = 50$, time horizon $T = 15$, with standard discretization ($L = 100$). **Right:** $N = 50$, $\beta = 4.5$, time horizon $T = 30$, with the same discretization. In both cases, an antipodal loop configuration is detected when all pairwise inner products satisfy $|\langle x_i, x_j \rangle| \ge 1 - \varepsilon$ and at least one pair satisfies $\langle x_i, x_j \rangle \le -1 + \varepsilon$, with $\varepsilon = 10^{-3}$.</figcaption></figure>