Here's a very simple model of a Faux Falcon 9 launched vertically, with no turn towards horizontal. That doesn't matter so much at the altitude at max-Q but the final altitude at MECO is higher than in the videos because it hasn't turned towards horizontal. There are several simplifications, but it should reproduce most things in a qualitative way.
The final velocity is a little high but that may be related to the model not throttling back near max Q, or to other approximations.
I chose a scale height model for density $\rho(h) = \rho_0 \exp(-h/h_{scale})$ and a trans-sonic drag coefficient $C_D$ of 0.6 (from here) which matters mostly near mach 1 when max-Q is happening. I assumed the first stage fuel is 70% of the total launch mass of 550,000 kg.
Answer: Max-Q happens around mach 1 because the Earth's atmosphere and gravity and structural materials are what they are. Rockets are designed to make due with our atmosphere and gravity to get the most mass to orbit or beyond, with the caveat that they don't fall apart under crushing forces at max-Q.
If we lived on a planet with a lower surface pressure, it would happen earlier. If we lived on planet with different mass or diameter, that would affect both gravity on the rocket and the scale height, and max-Q would also happen earlier or later.
Luckily we don't live here!

def deriv(X, t):
h, v = X
acc_g = -GMe / (h + Re)**2
m = m0 - mdot * t
acc_t = vex * mdot / m
rho = rho0 * np.exp(-h/h_scale)
acc_d = -0.5 * rho * v**2 * CD * A / m
return [v, acc_g + acc_t + acc_d]
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint
Re = 6378137. # meters
GMe = 3.986E+14 # m^3/s^2
rho0 = 1.3 # kg/m3
h_scale = 8500. # meters
# faux falcon-9 FT
vex = 3600. # m/s
tburn = 160. # sec
m0 = 550000. # kg
mdot = m0 * 0.70 / tburn # kg/s
CD = 0.6
A = np.pi * (0.5*3.66)**2 # m^2
times = np.arange(0, tburn+1, 1) # sec
X0 = np.zeros(2) # initial state vector
answer, info = ODEint(deriv, X0, times, full_output=True)
h, v = answer.T
hkm = 0.001 * h
vkph = 3.6 * v
mach = v / 330. # roughly
rho = rho0 * np.exp(-h/h_scale)
Q = 0.5 * rho * v**2
if True:
plt.figure()
plt.subplot(2, 2, 1)
things = (hkm, vkph, mach, rho, Q)
names = ('height (km)', 'velocity (km/h)', 'mach', 'density (kg/m^3)', 'Q')
for i, (thing, name) in enumerate(zip(things, names)):
plt.subplot(5, 1, i+1)
plt.plot(times, thing)
if i == 2:
plt.ylim(0, 3)
plt.plot(times, np.ones_like(times), '-k')
llim, ulim = plt.ylim()
plt.text(5, 0.7*ulim, name)
plt.xlabel('time (sec)', fontsize=16)
plt.show()