How Do Transformers Model Physics? Investigating the Simple Harmonic Oscillator

Abstract

How do transformers model physics? Do transformers model systems with interpretable analytical solutions or do they create an “alien physics” that is difficult for humans to decipher? We have taken a step towards demystifying this larger puzzle by investigating the simple harmonic oscillator (SHO), $\ddot{x}+2\gamma \dot{x}+{\omega }_0^{2}x=0$, one of the most fundamental systems in physics. Our goal was to identify the methods transformers use to model the SHO, and to do so we hypothesized and evaluated possible methods by analyzing the encoding of these methods’ intermediates. We developed four criteria for the use of a method within the simple test bed of linear regression, where our method was $y=wx$ and our intermediate was w: (1) Can the intermediate be predicted from hidden states? (2) Is the intermediate’s encoding quality correlated with the model performance? (3) Can the majority of variance in hidden states be explained by the intermediate? (4) Can we intervene on hidden states to produce predictable outcomes? Armed with these two correlational (1,2), weak causal (3), and strong causal (4) criteria, we determined that transformers use known numerical methods to model the trajectories of the simple harmonic oscillator, specifically, the matrix exponential method. Our analysis framework can conveniently extend to high-dimensional linear systems and nonlinear systems, which we hope will help reveal the “world model” hidden in transformers.

1. Introduction

Transformers are state-of-the-art models for a range of tasks [1,2,3,4], but our understanding of how these models represent the world is limited. Recent work in mechanistic interpretability [5,6,7,8,9,10,11,12] has shed light on how transformers represent mathematical tasks like modular addition [7,13,14], yet little work has been done to understand how transformers model physics. This question is crucial, as for transformers to build any sort of “world model” they must have a grasp of the physical laws that govern the world [15].

Our key research question was the following: How do transformers model physics? This question is intimidating, since even humans have many different ways of modeling the same underlying physics [16]. In the spirit of hypothesis testing, we reformulated the question as follows: Given a known modeling method g, does the transformer learn g? If a transformer leverages g, its hidden states must encode information about important intermediate quantities in g. We focused our study on the simple harmonic oscillator $\ddot{x}+2\gamma \dot{x}+{\omega }_0^{2}x=0$, where $\gamma$ and ${\omega }_0$ are the damping and frequency of the system, respectively. Given the trajectory points $\{(x_0,v_0),(x_1,v_1),\dots ,(x_n,v_n)\}$ at discrete times $\{t_0,t_1,\dots ,t_n\}$, we tasked a transformer with predicting $(x_{n+1},v_{n+1})$ at time $t_{n+1}$, as shown in Figure 1. In this setting, g could be a numerical simulation the transformer runs after inferring $\gamma ,{\omega }_0$ from past data points. We would then expect some form of $\gamma$ and ${\omega }_0$ to be intermediates encoded in the transformer. How, we asked, could we show that intermediates and the method g were being used?

We developed criteria to demonstrate that the transformer was using g by studying intermediates in a simpler setting: in-context linear regression, $y=wx$. As correlational evidence for the model’s internal use of w, we found that the intermediate w can be encoded linearly, nonlinearly, or not at all. We also linked the performance of models to their encoding of w, and we used this as an explanation for in-context learning. We generated causal evidence for the use of w by analyzing how much of the hidden states’ variance was explained by w and by linearly intervening in the network, to predictably change its behavior.

We used these developed criteria of intermediates to study how transformers model the simple harmonic oscillator (SHO), a fundamental model in physics. We generated multiple hypotheses for the method(s) that transformers use to model the trajectories of SHOs, and we used our criteria from linear regression to show correlational and causal evidence that transformers employ known numerical methods, specifically the matrix exponential, to model trajectories of SHOs. Although our analysis was constrained to the SHO in this paper, our framework naturally extends to some high-dimensional linear and nonlinear systems.

The organization of this paper is as follows. In Section 2, we overview related work. In Section 3, we define and investigate intermediates in the setting of linear regression, and we use this to develop criteria for transformers’ use of a method g. In Section 4, we hypothesize that transformers use numerical methods to model the SHO, and we use our criteria of intermediates to provide causal and correlational evidence for transformers’ use of the matrix exponential.

2. Related Work

Mechanistic interpretability. Mechanistic interpretability (MI), as a field, aims to understand the specific computational procedures machine learning models use to process inputs and produce outputs [5,6,7,8,9,10,11,12,13]. Some MI work focuses on decoding the purpose of individual neurons [17], while other work focuses on ensembles of neurons [11,12]. Our work is aligned with the latter.

Algorithmic behaviors in networks. A subset of MI attempts to discover the specific algorithms that networks use to solve tasks by reverse engineering weights. For example, it has been demonstrated that transformers use the discrete Fourier transform to model modular addition [13]. Reverse engineering has been studied extensively for neural circuits, which has, impressively, connected the behavior of neurons and brain areas [18,19]. We focused on transformers, and, instead of reverse engineering weights, we made use of linear probing [20] to discover byproducts of algorithms represented internally by transformers. Studies have found that algorithms in models are potentially an “emergent” behavior that manifests with size [21,22], which we also found.

AI and Physics. Many works design specialized machine learning architectures for physics tasks [23,24,25,26,27], but less work has been undertaken to see how well transformers perform on physical data out of the box. Recently, it was shown that LLMs can in-context learn physics data [15], which inspired the research question of this paper: how do transformers model physics?

3. Developing Criteria for Intermediates with Linear Regression

Our main goal was to determine which methods transformers use to model the simple harmonic oscillator. We aimed to do this by generating criteria based on the encoding of relevant intermediates. For this section, we developed our criteria of intermediates in a simpler setting: linear regression. Notably, linear regression is identical to predicting the acceleration from the position of an undamped harmonic oscillator ($\gamma =0$), making this setup physically relevant.

Setup. In our linear regression setup, we generated $\mathit{X}$ and $\mathit{w}$ between $[-0.75,0.75]$, where $\mathit{X}$ had size (5000,65) and $\mathit{w}$ had size (5000). We generated $\mathit{Y}=\mathit{wX}$, and we trained the transformers to predict $y_{n+1}$ given $\{x_1,y_1,...,x_n,y_n,x_{n+1}\}$.

Since in-context linear regression is well studied for transformers [28,29], we used this simple setting to ask and answer fundamental questions about intermediates, namely:

• What is an intermediate?
• How can intermediates be encoded and how can we robustly probe for them?
• When, or under what circumstances, are intermediates encoded?

All of these questions developed an understanding of intermediates that built up to the key question: How can we use intermediates to demonstrate that a transformer is actually using a method in its computations? By answering this question for linear regression, we generated four correlational and causal criteria to demonstrate that a transformer is using a method in its computations, which we could then apply to understanding the simple harmonic oscillator, as shown in Figure 1.

3.1. What Is an Intermediate?

We define an intermediate as a quantity that a transformer uses during computation, but which is not a direct input/output to/of the transformer. More formally, if the input to the transformer is $\mathit{X}$ and its output is $\mathit{Y}$, we can model the transformer’s computation as $\mathit{Y}=g(\mathit{X},I)$, where g is the method used and I is the intermediate of that method. For example, if we want to determine if the transformer is computing the linear regression task using $\mathit{Y}=\mathit{wX}$ then $I=\mathit{w}$, $g(\mathit{X},I)=g(\mathit{X},\mathit{w})=\mathit{wX}$.

3.2. How Can Intermediates Be Encoded and How Can We Robustly Probe for Them?

We wanted to understand what form of the intermediate, $f\left(I\right)$, is encoded in the network’s hidden states. For example, while it may be obvious to humans to compute $y=wx$, perhaps transformers prefer $exp(log(w)+log(x\left)\right)$ or $\sqrt{w^2{x}^2}$. We wanted to develop a robust probing methodology that captures these diverse possibilities. We identified three ways an intermediate I can be represented: linearly encoded, nonlinearly encoded, and not encoded at all. We used $HS$ to mean hidden state.

Linearly encoded. We say I is linearly encoded in a hidden state $HS$ if there is a linear network that takes $I=\mathrm{Linear}\left(HS\right)$. We determine the strength of the linear encoding by evaluating how much of the variance in I can be explained by $HS$, i.e., the $R^2$ of the probe.

Nonlinearly encoded. To probe for an arbitrary $f\left(I\right)$, we define a novel Taylor probe, which finds coefficients $a_i$, such that $f\left(I\right)=a_{1}I+a_2{I}^2+...+a_n{I}^n$, and $f\left(I\right)=\mathrm{Linear}\left(HS\right)$.

To actually implement this probing style, we use Canonical Correlation Analysis probes, which, given some multivariate data X and Y, find directions within X and Y that are maximally correlated [30]. Here, $X=[I,I^2,I^3,...,I^n]$ and $Y=HS$. If I is of bounded magnitude and n is sufficiently large, we are able to probe the transformer for any function $f\left(I\right)$. In practice, we used $n\le 5$.

Not encoded. If I fails to be linearly or nonlinearly encoded, we say that it is not encoded within the network. For example, there are at least two ways to predict $y_2$ from $\{x_1,y_1,x_2\}$, such that $y_2={\displaystyle \frac{y_1}{x_1}}{x}_2$: (1) $w=y_1/x_1$ is encoded, and $y_2=w{x}_2$; (2) $w^{\mathbf{\prime }}=x_2/x_1$ is encoded (so $w=y_1/x_1$ is not encoded) and $y_2=w^{\mathbf{\prime }}{y}_1$. Thus, it is not guaranteed that w is encoded.

3.3. When, or Under What Circumstances, Are Intermediates Encoded?

We wanted to apply our probing techniques to better understanding what type of models generate intermediates. Under the described setting of linear regression, we trained GPT style transformers of size $L=[1,2,3,4,5]$ and $H=[2,4,8,16,32]$, where L was the number of layers and H was the hidden size of the transformer. All the transformers trained in this study used one attention head and no LayerNorm to aid interpretability, and they were trained on a NVIDIA Volta GPU with the hyperparameters $\mathrm{epochs}=20,000, \mathrm{lr}={10}^{-3}, \mathrm{batchsize}=64$, using the Adam optimizer [31]. We did not use a token embedding—instead, the inputs were raw numbers representing the position of the SHO. We found that these models generalized to the out-of-distribution test data ($0.75\le \left|\mathit{w}\right|\le 1$) in Appendix A Figure A1, but we focused on investigating the intermediates on in-distribution training data.

Larger models have stronger encodings of intermediates. We found that the smaller models often did not have $\mathit{w}$ encoded, while the larger models encoded $\mathit{w}$ linearly, as evidenced by Figure 2. We formalized this further by defining $max\left(\overline{R^2}\right)$ as the maximum value taken over the depth positions of the mean $R^2$ of the $\mathit{w}$ probes taken over the context length. As shown in Appendix A Figure A2, we observed a clear phase transition in encoding across model size, and we also found that $max\left(\overline{R^2}\right)$ did not significantly improve if we extended the degree of the Taylor probes to $n>2$. Thus, in the case of linear regression, we found that the models represented $\mathit{w}$ linearly, quadratically, or not at all.

We attributed the stronger encoding of $\mathit{w}$ in the larger models to the “lottery ticket hypothesis”—larger models have more “lottery tickets” in their increased capacity to find a “winning” representation of $\mathit{w}$ [32,33]. Interestingly, the intuitive understanding that larger models have $\mathit{w}$ better encoded led us to the counterintuitive conclusion that larger models are actually more interpretable for our purposes.

Encoding quality is tied to model performance. As shown in Appendix A Figure A3, we found that the better-performing models generally had stronger encodings of $\mathit{w}$. As shown in Figure 3, we also found that the improvements in model prediction as a function of context length, or in-context learning, were correlated to improvements in $\mathit{w}$’s encoding, which we would expect if our models were using $\mathit{w}$ in their computations.

3.4. Key Question: How Can We Use Intermediates to Demonstrate That a Transformer Is Actually Using a Method in Its Computations?

So far, we had discovered that models encode $\mathit{w}$, either linearly or nonlinearly, and we had found relationships between model size, performance, and encoding strength. But how, we asked, could we ensure that the model was actually using $\mathit{w}$ in its computations and that the encoding of $\mathit{w}$ was not just a meaningless byproduct [34]?

Reverse Probing. To ensure that $\mathit{w}$ was not encoded in some insignificant part of the residual stream, we set up probes going from $[w,w^2]\to HS$, as opposed to $HS\to f\left(w\right)$. As shown in Figure 4, we often found that $\mathit{w}$ could explain large amounts of variance in model hidden states, implying that these hidden states were dedicated to representing $\mathit{w}$. We took this as weak causal evidence that $\mathit{w}$ was being used by the model—otherwise, it was unclear why a part of the model would be dedicated to storing $\mathit{w}$.

Intervening. We could also use reverse probes to intervene in the models’ hidden states and predictably change their output from $\mathit{w}\to {\mathit{w}}^{\mathbf{\prime }}$. As shown in Figure 4, we attempted to make ${\mathit{w}}^{\mathbf{\prime }}=0.5$ for all series and then measure the observed $\widehat{\mathit{w}}$ from the models’ outputs ($\widehat{w}={\widehat{y}}_n/x_n$). For 4 out of 25 models the intervention worked, providing strong causal evidence that the model uses its internal representation of $\mathit{w}$ in computations. For models where we identified a quadratic representation of $\mathit{w}$, we see that $\mathit{w}=0.5,-0.5$ were both represented in the observed intervention.

Putting it all together. We were able to generalize our understanding of intermediates from linear regression, to create criteria for a transformer’s use of a method g in its computations.

• Criteria for use of a method g with an associated unique intermediate I:

• If a model uses a method g then its hidden states should encode I (shown in Figure 2).
• If a model uses a method g then the model performance should improve if I is better represented (shown in Figure 3).
• If and only if the model uses g then we expect some hidden state’s variance to be almost fully explained by I (shown in Figure 4).
• If and only if the model uses g then we can intervene with hidden states, to change $I\to I^{\prime }$ and predictably change the model output from $g(\mathit{X},I)\to g(\mathit{X},I^{\prime })$ (shown in Figure 4).

The first two criteria for a transformer’s use of g are correlational, and the last two are weak and strong causal. Using these criteria (summarized in Figure 1), we could then investigate how transformers model more complex systems like the simple harmonic oscillator.

4. Investigating the Simple Harmonic Oscillator

We next applied our developed criteria of intermediates to investigating how transformers represent physics, specifically the methods they use to model the simple harmonic oscillator (SHO). The simple harmonic oscillator is ubiquitous in physics: it is used to describe phenomena as diverse as the swing of a pendulum, molecular vibrations, the behavior of AC circuits, and quantum states of trapped particles. Given a series of position and velocity data for a simple harmonic oscillator at a sequence of timesteps, we asked

• Can a transformer successfully predict the position/velocity at the SHO’s next timestep?
• Can we determine what computational method the transformer is using in this prediction?

4.1. Mathematical and Computational Setup

The simple harmonic oscillator is governed by the linear ordinary differential equation (ODE):

$$\ddot{x}+2\gamma \dot{x}+{\omega }_0^{2}x=0.$$

(1)

The two physical parameters of this equation are $\gamma$, the damping coefficient, and ${\omega }_0$, the natural frequency of the system. An intuitive picture for the SHO is a mass on a spring that is pulled from its equilibrium position by some amount $x_0$ and let go, as visualized in Figure 1; ${\omega }_0$ is related to how fast the system oscillates, and $\gamma$ is related to how soon the system decays to equilibrium from the internal resistance of the spring. We focused on studying how a transformer modeled the undamped harmonic oscillator, where $\gamma =0$. Given some initial starting position ($x_0$), velocity ($v_0$), and timestep $\mathsf{\Delta }t$, the time evolution of the undamped harmonic oscillator was

$$\begin{array}{cc}\hfill x_k& =x_{0}cos\left(k{\omega }_0\mathsf{\Delta }t\right)+{\displaystyle \frac{v_0}{{\omega }_0}}sin\left(k{\omega }_0\mathsf{\Delta }t\right)\hfill \\ \hfill v_k& =v_{0}cos\left(k{\omega }_0\mathsf{\Delta }t\right)-{\omega }_0{x}_{0}sin\left(k{\omega }_0\mathsf{\Delta }t\right),\hfill \end{array}$$

(2)

where $v={\displaystyle \frac{dx}{dt}}$. We generated 5000 timeseries of 65 timesteps for various values of ${\omega }_0,\mathsf{\Delta }t,x_0,$ and $v_0$, described in Appendix C. Following the procedure for linear regression, we trained transformers of size $L=[1,2,3,4,5]$ and $H=[2,4,8,16,32]$ to predict $(x_{n+1},v_{n+1})$, given $\{(x_0,v_0),(x_1,v_1),\dots (x_n,v_n)\}$. In Appendix C Figure A4, we see that our transformers were able to accurately predict the next timestep in the timeseries of out-of-distribution test data, and this prediction became more accurate with context length (i.e., in-context learning). But how was the transformer modeling the simple harmonic oscillator internally?

4.2. What Methods Could the Transformer Use to Model the Simple Harmonic Oscillator?

Human physicists would model the simple harmonic oscillator with the analytical solution to Equation (1), but it is unlikely that a transformer does so. Transformers are numerical approximators that use statistical patterns in data to make predictions, and, in that spirit, we hypothesize that transformers use numerical methods to model SHOs. There is a rich literature on numerical methods that approximate solutions to linear ordinary differential equations [35,36,37], and we highlight three possible methods the transformer could be using in our theory hub. For notation, we note that Equation (1) can be written as

$$\left[\begin{array}{c}\dot{x}\\ \dot{v}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -{\omega }_0^2& -2\gamma \end{array}\right]\left[\begin{array}{c}x\\ v\end{array}\right]=\mathit{A}\left[\begin{array}{c}x\\ v\end{array}\right].$$

(3)

Linear Multistep Method. Our model could have been using a linear multistep method, which uses values of derivatives from several previous timesteps to estimate the future timestep. We describe the kth order linear multistep method in

Table 1. Our theory hub of numerical methods and relevant intermediates that transformers could be using to model the simple harmonic oscillator.

Method	$\mathit{g}(\mathit{X},\mathit{I})$	I
Linear Multistep	${\displaystyle \left[\begin{array}{c}{x}_{n+1}\\ v_{n+1}\end{array}\right]=\sum _{j=0}^k{\alpha }_j\left[\begin{array}{c}{x}_{n-j}\\ v_{n-j}\end{array}\right]+\sum _{j=-1}^k{\beta }_j\mathit{A}\mathsf{\Delta }t\left[\begin{array}{c}{x}_{n-j}\\ v_{n-j}\end{array}\right]}$	$\mathit{A}\mathsf{\Delta }t$
Taylor Expansion	${\displaystyle \left[\begin{array}{c}{x}_{n+1}\\ v_{n+1}\end{array}\right]=\sum _{j=0}^k{\mathit{A}}^j\frac{\mathsf{\Delta }{t}^j}{j!}\left[\begin{array}{c}{x}_n\\ v_n\end{array}\right]}$	${\left(\mathit{A}\mathsf{\Delta }t\right)}^j$
Matrix Exponential	${\displaystyle \left[\begin{array}{c}{x}_{n+1}\\ v_{n+1}\end{array}\right]=e^{\mathit{A}\mathsf{\Delta }t}\left[\begin{array}{c}{x}_n\\ v_n\end{array}\right]}$	$e^{\mathit{A}\mathsf{\Delta }t}$

with coefficients ${\alpha }_j$ and ${\beta }_j$.

Taylor Expansion Method. The model could also have been using higher-order derivatives from the previous timestep to predict the next timestep (this is equivalent to the nonlinear single-step Runge–Kutta method for a homogeneous linear ODE with constant coefficients). We describe the kth order Taylor expansion in Table 1.

Matrix Exponential Method. While the two methods presented above are useful approximations for small $\mathsf{\Delta }t$, the matrix exponential uses a $2\times 2$ matrix to exactly transform the previous timestep to the next timestep. We describe this in Table 1. This method is the ${lim}_{k\to \infty }$ of the Taylor expansion method.

In order to use the criteria described in Section 3 to figure out which method(s) our model was using, we needed to define the relevant intermediates for each method g. Similarly to the linear regression, the intermediates were the coefficients of the input, but were now $2\times 2$ matrices and not a single value. We summarize our methods and intermediates in our theory hub in Table 1. Notably, these methods are viable for any homogeneous linear ordinary differential equation with constant coefficients and potentially for nonlinear differential equations as well (see Appendix B).

4.3. Evaluating Methods for the Undamped Harmonic Oscillator

We applied the four criteria established for linear regression (Figure 1) to evaluating if transformers use the methods in Table 1. For the Taylor expansion intermediate, we used $j=3$ to distinguish it from the linear multistep method, although our results were generally robust for $j\le 5$ (Appendix C Figure A6). We summarize our evaluations across the methods and criteria in Table 2.

Criterion 1: Is the intermediate encoded? In Figure 5, we see that all three intermediates were well encoded in the model, with the matrix exponential method especially prominent. This provides initial correlational evidence that the models were learning numerical methods. The magnitude of the encodings was generally smaller than the linear regression case, which we attribute to the increased difficulty of encoding $2\times 2$ matrices compared to a single weight value $\mathit{w}$. Notably, we only probed for linear encodings, given that $\mathit{w}$ was most often encoded linearly in the linear regression case study:

Criterion 2: Is the intermediate encoding correlated with the model performance? In Figure 6, we see that for all three methods the better-performing models generally had stronger encodings and the worse-performing models had weaker encodings. This correlation was strongest for the matrix exponential method. This provided more correlational evidence that our models were using the described methods.

Criterion 3: Can the intermediates explain the models’ hidden states? As shown in Figure 7, we reverse probed from the intermediates to the models’ hidden states, and we found that all the methods explained non-trivial variance in the model hidden states, while the matrix exponential method consistently explained the most variance by a sizable margin. This provided a little weak causal evidence that the models were using the linear multistep and Taylor expansion methods and stronger weak causal evidence that the models were using the matrix exponential method.

Criterion 4: Can we predictably intervene on the the model? Criterion 4.1: To intervene on the model, we used the reverse probes from Figure 7 to generate predicted hidden states from each intermediate. As shown in Figure 8, we then inserted these hidden states back into the model, to see if the model was still able to model the SHO. The matrix exponential method had the most successful interventions by an order of magnitude, and 18/25 of these intervened models performed better than guessing. This implies that the information the transformer uses to model the SHO is stored in the matrix exponential’s intermediate.

Figure 7. We found that the intermediates from all three methods could explain some variance in the undamped harmonic oscillator model hidden states, but that the matrix exponential method was the most consistent and successful by a wide margin.

Figure 7. We found that the intermediates from all three methods could explain some variance in the undamped harmonic oscillator model hidden states, but that the matrix exponential method was the most consistent and successful by a wide margin.

Figure 8. For each undamped harmonic oscillator model and method, we replaced the hidden state in Figure 7 with the reverse probe of the intermediate. We can see that this intervention was consistently the best performing for the matrix exponential method by an order of magnitude, and that 18/25 models performed better than our baseline of guessing.

Figure 8. For each undamped harmonic oscillator model and method, we replaced the hidden state in Figure 7 with the reverse probe of the intermediate. We can see that this intervention was consistently the best performing for the matrix exponential method by an order of magnitude, and that 18/25 models performed better than our baseline of guessing.

Criterion 4.2: We could also vary $\mathsf{\Delta }t\to \mathsf{\Delta }{t}^{\prime },{\omega }_0\to {\omega }_0^{\prime }$, regenerate intermediates and then hidden states, insert these modified hidden states into the model, and see if the model made predictions as if it “believed” that the input SHO data used $\mathsf{\Delta }{t}^{\prime },{\omega }^{\prime }$. As shown in Figure 9, we performed this intervention on $\mathsf{\Delta }t$, but our results were sufficiently robust to intervene on ${\omega }_0$ as well (Appendix C Figure A7). Even for the model with the best reverse probe quality for the linear multistep/Taylor expansion intermediates ($L=4,H=4$), the intervention with the matrix exponential method was the most successful. Combined with our previous intervention (4.1), we now had strong causal evidence for the matrix exponential method.

The transformer likely uses the matrix exponential to model the undamped harmonic oscillator. We now had correlational evidence that the model was using all three methods in our theory hub, with little causal evidence for the linear multistep and Taylor expansion methods, and strong causal evidence for the matrix exponential method. We suspected the model was only using the matrix exponential method in its computations, and that the evidence we had for the other two methods was a byproduct of the use of the matrix exponential. In Appendix C Figure A8, we give correlational evidence for this claim by generating synthetic hidden states from $e^{\mathit{A}\mathsf{\Delta }t}$ and showing that in this synthetic setup we retrieved values for criterions 1, 3 for linear multistep and Taylor expansion that were close to those we observed in Table 2.

Table 2. We summarize the evaluation of the methods and criteria for the undamped/underdamped models. For each criterion, we list a single quantity for readability: criterion 1 is the largest value in Figure 5, criterion 2 is the correlation in Figure 6, criterion 3 is the largest value in Figure 7, and criterion 4 is the ratio in the legend of Figure 8. We bold the best performing hypothesis for each criteria. The matrix exponential performed best across the criteria.

Criterion	Linear Multistep	Taylor Expansion	Matrix Exponential
1. Intermediate encoding	0.66/0.51	0.67/0.25	0.84/0.54
2. Performance, encoding correlation	0.73/0.44	0.74/0.39	0.89/0.44
3. Intermediate’s explanatory power	0.42/0.15	0.53/0.11	0.78/0.16
4. Intervention success	0.44/X	0.44/X	0.72/X

Table 2. We summarize the evaluation of the methods and criteria for the undamped/underdamped models. For each criterion, we list a single quantity for readability: criterion 1 is the largest value in Figure 5, criterion 2 is the correlation in Figure 6, criterion 3 is the largest value in Figure 7, and criterion 4 is the ratio in the legend of Figure 8. We bold the best performing hypothesis for each criteria. The matrix exponential performed best across the criteria.

Criterion	Linear Multistep	Taylor Expansion	Matrix Exponential
1. Intermediate encoding	0.66/0.51	0.67/0.25	0.84/0.54
2. Performance, encoding correlation	0.73/0.44	0.74/0.39	0.89/0.44
3. Intermediate’s explanatory power	0.42/0.15	0.53/0.11	0.78/0.16
4. Intervention success	0.44/X	0.44/X	0.72/X

Thus, we concluded that the transformer was likely using the matrix exponential method. This made sense, given the problem setting—both the linear multistep and Taylor expansion methods were only accurate for small $\mathsf{\Delta }t$, while our bound of $\mathsf{\Delta }t=U[0,2\pi /{\omega }_0]$ violated this assumption for some timeseries. On the other hand, the matrix exponential makes no such assumptions for the timestep, and it is, thus, a more general method that the transformer can employ to flexibly model a variety of situations. Still, it is remarkable that transformers use a known numerical method to model the undamped harmonic oscillator, and that we can provide evidence for its use, although our experiments do not rule out the possibility of other methods being used in conjunction with the matrix exponential.

4.4. Extension to the Damped Harmonic Oscillator ($\gamma \ne 0$)

We wanted to understand the generality of our finding by extending our problem space to the damped harmonic oscillator, where $\gamma \ne 0$. We have left the relevant details about our procedure to Appendix D, but, as shown in Table 2, we found that our intermediate analysis performed much more poorly on the underdamped case than on the undamped. We describe possible explanations in Appendix D, but because of this we temper our finding from the undamped harmonic oscillator with caution about its generality.

5. Discussion

After developing criteria for intermediates in the toy setting of linear regression, we found that transformers use known numerical methods for modeling the simple harmonic oscillator, specifically the matrix exponential method. We leave the door open for researchers to better understand the methods transformers use to model the damped harmonic oscillator and to use the study of intermediates to understand how transformers model other systems in physics.

Limitations. We analyzed relatively small transformers with only one attention head and no LayerNorm. While we demonstrated strong results for the undamped harmonic oscillator, our results for the underdamped harmonic oscillator were more mild. We only used noiseless data.