### Theoretical analysis

We start by considering a spinor BEC, which is well described by the following Hamiltonian under single-mode approximation

$$\widehat{H}(q)={c}_{2}\frac{{\widehat{L}}^{2}}{2N}+{\displaystyle \sum _{{m}_{F}=-1}^{1}}({\mathit{qm}}_{F}^{2}-{\mathit{pm}}_{F}){\widehat{a}}_{{m}_{F}}^{\u2020}{\widehat{a}}_{{m}_{F}}$$(1)where

${\widehat{a}}_{{m}_{F}}$(

${\widehat{a}}_{{m}_{F}}^{\u2020}$) is the annihilation (creation) operator for the spin *m _{F}* component corresponding to the hyperfine level ∣

*F*= 1,

*m*〉,

_{F} is the condensate’s total spin operator along μ (μ = *x*, *y*, *z*) with *f*_{μ} being the corresponding spin-1 angular momentum matrix, *c*_{2} is the spin-dependent interaction (*c*_{2} > 0 for the AFM sodium atoms), *N* is the total atom number, and *q* (*p*) is the quadratic (linear) Zeeman energy.

In the absence of the linear Zeeman energy (*p* = 0), there are two phases for the ground state: a polar phase with atoms all occupying the *m _{F}* = 0 level for

*q*> 0 and the AFM phase with atoms equally occupying the

*m*= ± 1 levels for

_{F}*q*< 0 (

*34*). If we take mean value 〈ρ

_{0}〉 with

as an order parameter, we can clearly see that 〈ρ_{0}〉 abruptly drops from one to zero at *q* = 0, showing the first-order quantum phase transition there (see Fig. 1A). At the transition point *q _{c}* = 0 Hz, these two phases coexist. Near this point, we can observe the existence of the polar phase for

*q*< 0 and AFM phase for

*q*> 0 as metastable states, which is the characteristic of the first-order phase transitions. In real experiments,

*p*is nonzero. However, because the Hamiltonian commutes with the total magnetization

, i.e.,

$[\widehat{H}(q),{\widehat{L}}_{z}]=0$, the quench dynamics is restricted in the subspace with zero magnetization if we prepare the initial state in the polar phase and the linear Zeeman term therefore becomes irrelevant.

To simulate the scaling in the dynamics across the first-order quantum phase transition, we start with the ground state of a spinor condensate in the polar phase for positive *q _{i}* and then linearly vary the quadratic Zeeman energy

*q*by

*q*(

*t*) =

*q*−

_{i}*vt*with

*q*>

_{i}*q*,

_{c}*q*<

_{f}*q*, and

_{c}*v*= (

*q*−

_{i}*q*)/τ

_{f}*characterizing the quench rate with τ*

_{q}*being the total time as*

_{q}*q*changes from

*q*to

_{i}*q*. To numerically simulate the dynamics, we solve the Schrödinger equation

_{f} by directly diagonalizing the many-body Hamiltonian with Fock state basis ∣*N*_{+1}, *N*_{0}, *N*_{−1}⟩ = {∣*N*/2,0, *N*/2⟩,∣*N*/2 − 1,2, *N*/2 − 1⟩, ⋯∣0, *N*,0⟩} (we will take *h* = 1 for simplicity hereafter). The time evolution of ρ_{0} can be obtained by 〈ρ_{0}〉(*t*) = 〈ψ(*t*)∣ρ_{0}∣ψ(*t*)〉 for distinct *v*. In the dynamics across the transition point, spin excitations from the polar state emerge, which can be reflected by the decrease of 〈ρ_{0}〉(*t*) from one. Let *t _{a}* be the temporal onset of the spin excitations and

*q*=

_{a}*q*(

*t*=

*t*) be the critical quadratic Zeeman energy at which 〈ρ

_{a}_{0}〉(

*t*) begins to change. In Fig. 1B, we show the presence of a power-law scaling for

*q*with respect to the quench rate

_{a}*v*(see the orange squares).

To delve into the reason underlying the presence of the scaling, let us show the presence of impulse and adiabatic evolution regions. It is well known that a metastable phase exists across a first-order quantum phase transition, as shown in Fig. 1A. Intuitively, when we vary the system parameter *q* across zero, the evolving state should stay around this metastable state if the energy gap relative to this metastable state is sufficiently small, suggesting the presence of an impulse region. Yet, when the energy gap becomes sufficiently large, the state cannot jump to the metastable state of the following *q* so that ρ_{0} begins to decrease, suggesting an entrance into an adiabatic region. Specifically, in the impulse region, an evolving state remains frozen in the initial state as time progresses. In other words, if the system remains in the initial state, its maximally occupied level for the evolving state is the same as the maximally occupied level for the initial state. Here, the maximally occupied energy level of the evolving state is defined as the *n*_{max}(*t*)th eigenstate ∣ψ_{nmax}(*q*)⟩ satisfying ∣〈ψ_{nmax}(*q*)∣ψ(*t*)〉∣≥∣〈ψ* _{n}*(

*q*)∣ψ(

*t*)〉∣ for all

*n*with ∣ψ

*(*

_{n}*q*)⟩ being the eigenstate of

, and the maximally occupied energy level for the initial state is defined as the *n*_{smax}th energy level that has the maximal overlap with the initial state, i.e., ∣〈ψ_{nsmax}(*q*)∣ψ_{0}〉∣≥∣ 〈ψ* _{n}*(

*q*)∣ψ

_{0}〉∣ for all

*n*. The latter level coincides with the metastable polar phase with respect to

*q*(see Fig. 1, A and C), which is consistent with the first-order quantum phase transition.

We find that when *q* is varied across zero, the former maximally occupied level *n*_{max}(*t*) rapidly increases by following the latter maximally occupied one *n*_{smax}, as shown in Fig. 1C, suggesting the existence of an impulse region where the state remains frozen. In contrast, when the system leaves this region, the maximally occupied level *n*_{max}(*t*) begins approaching a fixed level, suggesting the presence of an adiabatic evolution. For instance, when *v* = 260 Hz/s, the maximally occupied level *n*_{max}(*t*) follows *n*_{smax} inside the blue region and then converges to around the 2510th level in the long time limit (see the green line in Fig. 1C).

To further demonstrate the existence of impulse and adiabatic regions in the dynamics, we compute the evolution of the probability of atoms occupying the maximally occupied level, i.e., *P _{m}*(

*q*) = ∣〈ψ

_{nmax}(

*q*)∣ψ(

*t*)〉∣

^{2}. As shown in Fig. 1D, we find that the probability changes rapidly near the transition point, consistent with the prediction of an impulse evolution, and remains almost constant in other regions, consistent with the prediction of an adiabatic evolution. In addition, we mark out

*q*as diagonal crosses determined by the numerical simulation, which agrees well with the

_{a}*q*where the system leaves the impulse region and enters the adiabatic region (see Fig. 1, C and D).

The presence of the impulse and adiabatic regions suggests that the scaling law may be accounted for by the KZM. Suppose that at *t* = 0, *q* = *q _{c}* = 0 and the system is in the polar phase.

*q*is then linearly varied by

*q*= −

*vt*. On the basis of the KZM, the critical time when the system begins to respond is determined by τ(

*t*) =

_{a}*t*, where τ(

_{a}*t*) is the relaxation time proportional to 1/Δ

_{a}*E*(

*t*), with Δ

*E*(

*t*) being the energy gap near the transition point. We can also determine the critical time

*t*by 1/∣Δ

_{a}*E*(

*t*)∣ = ∣Δ

*E*(

*t*)/(

*d*Δ

*E*(

*t*)/

*dt*)∣, after which the adiabaticity is restored. If the energy gap Δ

*E*∝ ∣

*q*−

*q*∣

_{c}*with ν being a positive real number, then the critical time is given by*

^{ν}*t*∝

_{a}*v*

^{−ν/(ν + 1)}, yielding

*q*∝

_{a}*v*

^{1/(ν + 1)}. This shows a power-law scaling of

*q*with respect to

_{a}*v*and the scaling exponent is determined by the energy gap. At the second-order phase transition, the relevant energy gap is the gap between the ground state and the first excited state labeled Δ

*E*

_{12}. In our system, this energy gap Δ

*E*

_{12}∝

*q*

^{1/2}is contributed by the Bogoliubov spin excitations as

*q*→ 0 (

*34*). This gives us

*q*∝

_{a}*v*

^{2/3}, consistent with our numerical result

*q*∝

_{a}*v*

^{0.662}(see Fig. 1B). It also tells us that the finite-size effects are very small when

*N*= 1 × 10

^{4}(see the Supplementary Materials for details about finite-size effects). However, at the first-order transition point, the numerical evolution gives us the exponent of 0.740, which is larger than the value predicted by the KZM by more than 10%. In stark contrast, if the energy gap is taken as the gap (dubbed the generalized energy gap) between the maximally occupied energy level, i.e., the

*n*

_{smax}th level and the corresponding first excited state relative to it, i.e., the (

*n*

_{smax}+ 1)th level, we find the exponent of 0.734, which agrees well with our numerical result. This is due to the different energy gap scaling as shown in the inset of Fig. 1B (the scaling exponent for the generalized energy gap is ν = 0.371). For the maximally occupied energy level, while there are two gaps relative to it, one with the next level and the other with the previous level, only the former is relevant because it determines the impulse and adiabatic regions (see Materials and Methods for details). We call this method the generalized KZM. Yet, when we apply the generalized KZM to the second-order quantum phase transition, we find that the result is not as good as the one predicted by the first one, suggesting the difference between the first-order and second-order quantum phase transitions (see the Supplementary Materials for details about the KZM across the second-order quantum phase transition). While the energy gap for a second-order quantum phase transition generically exhibits a power-law scaling near the critical point, whether the power-law scaling of the generalized energy gap for a first-order quantum phase transition is universal is still an open question.

### Experimental results

In experiments, we prepare a sodium BEC in the 3^{2}*S*_{1/2}∣*F* = 1⟩ hyperfine state by evaporation of atoms in an all-optical trap (*31*) and then apply a strong magnetic field gradient to kick the atoms on the *m _{F}* = ± 1 levels out of the optical trap, leaving all atoms on the

*m*= 0 level. After that, we hold the BEC atoms in a uniform magnetic field for 3 s to obtain a polar phase under the quadratic Zeeman energy of

_{F}*q*= 42 Hz induced by the magnetic field. At the end of the holding, we turn on the microwave pulse with a frequency of 1.7701264 GHz (with a detuning of −1500 kHz from ∣

_{B}*F*= 1,

*m*= 0⟩ to ∣

_{F}*F*= 2,

*m*= 0⟩) to change the quadratic Zeeman energy to

_{F}*q*∼ 15 Hz (this time is defined as

_{i}*t*= 0). Subsequently, we linearly vary the quadratic Zeeman energy from

*q*≃ 15 Hz to

_{i}*q*≃ − 38 Hz by ramping up the amplitude of the microwave field. During the entire ramping time, we control the microwave power by a proportional-integral-derivative (PID) system according to the calibration of the quadratic Zeeman energy (see Materials and Methods for details about the

_{f}*q*calibration). As time progresses, we apply the Stern-Gerlach fluorescence imaging to measure ρ

_{0}(

*t*). At each time

*t*, we repeat 15 to 20 measurements to obtain the average 〈ρ

_{0}〉 over the ensemble. Figure 2 displays the observed 〈ρ

_{0}〉 as time evolves for a number of ramping rates

*v*.

*q*is taken as the value when 〈ρ

_{a}_{0}〉 drops below the threshold ρ

_{0c}= 0.98. Evidently,

*q*approaches zero as

_{a}*v*is decreased.

To experimentally measure the power-law scaling, *q _{c}* should be precisely probed. We here use the quench dynamics to identify the error of transitions point (

*31*) in our calibration.

Besides, we also use the result to evaluate the error of *v* (see Materials and Methods for details about the *q* error evaluation).

In Fig. 3 (A and B), we plot the observed ∣*q _{a}* −

*q*∣ with respect to

_{c}*v*in the logarithmic scale, showing the existence of a power-law scaling, i.e., ∣

*q*−

_{a}*q*∣∝

_{c}*v*

^{β}. The fitting of the experimental data gives the exponent of β = 0.728 ± 0.20 when

*c*

_{2}= 25.5 ± 1.5 Hz and β = 0.723 ± 0.25 when

*c*

_{2}= 23.5 ± 0.7 Hz, which are slightly different for different

*c*

_{2}owing to the finite-size effect. The experimental results are also in good agreement with the numerical simulation results: β = 0.739 for the former and β = 0.744 for the latter (see the insets in Fig. 3). We also calculate the scaling law determined by the KZM and find that the exponents predicted by the generalized KZM are 0.733 for (A) and 0.740 for (B), which are closer to the simulation and the experimental results than the exponents of 0.662 and 0.657 predicted by the KZM. In addition, we find that the scaling is not sensitive to the atom loss as we can still observe it in the presence of 18% atom loss (see Materials and Methods for details about atom loss).

The scaling can also be observed when the system is driven from the AFM phase to the polar phase. In experiments, we prepare the initial state of the spinor BEC in a nearly AFM state by shining a π/2-pulse radio frequency radiation to the BEC on the ∣*m _{F}* = 0⟩ level. We then shine a resonant microwave pulse with a frequency of 1.7716264 GHz on the atoms for 300 ms to remove the remaining atoms on the ∣

*m*= 0⟩ level to obtain an AFM state. After that, we suddenly switch off this microwave pulse and switch on another one with a frequency of 1.7701264 GHz. By controlling the amplitude of a microwave field, we are able to linearly vary the quadratic Zeeman energy from around −12 Hz to around 28 Hz. In Fig. 3C, we show the experimentally measured relation between ∣

_{F}*q*−

_{a}*q*∣ and

_{c}*v*, illustrating a power-law scaling with an exponent of 0.724 ± 0.32, which agrees very well with the numerical simulation result of 0.734 and the result of 0.730 predicted by the generalized KZM.