Eigenia Monte Carlo Engine: Importance Sampling & Graph Random Walks

J. McKenney

This is the Monte Carlo engine paper of the WG-08-MO series, which covers Monte Carlo and annealing methods. Its three siblings are WG-08-MO-02 on stratified Latin hypercube sampling over attack graphs, WG-08-MO-03 on quantum annealing and the Ising formulation of grid islanding, and WG-08-MO-04 on thermodynamic annealing under Levy flights. This paper is the engine those three sample with: it defines the unified graph, the importance-weighted traversal, the Boltzmann dynamics and the physical solvers, and it is the document to read first in the series.

Licence: CC BY 4.0. 17 September 2026.

Executive Abstract#

An asset scanner that finds an unpatched vulnerability on an engineering workstation and a network path to a controller reports that the controller is reachable. That answer is binary, and in an operational plant it is the wrong shape of answer. It carries no traversal time, no adversary skill, no probability of detection, and no physics: a controller regulating an office ventilation baffle and one regulating a high-pressure boil-off gas compressor score the same.

This paper specifies the engine that replaces the binary answer with a distribution. It binds a DEXPI 2.0 process graph to a CycloneDX multi-bill-of-materials layer so one graph carries both the plant and the software inside it, extracts the subgraph that matters by importance-weighted breadth-first search, traverses it under a Boltzmann rule in which an energy barrier is a real physical quantity rather than a difficulty score, and couples three physical solvers where a traversal reaches equipment: convective cooling collapse, cryogenic acoustic slugging, and electrochemical thermal runaway in a battery installation.

The output is an actuarial distribution with a fat tail rather than a risk matrix cell, streamed over server-sent events and reproducible bit for bit from a seeded generator, so a run can be repeated and argued with. The engine's own conclusion is the barbell: capital spent on deterministic analog protection and on high-velocity monitoring does more than capital spent on the software inspection tier between them, because that tier is what an adversary who already has root execution is standing on.

Abstract#

This specification sets out stochastic graph physics, importance-weighted breadth-first search, and spectral path traversal over unified DEXPI 2.0 and CycloneDX multi-BOM graphs. Traditional industrial cybersecurity evaluates vulnerability management on static, unweighted graphs, and the binary reachability it produces fails in operational technology for three reasons: it cannot express traversal dynamics, it does not encode conservation of mass, momentum and energy, and it assumes losses cluster around a Gaussian mean when common-cause software vulnerabilities and shared fieldbuses put them in a power-law tail. This specification defines the engine that answers those failures: a unified graph schema binding DEXPI 2.0 to CycloneDX through a cross-domain property extension, a substation and battery reference architecture, importance-weighted breadth-first extraction with anchor and fallback mechanics, spectral decomposition and the temperature-modulated Boltzmann traversal rule, three multi-physics consequence solvers, a deterministic telemetry pipeline, and the actuarial metrics and barbell principle that follow, traced against the standards and regulations it implements.

1. Executive Problem Statement: The Failure of Binary IT Reachability#

Traditional industrial cybersecurity tools evaluate vulnerability management on static, unweighted graphs. An asset scanner discovers an unpatched Common Vulnerabilities and Exposures (CVE) identifier on an engineering workstation, queries a network access control list, and concludes that a path to a Programmable Logic Controller (PLC) exists. This binary representation (reachable versus unreachable) fails in mission-critical operational technology (OT) environments for three concrete reasons:

1. Absence of Traversal Dynamics: Binary reachability cannot calculate the traversal latency, the required adversary skill, or the probability that an intrusion is detected across an operational conduit before reaching physical equipment.
2. Disconnection from Physical Laws: A flat network topology does not encode conservation of mass, momentum, and energy. A compromised controller is treated identically whether it regulates a benign office ventilation baffle or a high-pressure $250\,\text{bar}$ boil-off gas compressor.
3. Gaussian Risk Distortions: Traditional risk matrices multiply subjective probability by static impact (1 to 5), assuming losses cluster around a Gaussian mean. In cyber-physical systems, common-cause software vulnerabilities and shared fieldbuses create extreme coupling. Losses follow an Extremistan power law where catastrophic single-event consequences dominate the cumulative tail risk.

The Eigenia Cyber Digital Twin (CDT) Monte Carlo Engine resolves these deficiencies. The Cyber Digital Twin models the entire world across seven architectural layers: from silicon roots and physical processes up through cognitive operator dynamics and actuarial risk. The engine transforms the physical and digital infrastructure of a facility into a unified, attributed directed graph. By executing real-time, Boltzmann-weighted random walks over live graph topologies enriched by spectral graph theory and real-time physical dissipation solvers, the engine calculates the empirical probability density function of cyber-physical catastrophe.

2. The Unified Graph Schema: Binding DEXPI 2.0 to CycloneDX 1.6+#

2.1 Formal Graph Definition ($G_{\text{CPDT}}$)#

The Cyber-Physical Digital Twin is modeled as a unified attributed directed multigraph:

The vertex set $V$ represents the union of physical engineering assets ($V_{\text{phys}}$) and cyber components ($V_{\text{cyber}}$):

Where:

• $V_{\text{phys}} = V_{\text{equipment}} \cup V_{\text{piping}} \cup V_{\text{valve}} \cup V_{\text{instrument}} \cup V_{\text{terminal}}$
• $V_{\text{cyber}} = V_{\text{hardware}} \cup V_{\text{software}} \cup V_{\text{crypto}} \cup V_{\text{operations}} \cup V_{\text{service}}$

The edge set $E$ represents physical fluid/power connections, digital communications conduits, and cross-domain bindings:

2.2 Semantic Ingestion of DEXPI 2.0 XML#

DEXPI 2.0 unifies the P&ID Specification 1.4 (Plant Model) and the Process Specification 1.0 (Process Model) into standardized DEXPI XML, replacing the legacy Proteus schema. The ingestion parser extracts:

1. Piping and Instrumentation Topology: Line numbers, nominal pipe diameters ($DN$), piping materials, valve discharge coefficients ($C_v$), and fail-safe positions (Fail-Closed, Fail-Open, Fail-Locked).
2. Process State Parameters: Fluid composition, operating temperatures ($T$), operating pressures ($P$), volumetric flow rates ($Q$), and enthalpy balances.
3. Instrument Functions: Sensing elements (temperature transmitters, differential pressure cells), actuators (pneumatic diaphragms, motorized gearboxes), and hardwired trip switches.

2.3 CycloneDX 1.6+ Multi-BOM Layering#

The digital properties of each asset are captured using a 5-BOM architecture within CycloneDX 1.6:

• HBOM (Hardware BOM): Silicon processors, microcontrollers, field-programmable gate arrays (FPGAs), variable frequency drive (VFD) power stages, printed circuit board (PCB) revisions, and hardware roots-of-trust (Caliptra RoT).
• SBOM (Software BOM): Real-time operating system (RTOS) kernels, embedded firmware binaries, open-source libraries (OpenSIL, OpenBMC), and programmable logic controller (PLC) application programs.
• CBOM (Cryptographic BOM): Device Identity Composition Engine (DICE) certificates, asymmetric key hierarchies, symmetric session keys, and post-quantum algorithms (ML-DSA, ML-KEM).
• OBOM (Operational BOM): Maximum allowable operating pressure (MAOP), critical trip setpoints, permitted actuator slew rates ($|dx/dt| \le \text{limit}$), and safety instrumented function (SIF) response times.
• VEX (Vulnerability Exploitability eXchange): Machine-readable vulnerability status records correlating open CVEs with runtime exploitability, CISA Known Exploited Vulnerabilities (KEV) status, and EPSS scores.

2.4 Cross-Domain Property Extension (dexpi:)#

To bind mechanical equipment to executing firmware, CycloneDX component records incorporate the standardized dexpi: property namespace:

{
  "type": "device",
  "bom-ref": "cdu-pump-vfd-01",
  "name": "Coolant Distribution Unit Primary Pump VFD",
  "version": "Rev-4.2",
  "properties": [
    {
      "name": "dexpi:plant:equipmentTag",
      "value": "PMP-CDU-101A",
      "description": "Equipment Tag in DEXPI 2.0 P&ID XML"
    },
    {
      "name": "dexpi:process:designFlowRateLpm",
      "value": "122",
      "description": "Nominal Volumetric Flow of PG25 Coolant"
    },
    {
      "name": "dexpi:safety:failState",
      "value": "Fail-Locked",
      "description": "Actuator mechanical state on control loop loss"
    },
    {
      "name": "dexpi:iec62443:zone",
      "value": "Zone-1-Cell-Level",
      "description": "Purdue Model Zone Assignment"
    },
    {
      "name": "dexpi:iec62443:slt",
      "value": "SL-3",
      "description": "Target Security Level per IEC 62443-3-3"
    }
  ]
}

3. Substation & BESS Electrical/OT Reference Architecture#

To demonstrate the concrete application of the unified graph schema, the engine maps the Endeavour Energy substation and battery energy storage system (BESS) network specification (Substation_spec_design.md) into a topological graph.

3.1 Substation & BESS Node Inventory (Table 1)#

3.2 Graph Edge Traversal Matrix (Table 2)#

Table 2: Substation and BESS attack conduit topology.

4. Subgraph Extraction: The Importance-Weighted BFS Algorithm#

Industrial facilities often contain more than $250{,}000$ logical and physical nodes. Running Monte Carlo simulations across an entire multi-million element matrix violates operational response constraints. To extract the critical operational subgraph in real time without losing systemic hazard fidelity, the engine uses an Importance-Weighted Breadth-First Search (IW-BFS).

4.1 Scoring Formulation#

The extraction frontier scores candidate nodes by topological centrality, active exploitability, and spectral connectivity:

Where:

• $\text{degree}(u)$ is the total degree of node $u$, serving as a computationally efficient proxy for betweenness centrality.
• $\text{EPSS}(u) \in [0, 1]$ is the Exploit Prediction Scoring System probability attributed to known CVEs residing on that asset.
• $\text{SpectralBoost}(u) \in [0, 1.8]$ is the normalized eigenvector centrality extracted from PostgreSQL (seldon.spectral_analysis), prioritizing topological bridgeheads.

// Production IW-BFS Scoring Implementation
export function computeImportanceScore(
  degree: number, 
  epssScore: number = 0, 
  spectralBoost: number = 0
): number {
  return degree * (0.3 + Math.max(0, epssScore)) * (1 + Math.max(0, spectralBoost));
}

4.2 Anchor Nodes and SLA Fallback Mechanics#

1. Top-20 Anchor Guarantee: The top-20 nodes ranked by global degree centrality are automatically injected into the simulation kernel, ensuring global structural choke points are preserved.
2. Three-Second SLA Fallback: If graph scoring and spectral join execution exceed a 3.0-second computational window, the engine gracefully drops spectral joins and executes a uniform topological BFS. This guarantees responsiveness in operational security consoles.

5. Spectral Graph Physics & Boltzmann Traversal Dynamics#

5.1 Spectral Graph Theory & Laplacian Decomposition#

Let $A$ be the adjacency matrix of graph $G$, and $D$ be the diagonal degree matrix where $D_{ii} = \sum_j A_{ij}$. The unnormalized graph Laplacian is defined as:

The normalized symmetric graph Laplacian is formulated as:

Eigenvalue decomposition of $\mathcal{L}$ yields:

• Algebraic Connectivity ($\lambda_2$): The Fiedler value ($\lambda_2$) measures the structural robustness of the network. A small $\lambda_2$ indicates sparse conduits connecting IT to OT zones, revealing natural choke points where defensive monitoring must be concentrated.
• Eigenvector Centrality: The principal eigenvector $v_n$ of the adjacency matrix provides the relative centrality score for every node: $$A v = \lambda_{\max} v \implies v_i = \frac{1}{\lambda_{\max}} \sum_{j \in \mathcal{N}(i)} v_j$$ Nodes in the top 50 ranks of eigenvector centrality in PostgreSQL (seldon.spectral_analysis) receive up to a $1.8\times$ weight boost, identifying critical lateral pivot points.

5.2 Physical Energy Barriers & Boltzmann Distribution#

Adversary traversal is modeled as a discrete-time Markov random walk governed by a Boltzmann probability distribution. Every network edge $e = (u, v)$ presents an energy potential barrier $\Delta E(u, v)$ determined by defensive friction:

Where:

• $\Delta E_{\text{base}} = 1.0\,\text{eV}$ baseline protocol hop friction.
• $\Delta E_{\text{firewall}} = 3.5\,\text{eV}$ for stateful deep packet inspection (DPI) matching IEC 62443 conduits.
• $\Delta E_{\text{crypto}} = 4.0\,\text{eV}$ for mutual TLS (mTLS) with hardware-enforced DICE certificates.
• $\Delta E_{\text{airgap}} = 12.0\,\text{eV}$ for physically isolated or serial-switched conduits. Unidirectional hardware data diodes set $\Delta E \to \infty$ in the reverse direction.

The probability of an adversary $a$ transitioning from node $u$ to adjacent node $v$ is:

Where:

• $k_B = 1.0$ is the normalized computational Boltzmann constant.
• $\mathcal{T}_{\text{eff}}(a)$ is the adversary's effective excitation temperature: $$\mathcal{T}_{\text{eff}}(a) = \mathcal{T}_0 \cdot \left( \frac{\text{ATQ}_a}{100} \right)^{1.85} \cdot \prod_s \mu_s(a)$$
• $\text{ATQ}_a \in [0, 100]$ is the Actor Threat Quotient.
• $\mu_s(a) \in [1.0, 2.5]$ is the TACAM sector-affinity multiplier.

5.3 TACAM Product Exposure ($B_2$ Edge Boost Map)#

During random walks, the engine correlates the adversary's TACAM product exploit portfolio against the destination node's Common Platform Enumeration (CPE). If the adversary possesses confirmed exploit capability against the target equipment (e.g. Volt Typhoon against Siemens S7, or Sandworm against Schneider Electric Modicon), the $B_2$ edge boost reduces the effective energy barrier:

Where $B_2 \in [1.5, 3.2]$, substantially increasing the probability of lateral penetration across that conduit.

5.4 Temperature Modulated Sweeps (Black Swan Mode)#

By modulating simulation temperature $T$, engineers evaluate two distinct failure regimes:

• Low Temperature ($T \to 0$): The adversary follows the path of least resistance. Models opportunistic, automated malware campaigns.
• Elevated Temperature ($T \ge 3.0$ / Black Swan Mode): Stochastic exploration forces paths across high-barrier conduits ($\Delta E > 8.0\,\text{eV}$), discovering non-intuitive lateral attack combinations that standard vulnerability scanners overlook.

6. Multi-Physics Physical Plant Coupling Solvers#

A cyber penetration in critical infrastructure is consequential only when it alters physical state variables. The engine couples graph step traversal to three first-principles physical solvers:

6.1 Convective Cooling Thermal Collapse Solver#

In high-density liquid cooling ($120\,\text{kW/rack}$ direct-to-chip AI clusters), shutting down secondary coolant distribution unit (CDU) pumps starves liquid flow. Silicon junction temperature $T_j(t)$ is governed by the energy balance:

Where:

• $P_{\text{die}} = 1{,}200\,\text{W}$ per accelerator package, the configurable maximum NVIDIA publishes for a GB200-class Blackwell GPU.
• $\dot{Q}_{\text{vol}}$ is the volumetric flow rate ($122\,\text{L/min}$ nominal per rack, collapsing to $0.0\,\text{L/min}$).
• $h_{\text{conv}}$ is the convective heat transfer coefficient, which drops exponentially under flow stagnation.
• $C_{\text{thermal}} = 800\,\text{J/K}$ is the lumped thermal capacitance of the stagnant cold plate assembly, dominated by the coolant retained in the channels once flow stops.

Under complete flow stagnation, silicon junction temperature climbs at $1.46^\circ\text{C/s}$ from a nominal $72.5^\circ\text{C}$, reaching the emergency hardware shutdown trip point ($94.0^\circ\text{C}$) at $14.8\,\text{seconds}$. Because SCADA polling and operator verification take $30\text{ to }90\,\text{seconds}$, software-only alerts cannot hold the compute load before the trip fires.

6.2 Cryogenic Joukowsky Acoustic Slugging Solver#

In LNG regasification terminals, spoofing suction pre-heater telemetry (+8.5°C false feedback) admits sub-cooled liquid methane droplets into reciprocating compressor cylinders running at $500\,\text{RPM}$. The acoustic pressure shock follows the Joukowsky equation:

With liquid methane density $\rho \approx 422\,\text{kg/m}^3$, acoustic velocity $c \approx 1{,}400\,\text{m/s}$, and piston impact velocity $\Delta v \approx 6.8\,\text{m/s}$:

Superimposed on the nominal $250\,\text{bar}$ discharge pressure, the peak internal stress exceeds $520\,\text{bar}$, blowing the cylinder head off the crankcase and causing catastrophic loss of containment.

6.3 Electrochemical BESS Thermal Runaway Solver#

In battery energy storage systems (BESS), overriding the battery management system (BMS) charging ceiling beyond $4.4\,\text{V/cell}$ triggers self-sustaining exothermic decomposition modeled by Arrhenius kinetics:

Once cell temperature exceeds the critical solid electrolyte interphase (SEI) decomposition point ($135^\circ\text{C}$), reaction heating outpaces convective cooling, driving cell temperature past $800^\circ\text{C}$ with explosive venting of hydrogen fluoride (HF) and flammable hydrocarbons.

7. Real-Time Telemetry Streaming & Mulberry32 Determinism#

7.1 Server-Sent Events Pipeline#

The simulation pipeline is exposed over HTTP Server-Sent Events at /api/mc-real/simulate/stream. Implemented as a TypeScript generator (sampleWalkSteps()), it yields fine-grained state telemetry for every discrete hop:

// Production SSE Walk Step Generator
export function* sampleWalkSteps(
  graph: UnifiedGraph,
  seedId: string,
  targetId: string,
  maxHops: number,
  temperature: number,
  rngSeed: number
): Generator<WalkStep, WalkSummary, void> {
  const prng = mulberry32(rngSeed);
  let currentId = seedId;
  let cumulativeCost = 0;
  let cumulativeProb = 1.0;
  const visited = new Set<string>([seedId]);

  for (let hop = 0; hop < maxHops; hop++) {
    const candidates = graph.getOutEdges(currentId);
    const selection = boltzmannSelect(candidates, temperature, visited, prng);
    if (!selection) break;

    const { edge, probability } = selection;
    cumulativeCost += edge.weight;
    cumulativeProb *= probability;
    visited.add(edge.targetId);

    const step: WalkStep = {
      hop,
      sourceId: currentId,
      targetId: edge.targetId,
      edgeType: edge.type,
      probability,
      cumulativeProb,
      cumulativeCost,
      zoneCrossed: edge.sourceZone !== edge.targetZone ? `${edge.sourceZone}->${edge.targetZone}` : null,
      timestamp: Date.now()
    };

    yield step;
    currentId = edge.targetId;
    if (currentId === targetId) break;
  }

  return {
    success: currentId === targetId,
    totalHops: visited.size - 1,
    finalProbability: cumulativeProb,
    totalCost: cumulativeCost
  };
}

7.2 Bit-for-Bit Determinism via Mulberry32 PRNG#

To support forensic audits, conformity assessment under Article 32 and Annex VIII of the EU CRA, and claims under policies carrying the state-backed cyber-attack exclusion that Lloyd's Market Bulletin Y5381 requires, every random walk is bit-for-bit deterministic. Passing a 32-bit unsigned integer rngSeed produces identical traversal trajectories across distributed test runners:

export function mulberry32(seed: number) {
  return function(): number {
    let t = (seed += 0x6D2B79F5);
    t = Math.imul(t ^ (t >>> 15), t | 1);
    t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
    return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
  };
}

8. Actuarial Risk Metrics & Fat-Tail Mathematics#

Traditional IT risk models calculate risk using Gaussian mean and variance. In operational technology, software vulnerabilities exhibit common-cause coupling: a single unauthenticated Modbus command shuts down all redundant cooling pumps or battery cooling loops simultaneously. Losses follow a Pareto fat-tail distribution ($P(L > x) \sim x^{-\alpha}$, where $1 < \alpha < 2$).

The engine computes five rigorous actuarial metrics across $10{,}000+$ simulation runs:

Table 3: Fat-Tail Actuarial Metrics#

8.1 The Barbell Defense Principle#

Because software security layers are susceptible to zero-day bypasses, supply chain tampering, and credential theft, relying solely on software firewalls increases tail risk. The Barbell Defense Principle mandates allocating security capital to two extremes:

1. Low-Risk Deterministic Analog Layer (80 to 90 percent of budget): Hardwired bimetallic switches, spring-loaded mechanical relief valves, pneumatic fail-safe actuators, and physical rupture discs completely decoupled from digital networks.
2. High-Velocity Digital Monitoring Layer (10 to 20 percent of budget): Real-time network telemetry, out-of-band optical sensors, and cryptographic hardware roots-of-trust (Caliptra).

Eliminating reliance on complex middle-tier software inspection firewalls ensures that even if an adversary gains root execution on supervisory servers, physical destruction remains impossible.

9. Regulatory and Standards Traceability Matrix#

The engine is not a research prototype and the standards below are not a reading list. Each one fixes something the engine has to get right to produce an answer an underwriter or a safety engineer can use. DEXPI fixes what a process graph is, CycloneDX fixes what the software inside it is, and the cross-domain property extension of section 2.4 is the only thing this specification adds to either. IEC 62443-3-3 fixes what a security level means, and IEC 61511-1 fixes what a safety instrumented function is, which is what the solvers of section 6 are computing against when they decide whether a protective action arrives in time. The Seveso III Directive and its Dutch implementation fix which installations are in scope for major-accident control, and the Cyber Resilience Act fixes what a manufacturer owes about the components the bill of materials lists. Lloyd's Market Bulletin Y5381 sits at the end because it decides what survives into a policy: it requires a state-backed cyber-attack exclusion, it binds managing agents rather than insureds, and no engineering evidence removes it, so what the distributions in section 8 move is the retained loss behind the exclusion rather than the exclusion itself.

• DEXPI 2.0 Specification (2025): Data Exchange in the Process Industry: Process and Plant Model Specification, DEXPI e.V.
• OWASP CycloneDX 1.6 / 1.7: Universal Software, Hardware, Operations, and Cryptography Bill of Materials Standard, ECMA-424.
• IEC 62443-3-3:2018: Industrial communication networks: Security for industrial automation and control systems: System security requirements and security levels.
• IEC 61511-1:2016: Functional safety: Safety instrumented systems for the process industry sector.
• Directive 2012/18/EU (Seveso III): Control of major-accident hazards involving dangerous substances, European Parliament.
• BRZO 2015: Besluit risico's zware ongevallen, Staatsblad van het Koninkrijk der Nederlanden.
• EU CRA (Regulation 2024/2847): European Cyber Resilience Act: Essential cybersecurity requirements for products with digital elements, Annex I & VII.
• Lloyd's (2022): Market Bulletin Y5381: Cyber-attack exclusions, Corporation of Lloyd's, 16 August 2022. The model wordings drafted to meet it are the Lloyd's Market Association's clauses LMA5564 to LMA5567, November 2021.

10. References#

The engine implements the standards and regulations listed below, each mapped to the clause it satisfies in the traceability matrix of section 9.

1. DEXPI e.V. DEXPI 2.0 Specification: Data Exchange in the Process Industry, Process and Plant Model Specification. 10 October 2025.
2. ECMA International. ECMA-424, CycloneDX Bill of Materials Specification. First edition, June 2024.
3. IEC. IEC 62443-3-3, Industrial communication networks, Security for industrial automation and control systems, System security requirements and security levels. Edition 1.0, 2013, with corrigendum 2014.
4. IEC. IEC 61511-1, Functional safety, Safety instrumented systems for the process industry sector. Second edition, 2016.
5. European Parliament and Council. Directive 2012/18/EU on the control of major-accident hazards involving dangerous substances. 4 July 2012.
6. Staatsblad van het Koninkrijk der Nederlanden. Besluit risico's zware ongevallen 2015.
7. European Parliament and Council. Regulation (EU) 2024/2847 on horizontal cybersecurity requirements for products with digital elements. 23 October 2024.
8. Corporation of Lloyd's. Market Bulletin Y5381: Cyber-attack exclusions. 16 August 2022. The model wordings drafted to meet it are the Lloyd's Market Association's clauses LMA5564 to LMA5567, 25 November 2021.

The per-accelerator power figure used in section 6.1 is NVIDIA Corporation's own published figure, given in its Datasheet for NVIDIA Blackwell Architecture, product datasheet.