Why substeps affect the speed of movement on a simple verlet physics simulation?

Question

I have made a simple Verlet solver which works really well, but changing the substeps for the simulation results in forces being stronger the more substeps there are.

Important pieces in pseudo code:

// Every frame
let subDeltaTime = deltaTime / substeps;// In seconds
for (let i = 0; i < substeps; i++) {
    for (let body of simulation.bodies) {
        // Gravity
        body.acceleration.y += 9.8 * body.mass;
    
        // Solve constraints and collisions
    
        // Update position using Verlet
        let displacement = body.position - body.oldPosition;
        body.oldPosition = body.position;
        body.position += displacement + body.acceleration * (subDeltaTime * subDeltaTime);
        body.acceleration *= 0; // Reset acceleration
    }
}

Relevant code and repository is here, if you want to take a look deeper.

I have tried:

• Dividing the forces being applied by the amount of substeps.
• Accumulating the acceleration and applying it only once.
• Fixing the frame rate and delta time to respectively 60FPS and 16ms.
• Looking at other implementations but couldn't really find anything else.

Minimal example demoing the issue:

const canvas = document.getElementById("canvas");
canvas.width = 400;
canvas.height = 150;
const ctx = canvas.getContext("2d");

let lastFrameTime = 0;
let frameCounter = 0;
let frameTimer = 0;
let fps = 0;

let substeps = 4;

let bodies = [];
class Body {
  constructor(x, y, radius) {
    this.x = x;
    this.lastX = x;
    this.y = y;
    this.lastY = y;
    this.accX = 0;
    this.accY = 0;
    this.radius = radius;
  }
}

document.getElementById("substeps-slider").oninput = function() {
  substeps = this.value;
}

function setup() {
  for (let i = 0; i < 50; i++) {
    bodies.push(new Body(Math.random() * canvas.width, Math.random() * canvas.height, 7.0));
  }
}

function render(now) {
  // Clear
  ctx.fillStyle = "black";
  ctx.fillRect(0, 0, canvas.width, canvas.height);

  // Delta time
  let elapsed = now - lastFrameTime;
  let dt = elapsed / 1000.0;
  lastFrameTime = now;

  // Bodies - update
  let delta = dt / substeps;
  for (let i = 0; i < substeps; i++) {
    for (const body of bodies) {
      // Apply constraints
      if (body.x < body.radius) body.x = body.radius;
      if (body.y < body.radius) body.y = body.radius;
      if (body.x > canvas.width - body.radius) body.x = canvas.width - body.radius;
      if (body.y > canvas.height - body.radius) body.y = canvas.height - body.radius;

      // Collision
      for (const other of bodies) {
        if (other === body) continue;

        const distance = Math.sqrt(Math.pow(body.x - other.x, 2) + Math.pow(body.y - other.y, 2));
        const max = other.radius + body.radius;
        if (distance < max) {
          let normalX = (body.x - other.x) / distance;
          let normalY = (body.y - other.y) / distance;

          let overlap = max - distance;

          body.x += normalX * overlap * 0.5;
          body.y += normalY * overlap * 0.5;
          other.x -= normalX * overlap * 0.5;
          other.y -= normalY * overlap * 0.5;
        }
      }

      // Apply gravity
      body.accY += 9.8 * delta * delta;

      // Apply verlet
      let displacementX = body.x - body.lastX;
      let displacementY = body.y - body.lastY;
      body.lastX = body.x;
      body.lastY = body.y;

      body.x += displacementX + body.accX * (delta * delta);
      body.y += displacementY + body.accY * (delta * delta);

      body.accX = 0;
      body.accY = 0;
    }
  }

  // Bodies - render
  for (const body of bodies) {
    ctx.fillStyle = "red";
    ctx.beginPath();
    ctx.arc(body.x, body.y, body.radius, 0, 2 * Math.PI, false);
    ctx.fill();
  }

  // FPS
  frameTimer += dt;
  frameCounter++;
  if (frameTimer >= 1.0) {
    frameTimer -= 1.0;
    fps = frameCounter;
    frameCounter = 0;
  }

  // HUD
  ctx.fillStyle = "white";
  ctx.fillText("Frame time: " + Math.round(elapsed).toString().padStart(2, "0") + "ms", 10, 10);
  ctx.fillText("FPS: " + fps, 10, 20);
  ctx.fillText("Substeps: " + substeps, 10, 30);
  ctx.fillText("Bodies: " + bodies.length, 10, 40);

  requestAnimationFrame(render);
}

setup();
lastFrameTime = performance.now();
requestAnimationFrame(render);

div {
  display: flex;
  flex-direction: column;
  place-content: center;
  place-items: center;
  font-family: "Roboto", "Arial", "sans-serif"
}

input {
  margin-top: 1pt;
}
p {
  margin: 1pt 0;
}

#root {
  height: 500px;
}

<div>
  <small>You may need to scroll down</small>
  <canvas id="canvas"></canvas>
  <input type="range" min="1" max="20" value="1" class="slider" id="substeps-slider" style="width: 400px">
  <p>Move the slider to change the substeps of the simulation</p>
  <small>Notice that the more substeps, the faster the simulation runs</small>
</div>