You can pass cmin = min and cmax = max inside add_surface() to get a consistent color scale across all plots:

add_surface(
    z = z_matrix,
    x = p1,
    y = p2,
    colorscale = "Viridis",
    cmin = 0,
    cmax = 4
)