PhySim

A virtual science environment

A central challenge in using AI systems to solve problems in science is the need to interact with the world through experiments. Compared to fields like mathematics, it is likely not possible to make fundamental progress in science through thought alone. The cost of experiments, either economic or time, constrains the rate at which AI can build its understanding of physics, chemistry, and biology.

While improving experimental throughput in reality is a monumental task, we can ask what AI can learn about a virtual world where experiments can be done cheaply in-silico. PhySim was built to enable this. It defines specifications for

Virtual worlds in which complex emergent phenomena can arise from simple rules
An agent evaluation suite that provides tools for experimentation and tests agents’ ability to forecast world dynamics

While these worlds can be designed by hand, we also provide a package for breeding interesting worlds with an evolutionary algorithm. This combination enables us to measure model performance, use RL to improve this performance, and generate increasingly complex worlds that push models to their limits.

A PhySim virtual world with complex emergent behavior

PhySim Worlds

PhySim worlds are built from two-dimensional fields. These fields are functions with two-dimensional domains; for every point in two-dimensional space, they define a single numerical value.

The four BF fields at one moment: round blobs in u₀ with matching halos in x₀ and x₁, and winding trails in x₂ — An example world called “BF” is defined as a collection of four fields. Each panel shows one field over the same region, on that field’s own fixed color scale.

While fields are defined everywhere, when large field values concentrate in one spatial location we can think of the field as describing a localized object. In PhySim we call these objects “blobs”. To highlight this behavior, we refer to a world called “BF” throughout this post. The fields in BF all exhibit this behavior, with some blobs appearing as spots constrained to a single locus and others tracing out elongated paths through the domain.

Given a set of initial conditions, we evolve all the fields in time according to diffusion laws and couplings, or interactions, between fields. The laws of this evolution come from field equations which describe the rate of change of each field at each moment in space.

BF fields evolve in time according to their field equations and a choice of initial conditions.

View the full BF field equations

∂_tu₀ = ∇²u₀ + 2 · u₀ − u₀ ³ + 1.059 − x₀ − 1.5 · x₁ − x₂ · x₁

∂_tx₀ = 0.7018 · ∇²x₀ + (u₀ + 0.7035 − x₀) / 5.7

∂_tx₁ = 20 · ∇²x₁ + (u₀ + 0.7035 − x₁) / 0.7

∂_tx₂ = (0.05 · tanh(max(u₀ − 0.2494, 0) / 0.4) − x₂) / 200

Download the world JSON · Published source

It is convenient to think of worlds and their fields as continuous objects, but practically we simulate dynamics on a grid of discrete points. Simulation additionally depends on a choice of noise realization. That means that, along with initial conditions, choices of grid size/spacing and noise also affect the dynamics of a given world simulation. The simulator source in the repository shows exactly how the noise enters each step.

The whole periodic domain of u₀, a 16 by 16 unit zoom on one blob, and a 6 by 6 block of grid points with their stored values — BF is simulated on a 256 × 256 grid. Zooming in on one blob in u₀ reveals the discrete values the simulator stores.

A set of fields and field equations defines the world, and a set of initial conditions and simulator settings provides a single realization of the world’s dynamics. From these basic descriptions, worlds can exhibit complex emergent phenomena such as interaction between blobs and the rise or decay of blob populations.

Experimenting with virtual worlds

In science we do not have the privilege of knowing the underlying framework behind mechanisms we observe and must develop theories which explain empirical data we collect. When building the agent evaluation suite for PhySim, we likewise chose to keep most details of the virtual worlds hidden from the agent.

Running experiments

Agents are given the following experimental apparatus:

Device 0 with 13 sensors and device 1 with 19 sensors, both centered on a blob and drawn on all four BF fields — The apparatus placed in BF. Both devices are centered on a blob, with the source at their shared center. Every sensor reads every field, so the same markers appear on each panel. The agent never sees this picture.

Each device consists of a lattice of sensors, an input source, and controls for adjusting positioning.

Sensors: record the value of each of the fields at the sensor’s location, interpolating between nearby simulation grid points.

A sensor between four grid points in u₀: each point's value, its weight, and the weighted sum that becomes the sensor's reading — A sensor between grid points reads a weighted average of the field at the four nearest points; the closer a point, the larger its weight. The same weights apply to every field.

Sources: allow injecting pulses to each of the fields. The pulses are emitted from the center of the device and have a gaussian distribution.

A pulse into the trail field x₂ at device 0. While the source is on, it adds a Gaussian bump with a standard deviation of 2 units, centered on the device. Right: x₂ along the dashed line through the device’s center, before the pulse (gray) and as it builds, with the five sensors on that line marked.

Adjustments: the entire device can be moved in either direction within the world grid, as well as dilated.

Device 0 shifts one unit down at t = 0, 5 and 10, one unit right at t = 15, 20 and 25, and dilates by about 1.5× at t = 30. Right: the 13 sensors’ readings of x₂, which change as the device moves. A faint outline shows the device just before each move.

However, the agent only sees a nondescript experimental interface:

actions = [
    {"t": 1.0, "kind": "inject", "device": 0, "port": 2, "amp": 0.7, "dur": 4.0},
    {"t": 2.0, "kind": "adjust", "device": 0, "u": [0.2, -0.4, 0.1]},
]
queries = [
    {"sensor": "device0", "t": [0, 3, 12]},
    {"sensor": "device1", "t": [12, 20]},
    {"sensor": "global", "t": [20]},
]

Each experiment consists of actions (changes to the world the agent can specify, these can be either source injections or device adjustments) and queries (readings recorded from the sensors at specified times). The interface is intentionally opaque and agents need to learn the meaning of parameters like t (time), port (which field the source affects), and the adjustment coordinates u (which correspond to left-right, up-down, and dilation inputs, respectively). In addition to the two devices shown above, agents can query sensors from a global device, which records the spatial mean and variance of each field. Results from each experiment are returned in arrays. Below we see what those results look like for the actions and queries shown above.

# After laboratory_experiment(actions, queries):
# Example output in BF (four ports).
observations = np.load("/observations/experiment_001.npz")
# Axes: realization, times, ports, slots
observations["query0"].shape  # (1, 3, 4, 13)
observations["query1"].shape  # (1, 2, 4, 19)
observations["query2"].shape  # (1, 1, 4, 2)
# 13 readings at time 3 on port 2:
observations["query0"][0, 1, 2, :]

Crucially, the agent is not told anything about the structure of the world; that it exists on a 2d grid, how the devices and their sensors are set up, or even what the inputs and adjustments actually do. Likewise, sensor readouts are just opaque arrays of data with no explicit meaning attributed to them. This means that agents need to learn how the experimental apparatus works by trial and error. Think of this as being born with sense organs, but needing to learn through experience how to make use of their data streams.

This has interesting implications for interpreting agent rollouts. Agents are free to build their own world models, and the exact shape of those models does not affect how we score agents’ predictions (see below). In particular, agents may come up with diverse theoretical frameworks for describing each world, and we do not judge those frameworks on how close they are to the underlying world description we have chosen. Still, we must grade agents somehow, and choose the following general approach.

Making predictions

Agents must submit a prediction function, which is able to forecast the world response to an arbitrary experiment (that is, a set of actions and queries). Agents do not know ahead of time which actions and queries we will choose to evaluate these predictions, and must build general models of world dynamics.

def predict(actions, queries, n_samples=64, seed=0):
    # One array per query, in the same order.
    # Shape: (n_samples, requested times, n_ports, sensor slots)
    return {"samples": arrays}

We allow agents to predict a set of sample outputs rather than a single output, which effectively makes predictions distributional. In addition to running experiments and submitting their final prediction function, we also supply agents with a validation tool which confirms their prediction adheres to the interface we expect. It checks execution, output format, finite values, and repeatability on public example requests. Passing those checks does not guarantee valid predictions for every possible experiment.

Evaluating agents

The full agent evaluation loop looks like this:

Agents experiment with the world
Agents validate (optional) and submit their final prediction function
We measure how closely the submitted prediction function agrees with simulation for a fixed (and a priori unspecified) set of experiments

In this section, we explain how the third step works in more detail.

Even under finite time constraints simple combinatorics tell us that there are a vast number of possible experiments we could use to judge predictions. However, physical intuition tells us that the majority of these will be uninteresting. After all, worlds may exhibit localized excitations, both spatially (blobs surrounded by empty space) and temporally (worlds that undergo interesting dynamics and eventually settle into stable equilibria). For this reason, we require that any world we use for agent evaluation come equipped with a curated set of experiments which are designed to capture interesting behavior in that world’s dynamics.

For example, the blobs in BF’s field u₀ are subject to a feedback loop mediated by the trail field x₂. When blobs in the field u₀ move they leave trails of elevated values in the field x₂. These trails then feed back into u₀’s field equations through an interaction term and affect the trajectory of the u₀ blobs. One of the BF evaluation experiments, shown below, tests this directly. It injects a pulse into x₂ at the center of device 0 for five time units, and the sensors record the response for 50 time units. The pulse leaves a lasting spot in the trail, and about 20 time units after the pulse ends, the nearby blob in u₀ is pushed off that spot. To forecast this delayed response, a model must capture both what the pulse does to x₂ and how u₀ is affected by this later.

The high trail pulse experiment, seen through device 0’s 13 sensors. Dots mark the readings of u₀ and x₂ that are scored, at the times they are taken; the shaded band is the pulse.

Other experiments in the BF evaluation include pulse strengths in the trail field x₂, pulses in other fields, and changes in device positions. The full set of experiments is described below.

The 15 BF evaluation experiments
Experiment	Actions
No actions	None; the world runs on its own.
Weak u₀ pulse	A weak pulse into u₀ at device 0, from t = 0 to 5.
Strong u₀ pulse	The same pulse, six times stronger.
Low trail pulse	A weak pulse into the trail field x₂ at device 0, from t = 0 to 5.
Medium trail pulse	The same pulse, five times stronger.
High trail pulse	Five times stronger again; the experiment shown above.
u₀ pulse, then trail pulse	A strong u₀ pulse from t = 0 to 5, then a medium trail pulse from t = 10 to 15.
Delayed u₀ pulse	A strong u₀ pulse from t = 20 to 25.
x₁ pulse	A pulse into the fast inhibitor x₁ at device 0, from t = 0 to 5.
Shift device 0	Device 0 shifts one unit down at t = 0; nothing is injected.
Dilate device 1	Device 1 dilates by about 2.7× at t = 0; nothing is injected.
u₀ pulse, then shift	A strong u₀ pulse starts at device 0, which shifts one unit down at the same moment; the pulse stays at the old center.
Shift, then u₀ pulse	Device 0 shifts one unit down, then a strong u₀ pulse starts at its new center.
Two sources	Device 1 shifts one unit right, then both devices pulse u₀ from t = 0 to 5.
Shift during u₀ pulse	A strong u₀ pulse at device 0; at t = 2 the device shifts one unit down while the pulse continues where it started.

Once we have a set of experiments built for a world, we combine them into a single score using the following recipe.

A group selects readings to score jointly, such as one field at several sensors and times. Every BF experiment is scored on the same four groups, listed below. For a fixed experiment and group, let f_ik be reading k of forecast sample i and o_jk the same reading in simulated realization j. Divide both by predeclared positive scales s_k in their native units:

{\hat{f}}_{i k} = \frac{f_{i k}}{s_{k}}, {\hat{o}}_{j k} = \frac{o_{j k}}{s_{k}}

The scales s_k set the relative weight of the readings. A smaller scale makes the same absolute prediction error count more. Scales can account for expected differences in field magnitudes, or emphasize a field whose response an experiment is designed to detect. They are chosen before evaluating predictors and remain fixed for that experiment. In BF, an error of 0.05 in the trail field x₂ counts as much as an error of 1 in u₀, matching the trail’s much smaller values.

The four score groups of every BF experiment
Readings	Times	Scale
u₀ at device 0’s 13 sensors	5, 10, 20, 30, 40, 50	1
x₀ at device 0’s 13 sensors	5, 10, 20, 30, 50	0.5
x₂ at device 0’s 13 sensors	5, 10, 20, 30, 50	0.05
u₀ at device 1’s 19 sensors	10, 30, 50	1

x₁ and the global sensor are recorded in every experiment but not scored.

For these dimensionless vectors with d coordinates, define the root mean square distance:

D (a, b) = \sqrt{\frac{1}{d} \sum_{k = 1}^{d} {(a_{k} - b_{k})}^{2}}

With M predicted samples and N independent simulated realizations, the energy score for group g in experiment c is:

E_{c g} = \frac{1}{M N} \sum_{i = 1}^{M} \sum_{j = 1}^{N} D ({\hat{f}}_{i}, {\hat{o}}_{j}) - \frac{1}{2 M^{2}} \sum_{i = 1}^{M} \sum_{ℓ = 1}^{M} D ({\hat{f}}_{i}, {\hat{f}}_{ℓ})

\begin{aligned} E_{c g} & = \frac{1}{M N} \sum_{i = 1}^{M} \sum_{j = 1}^{N} D ({\hat{f}}_{i}, {\hat{o}}_{j}) \\ - \frac{1}{2 M^{2}} \sum_{i = 1}^{M} \sum_{ℓ = 1}^{M} D ({\hat{f}}_{i}, {\hat{f}}_{ℓ}) \end{aligned}

Average equally over the G_c groups in each experiment, then over the C experiments:

S_{c} = \frac{1}{G_{c}} \sum_{g = 1}^{G_{c}} E_{c g}, S = \frac{1}{C} \sum_{c = 1}^{C} S_{c}

For BF, M = 64, N = 2, G_c = 4 for every experiment, and C = 15. The quantity S defines the total energy, which can be interpreted as the prediction error. To shape this into a useful reward for training, we use the following formula:

R_{K} (S) = \frac{\max (0, \min (- \log_{10} S, K))}{K}

At S = 0, define R_K = 1. A missing predictor or a failed prediction contract receives zero reward.

The shaping of the reward requires some explanation. The idea is to reward successive orders of magnitude of reduction in normalized prediction error. Explicitly, the function −log₁₀(S) maps a prediction error of 10⁻ⁿ to n. So the reward is essentially picking out this exponent, clamping it to the window [0, K], and finally scaling so it lies between 0 and 1. K is a hyperparameter which sets the precision target for the reward. If we loosely think of 1 − S as the accuracy of the prediction, this reward shape is analogous to counting “9’s” of accuracy. Before normalization (dividing by K), an energy of 0.1 corresponds to an accuracy of 0.9, which gets an unnormalized reward of 1, while an energy of 0.01 corresponds to an accuracy of 0.99 and an unnormalized reward of 2.

For the results below, we use K ≈ 2.746305482, giving full reward at S ≤ 0.001793472. These rewards are recomputed from saved predictions; the agents received K = 2 in their prompts.

Precision calibration

We generated 64 forecast samples using the actual world equations, with independent realizations of the world’s noise. These were scored against two saved simulation realizations for every experiment in the BF, XV, and p4g2_044 suites. BF had the largest mean energy, 0.001793472. We set K to −log₁₀ of this energy, without integer rounding or a margin. This makes the full-reward threshold equal to the largest native-simulator suite energy. Even a correct stochastic predictor can have nonzero energy; this calibration is specific to these experiments and measurement scales. Native calibration measurements (including the earlier integer-K rule) · Current model rewards.

Preliminary Results

We performed rollouts with a range of models to study how agents behave in PhySim environments. All rollouts were conducted in the BF world described above. Native harnesses (Codex and Claude Code) were used for GPT and Claude models, respectively, while all other models used the Prime Agent harness. Overall performance is shown in the table below.

One BF investigation per model · 15 held-out experiments · 6–7 October 2026
Model / predictor source	Reward ↑	Energy ↓	Experiments	Inference cost
GPT-6 Astra Codex	0.796	0.0065	69	$23.09
GPT-6.1 Sol Codex	0.745	0.0090	142	$7.56
Claude Opus 5.5 Claude Code	0.707	0.0114	227	$90.99*
Claude Sonnet 5 Claude Code	0.324	0.1291	384	$11.42*
GPT-6 Luna Codex	0.291	0.1590	536	$1.14
GPT-5.6 Terra Codex	0.231	0.2327	100	$25.85*
Qwen3.5 397B A17B Prime Agent	0.164	0.3542	224	$1.57*
Qwen3.5 35B A3B Prime Agent	0.077	0.6128	26	$0.94*
Grok 4.6 Prime Agent	0.000	1.8576	104	$42.75*

* Estimated from token usage and pricing rates. Unmarked costs are reported directly by the API.

Experiment record and methodology

The table entries above average over different experimental predictions agents must make, each specific to the BF world. The heatmap below shows performance on each experiment for all models.

Energy scores for nine completed predictors across fifteen BF experiments, including Claude Opus 5.5, GPT-6 Astra, and GPT-6.1 Sol — Entries are joint energy per experiment, averaged over its four score groups; lower is better. Colors use a logarithmic scale.

What agents built

While all agents performed experiments and developed their own models for world dynamics, the complexity of those models varied across the rollouts. The top scoring agents all developed inductive biases which matched the hidden underlying physics, building spatial reaction-diffusion models. Other models took a black box approach, generating data through experiments and fitting neural networks or regression models. The simplest predictors relied on matching or interpolating observed responses, with little explicit theory or world model development. Full details of the submitted prediction functions can be seen here.

Discovering new worlds

While the task above is already non-trivial, we believe PhySim environments can be made far more challenging if we can generate more complex worlds. To discover equations that gave rise to more complex behavior we applied an evolutionary algorithm which selected for interesting behavior and bred worlds by combining their fields and adding couplings between them.

See details of the evolutionary world search

Initialize. Begin with a population of starting worlds and a seeded random generator.
Vary. Mutate numerical parameters or field structure. To combine two parents, first join their field systems, then share a feedback field, add cross-couplings, or introduce a slowly relaxing field driven by both parents.
Simulate and measure. Reject unsupported or unstable candidates. Simulate the remainder and calculate metrics in Python, such as persistence, movement, spatial organization, and interactions. Selection uses these measurements, not an LLM judge.
Select and repeat. A recipe converts the metrics into fitness and chooses survivors and parents for the next generation. Recipes can define their own metrics and selection rules; the packaged default retains the highest-scoring distinct worlds.
Preserve. Harvest selected worlds, recipe, seeds, parentage, metrics, and checkpoints in the registry.

Search implementation · Mutation and merge operations · Example recipe

We package the simulation tools along with the world breeding algorithm into a package called Blobkit. Since the evolutionary search is constrained by simulation throughput, we also include accelerated simulation code to run multiple worlds in parallel.

An example of what this process produced is the world we call “p4g2_044”. This world is described by twelve fields and produces dynamics far richer than BF. It was found in the second generation of the evolutionary search we ran. We evaluated GPT-6 Astra on this world; its aggregate reward was 0.573, corresponding to roughly four times its error on BF. Highlights of the world dynamics are illustrated below.

Four of p4g2_044’s twelve fields over 450 time units, with nothing injected. Blobs in u₀ travel, split and line up into a chain; the stripes in u₁ reorganize around them; a labyrinth in u₃ coexists with both; and x₇ keeps a slowly fading record of where the blobs have been. White marks each field’s resting value; u₃ rests high, so its scale runs the other way and dark marks its low channels.

View the full p4g2_044 field equations

∂_tu₀ = 1.083 · ∇²u₀ + 1.948 · u₀ − u₀ ³ + 1.01 − x₀ − 1.5 · x₁ − 0.01273 · x₇

∂_tu₁ = 0.65 · ∇²u₁ + 2 · u₁ − u₁ ³ + 1.082 − 2.05 · x₁ − x₂ − 0.1472 · x₃

∂_tu₂ = 0.4411 · ∇²u₂ + 2.852 · u₂ − u₂ ³ − 1.247 − 1.964 · x₄ − 0.8124 · x₅ − 0.3082 · x₆ − 0.01273 · x₇

∂_tu₃ = 0.5266 · ∇²u₃ + 2.91 · u₃ − u₃ ³ − 1.417 − 1.964 · x₄ − 0.8124 · x₅ − 0.3082 · x₆ − 0.01273 · x₇

∂_tx₀ = ∇²x₀ + (u₀ + 0.6798 − x₀) / 3.479

∂_tx₁ = 20 · ∇²x₁ + (u₀ + 0.5 · u₁ + 1.114 − x₁) / 0.7

∂_tx₂ = ∇²x₂ + (u₁ + 0.8676 − x₂) / 3

∂_tx₃ = 0.03183 · ∇²x₃ + (u₁ + 0.8676 − x₃) / 28.46

∂_tx₄ = 11.95 · ∇²x₄ + (u₂ + u₃ − 2.77 − x₄) / 78.89

∂_tx₅ = 5.09 · ∇²x₅ + (u₂ + u₃ − 2.77 − x₅) / 1.882

∂_tx₆ = 8.601 · ∇²x₆ + (u₂ + u₃ − 2.77 − x₆) / 6.085

∂_tx₇ = 0.4244 · ∇²x₇ + (0.1057 · tanh(max(u₀ + 0.02073, 0) / 0.3031) + 0.1057 · tanh(max(u₂ − 2.06, 0) / 0.3031) + 0.1057 · tanh(max(u₃ − 2.029, 0) / 0.3031) − x₇) / 52.28

Download the world JSON · Published source

What’s next?

The evaluations above serve as a proof of concept for PhySim. We think there are many ways to extend this work. We’d like to run more evaluations on the most capable models available now to see how well they do on worlds like BF. We would also like to put more effort and research into how to build complex worlds. In particular, it would be interesting to see what sorts of worlds can be produced with larger budgets for evolutionary search, and what sorts of selection metrics produce results that most closely align with what humans perceive as “interesting”.

The most interesting question may be what role environments for virtual science might play in the overall development of AI for science. While predicting the future of this technology is best left to the experts, we can try to reason from first principles about the inputs and outputs to scientific research and make the following conclusions:

Scientific progress will be bottlenecked by experimental throughput
To optimize this bottleneck, systems will need full flexibility to design new experimental methodologies
Our theoretical models for the natural world will be judged not only by their epistemic value, but also their ability to increase experimental throughput (simulation).

Get Involved

Interested in collaborating on this project?

World contribution guide if you want to build or grow your own PhySim worlds
PhySim GitHub if you want to explore the libraries supporting PhySim
Prime Intellect environments if you want to contribute to the agent evaluation and training environments
@the1pun1isher if you just think the project is cool and want to get in touch

Credit

This work was done as part of Prime Intellect’s RL residency. I would like to thank Prime Intellect for supporting the compute and token costs of this work, as well as @omouamoua for organizing this residency cohort and for all his input and feedback on this project. PhySim was partially inspired by @patience_cave’s work on MazeBench (be sure to check it out, it’s awesome). I’m also thankful to Yajit for many discussions on AI, agents, and science which have shaped parts of this work.