Absorbing phase transition with a continuously varying exponent in a quantum contact process: a neural network approach

Phase transitions in dissipative quantum systems are intriguing because they are induced by the interplay between coherent quantum and incoherent classical fluctuations. Here, we investigate the crossover from a quantum to a classical absorbing phase transition arising in the quantum contact process (QCP). The Lindblad equation contains two parameters, ω and κ, which adjust the contributions of the quantum and classical effects, respectively. We find that in one dimension, when the QCP starts from a homogeneous state with all active sites, there exists a critical line in the region 0 ≤ κ < κ∗ along which the exponent α (which is associated with the density of active sites) decreases continuously from a quantum to the classical directed percolation (DP) value. This behavior suggests that the quantum coherent effect remains to some extent near κ = 0. However, when the QCP in one dimension starts from a heterogeneous state with all inactive sites except for one active site, all the critical exponents have the classical DP values for κ ≥ 0. In two dimensions, anomalous crossover behavior does not occur, and classical DP behavior appears in the entire region of κ ≥ 0 regardless of the initial configuration. Neural network machine learning is used to identify the critical line and determine the correlation length exponent. Numerical simulations using the quantum jump Monte Carlo technique and tensor network method are performed to determine all the other critical exponents of the QCP.

Here, we aim to answer these questions by considering the quantum contact process [33][34][35][36][37][38][39][40] in one dimension (1D-QCP) and two dimension (2D-QCP). In the contact process (CP), each element of the system is in an active or inactive state, and its state changes according to the given CP rules [41][42][43][44][45][46][47]. When all the elements are in the inactive state, the system becomes trapped in a frozen configuration, and the dynamics do not proceed further. This absorbing phase transition of the classical CP (CCP) belongs to the directed percolation (DP) universality class. In addition, the DP transition appears in diverse nonequilibrium systems; however, experimental observation has been difficult, although it has been observed in turbulent liquid crystals [48] and Rydberg atoms [49]. Here, we focus on the recent observation of the 1D-QCP in a dissipative quantum system of Rydberg atoms. * bkahng@snu.ac.kr FIG. 1. Schematic phase diagram of the QCP in the parameter space (κ, ω) in the mean-field limit (inside, d ≥ d c , where d is the spatial dimension, and d c = 3 is the upper critical dimension), in two dimensions (middle, d = 2), and in one dimension (outside, d = 1). Here, "ab" and "ac" represent the absorbing and active phases, respectively. For d ≥ d c , discontinuous (dashed curve) and continuous transitions (solid line) occur, and they meet at a tricritical point (TP). For d = 2, a continuous DP transition occurs over the entire region [0, κ c ]. For d = 1, a continuous DP transition occurs in the region [κ * , κ c ]; however, in the interval [0, κ * ], the exponent α of the density of active sites n(t) ∼ t −α starting from a homogeneous initial state decreases continuously as κ is increased for the QDP values.
The dynamics of the 1D-QCP is described by the Lindblad equation, which consists of a Hamiltonian and dissipative terms. Their contributions to the overall dynamics are adjusted by the model parameters ω (for the coherent quantum effect) and κ (for the incoherent classical dynamics). Thus, the system can exhibit a quantum or classical phase transition in extreme cases. A previous result based on the semiclassical mean-field solution [33] showed that the QCP exhibits a continuous (discontinuous) phase transition when κ is large (small). Thus, a tricritical point exists, as shown in Fig. 1.
The continuous transition belongs to the mean-field DP (MF-DP) universality class, and the critical behavior at the tricritical point belongs to the mean-field tricritical DP (MF-TDP) class, which occurs at a tricritical point of the tricritical CP (TCP) model. Studies using the functional renormalization approach and tensor network approach showed that as the dimension is decreased from the upper critical dimension d c = 3 [34,35], the tricritical point shifts toward the quantum axis with κ = 0 [36,37]. This result indicates that a quantum phase transition would occur only for κ = 0. This behavior may be related to the fact that quantum phase transitions in equilibrium systems occur only at zero temperature. On the basis of this conjecture, a recent numerical study [39] of the 1D-QCP with κ = 0 using the tensor network approach [37][38][39] revealed that the QCP exhibits a continuous transition in a quantum DP (QDP) class. When the dynamics starts from a homogeneous state, that is, all the sites are in the active state, the density of active sites at time t, denoted as n(t), decays as n(t) ∼ t −α at a quantum critical point ω c . The exponent α QDP ≈ 0.32 obtained by the tensor network approach differs from the DP value, α DP ≈ 0.16. However, the other exponents, except for the hyperscaling exponent ν ⊥ of the spatial correlation length, have the DP values. In fact, the discrepancy in ν ⊥ is inevitable because the tensor network approach is not valid for strongly entangled quantum systems. Using the quantum jump Monte Carlo (QJMC) approach, we obtain ν ⊥ as the DP value, as described later. On the other hand, when the QCP starts from a heterogeneous state, that is, all the sites but one are in the inactive state, all the critical exponents have the classical DP values. Thus, the critical behavior of the 1D-QCP with κ = 0 depends on the initial configuration.
For the 1D-QCP, we are interested in how the quantum fluctuations responds to an classical fluctuations, that is, how the critical exponent α changes to the classical DP behavior as the strength κ of the classical fluctuations is increased from zero. We find that there exists an interval [0, κ * ] in which the exponent α decreases continuously from the quantum value, α ≈ 0.32, to the DP value, α ≈ 0.16, as κ is increased from κ = 0 to κ * . The phase diagram for the 1D-QCP is shown in Fig. 1. This result implies that the quantum effect remains to some extent in the region κ = [0, κ * ]. Consequently, the crossover from the 1D-QCP class to the classical DP class occurs in a soft-landing manner.
We remark that the critical behavior of the TCP model in the mean-field limit, equivalent to the MF-TDP class, exhibits a similar behavior to the 1D-QDP class [50,51]. Only the critical exponent α is also deviated from the DP value and the other exponent values remain the same as the DP values. One may guess that quantum fluctuations shift the tricritical point onto the quantum transition point (0, ω c ) [37]. Thus the first-order transition curve disappears and the exponent α QDP has a different value from the corresponding DP value. If this scenario is correct, then the exponent α would have the DP value for any κ > 0. However, we obtain that there exists an interval in which the soft crossover behavior occurs.
We extend our studies to the 2D-QCP. A previous study [34] showed that the 2D-QCP exhibits a discontinuous transition at κ = 0 using the functional renormalization approach [34]. However, we find that the transition is continuous, and it belongs to the classical DP class. The exponent α for homogeneous initial conditions has the DP value. Thus, the intermediate region [0, κ * ] is absent. This suggests that the quantum fluctuations are weaker than the classical fluctuations in two dimensions. For higher dimensions, d ≥ d c = 3, a discontinuous transition occurs at κ = 0.
The critical behaviors of the 1D-QCP and 2D-QCP are obtained as follows. The transition point ω c (κ) for each κ value is determined using the neural network (NN) machine learning algorithm. Then, by applying finite-size scaling (FSS) analysis for different system sizes, the exponent ν ⊥ at ω c (κ) is determined [52][53][54]. At the transition point ω c (κ), the other critical exponents are determined using the QJMC method [55,56] as well as the tensor network method to confirm the results for large systems.
This paper is organized as follows. First, we introduce the 1D-QCP model and specify the classical and quantum limits in Sec. II. The structure of our NN and the optimization scheme are presented in Sec. III A. In Sec. III C, we present the FSS behaviors obtained using the NN, QJMC, and tensor network methods; the universality class; and the crossover behavior. A summary and final remarks are presented in Sec. IV. Additionally, in the Appendices we present a method of experimental realization of QDP behavior in the Appendix A, the derivation of the model parameter κ for the classical dynamics in terms of the experimental parameters of the Rabi frequency Ω and dephasing rate Γ in the Appendix B, the technical details of the QCP in the classical limit using the quantum Monte Carlo method in the Appendix C, the use of the NN approach with different training regions in the Appendix D, the relation δ = α for the classical CP and classical tricritical CP is presented with other critical exponents in the Appendix E, the continuously varying exponent α is discussed more thoroughly in the Appendix F, and we estimate the critical exponent associated with the temporal correlation ν , which is helpful for confirming the critical exponent ν ⊥ = ν /z in the Appendix G.

II. MODEL
We consider a one-dimensional quantum spin chain with a periodic boundary condition, where the active and inactive states of a site represent the up and down spin states, denoted as |↑ and |↓ , respectively. A QCP consists of three incoherent and two coherent processes: i) decay, in which an active site is incoherently inactivated spontaneously at a rate γ; ii) incoherent branching or coagulation, in which an active particle incoherently activates or inactivates an inactive particle at the nearest-neighbor site at a rate κ, respectively; and iii) coherent branching or coagulation, which is a quantum counterpart of process ii) driven by a Hamiltonian at a rate ω. The classical ii) and quantum iii) rules are in competition, which may change the transition behavior.
The time evolution of the density matrixρ is described by the Lindblad equation, which consists of the Hamiltonian and dissipative terms [57]: The HamiltonianĤ S , which governs the branching and coagulation processes and represents coherent interactions, is expressed asĤ Here,σ i l denotes the Pauli matrix, where the superscript and subscript represent the spin axis and site index, respectively, andn l is the number operator for an up spin at the lth site. The Lindblad decay, branching, and coagulation operators are given byL respectively.σ + andσ − are the raising and lowering operators of the spin at site , which are defined in terms of the spin basis asσ + = |↑ ↓| andσ − = |↓ ↑|, respectively. In addition,n =σ +σ− andσ x =σ + +σ − are the number operator and spin flip operator, respectively. We rescale the time and quantum control parameters ω and κ, respectively, in units of γ; therefore, we set γ = 1. Quantum branching and coagulation occur at a rate ω, and the corresponding classical processes occur at a rate κ. When ω → 0, the model is reduced to the CCP, which belongs to the DP class. Here, we first consider the pure quantum limit κ → 0, but with finite ω. The opposite limit, ω → 0 with finite κ, and the case of finite ω and κ are discussed in the Appendix C.
When ω is small, inactive particles become more abundant with time, and eventually the system is fully occupied by inactive particles. Thus, the system is no longer dynamic and falls into an absorbing state, which is represented bŷ ρ ab = |↓ · · · ↓ ↓ · · · ↓|. When ω is large, the system remains in an active state with a finite density of active particles [ Fig. 2(a)]. Thus, the QCP exhibits a phase transition from an active to an absorbing state as ω is decreased.

A. NN approach
The NN approach has recently been used as a powerful tool [58,59] for understanding phase transitions in classical systems [52], which exhibit patterns involving many components. Each component has one of two values, for instance, the up and down spin states in ferromagnetic systems. By contrast, each component of a quantum system has a real value, and thus the patterns are much more complex. Nevertheless, the NN approach has reportedly been used successfully to determine the transition points of closed quantum systems on the basis of simulation data [60][61][62] and experimental images [63,64]. In addition, an unsupervised NN approach was used to generate the configurations of quantum dissipative systems in the steady state [65][66][67][68][69]. This approach was efficient, because it uses less computing resources to generate configurations compared to conventional simulations. However, determining the transition point accurately in dissipative quantum systems requires large system sizes, and thus the simulations require considerable computational resources. Here, we use the NN approach as an alternative method and obtain a transition point that is sufficiently accurate for investigation of the critical behavior of the 1D-QCP at the transition point.
For classical systems, the transition point of a continuous absorbing transition is usually indicated by the presence of power-law behavior of the order parameter with respect to time [41,51]. Consequently, a large system size is required to identify the transition point. Accurately identifying the transition point using QJMC simulations of the QCP is even more difficult and is thus a challenging problem. To overcome this difficulty, we note that the system is in the absorbing state for ω ω c and in the active state for ω ω c . Combining this observation with a recently proposed NN supervised learning concept, we identify the transition point as follows.
To implement the NN approach, we first organize a dataset of the occupation probability of site , which is denoted as Using the QJMC method, we generate a steady-state configuration and obtain the occupation probabilities of each site, {p }. We collect 5 × 10 3 configurations in ω ∈ [0, 12] at ∆ω = 0.04 intervals. To prepare the training dataset for supervised learning, we label the configurations using one-hot encoding [70], where the absorbing state (ω ∈ [0, 4]) is encoded as (0, 1), and the active state (ω ∈ [8,12]) is encoded as (1, 0) [see shaded regions in Fig. 3(a)].
After we collect the snapshots, we try to train the NN. The objective of the learning procedure is to optimize the NN to adjust the weights of the connections between neural units to achieve the variational minimization of a properly defined cost function. To this end, we construct the hidden layers of the NN, including one-dimensional convolutional layers, batch normalization layers [71], and fully connected layers, as shown in Fig. 4. We employ the framework of tensorflow [72] and use ReLU and tanh for the activation function in the hidden layer. Two neurons in the output layer are used, and a softmax function is used as the activation function in the output layer. We employ the cross-entropy or the meansquare error function as the cost (error) function of the NN, which is then optimized using Adam [73] or RMSProp. We vary the architecture and optimization algorithms in various ways. Regardless of these changes, the well-trained machines produce consistent results. To check the sensitivity to the positions of the left and right boundaries, the NN approach with different training regions is described, and its results are presented in the Appendix D. When the NN is well-trained with the labeled training dataset in the two regions, we obtain the outputs for the entire ω region.
Next, using the obtained transition points ω c (N) for selected system sizes, we perform FSS analysis and identify the transition point in the thermodynamic limit, ω c . We also determine the correlation length exponent ν ⊥ . Next, we determine the other critical exponents by performing extensive QJMC simulations up to system size N = 20 and applying the tensor network method to a large system (N = 80) at ω c . Specifically, the tensor network method is implemented by the QJMC method with the matrix product states and timeevolving block decimation.

B. FSS analysis of 1D-QCP
In Fig. 3(a), the two outputs of the NN indicate the probabilities that the system will fall into the absorbing state and remain in the active state, respectively. The crossing point of these outputs indicates a transition point ω c (N) for a given system size N [ Fig. 3(a)]. Several studies [52][53][54] have shown that the predictability exhibits FSS behavior. Using the obtained ω c (N) for different system sizes, we determine ω c in the thermodynamic limit by plotting ω c − ω c (N) versus N [ Fig. 3(b)], which is expected to behave as ω c − ω c (N) ∼ N −1/ν ⊥ . Indeed, the plot exhibits power-law behavior when an appropriate value of ω c is chosen, and the critical exponent ν ⊥ is obtained from the slope of the power-law curve. We obtain ω c ≈ 6.04 and ν ⊥ = 1.06 ± 0.04; the latter is consistent with the value of ν ⊥ ≈ 1.096 for the DP class in one dimension but differs from the value of ν ⊥ ≈ 0.5 ± 0.2 obtained using the tensor network approach. The scaling plot is drawn in the form of the output versus (ω − ω c )N 1/ν ⊥ for different N [ Fig. 3(c)]. The data for different system sizes seem to collapse.
Next, we measure the values of the other critical exponents using the QJMC method in the critical region around ω c with the tensor network method for larger system sizes. First, we consider an initial state in which a single active seed is present at = 0, and the remaining sites are inactive. This configuration is expressed asρ(0) =σ + 0ρ abσ − 0 . We measure the following quantities: i) the survival probability, that is, the probability that the system does not fall into an absorbing state, Tr[ρ(t)n ]; iii) the mean square distance of the active sites from the origin, and v) the seed-site density over surviving runs, ρ d,s (t) = ρ d /P(t). At the transition point, these quantities exhibit the following power-law behaviors: We estimate the exponents δ+δ , η, δ , z, and δ by direct measurement of the slopes in the double-logarithmic plots, as shown in Fig. 5. We estimate the exponent z using the data collapse technique. For instance, for the survival probability P(t), we plot P(t)t δ versus tN −z for different system sizes N. We determine z as the value at which the data for different system sizes collapse onto a single curve. The values of the critical exponents are in good agreement with the DP values within the error bars (Table I). Second, we take a homogeneous initial state in which the entire system is occupied by active sites at t = 0, which is expressed asρ(0) = |↑ · · · ↑ ↑ · · · ↑|. From this initial state, we measure vi) the density n(t) of active sites at time t averaged over all runs. This quantity is formulated as n(t) = ( Tr[ρ(t)n ])/N. We find that n(t) exhibits power-law decay as n(t) ∼ t −α with exponent α = 0.32 ± 0.01, as shown in Fig. 6(a). This value is consistent with the result obtained by applying the tensor network approach; however, it is not con-  sistent with the corresponding DP value, which was estimated as α DP = 0.16. Therefore, the 1D-QCP for κ = 0 creates a new type of universal behavior.
We note that ρ d (t) and n(t) are actually the same quantity even though they emerge from different initial states [74]. They exhibit the same critical behaviors in the CP class (See the Appendix E), but they exhibit different critical behaviors for the 1D-QCP. This behavior is unusual because the universality class is independent of the initial state according to the theory of critical phenomena. To understand the underlying mechanism, we increase the control parameter κ from zero to 0.6 in steps of 0.2 and explore the behavior of n(t) at each ω c (κ) [ Fig. 6(b)]. We find that the value of α decreases continuously from 0.32 for κ = 0 to α = 0.16 for κ = 0.6. Furthermore, we apply the tensor network method based on the matrix product states and time-evolving block decimation to confirm the results for a large system [ Fig. 6(c)]. Using FSS analysis, we determine the exponent α for each κ (See the Appendix F). These results suggest that α decreases continuously as κ is increased and reaches the DP value at κ * ≈ 0.58 (Table II).

C. FSS analysis of 2D-QCP
We investigate the critical behavior of the 2D-QCP. At κ = 0, we find that the 2D-QCP exhibits a continuous absorbingstate phase transition. By a process similar to that used for the 1D-QCP, we obtain the critical exponents of the 2D-QCP using QJMC simulations. For κ = 0, we obtain ω c ≈ 0.94. At this ω c , we find that n(t) exhibits power-law decay as n(t) ∼ t −α with α = 0.45 ± 0.03, as shown in Fig. 7(a). This value is in agreement with the corresponding DP value in two dimensions. In Fig. 7(b), the exponent ν is determined by rescaling the plot of n(t)t α versus t(ω c − ω) ν for different ω as ν = 1.30. Thus, ν ⊥ = ν /z 0.74, and the critical exponents obtained for a homogeneous initial state are in good agreement with the DP values within the error bars.
Next, we estimate the exponents η and δ by direct measurement of the slopes of the double-logarithmic plots, as shown in Fig. 8. We estimate the exponent z using the data collapse technique. For instance, for the survival probability P(t), we plot P(t)t δ versus tN −z for different system sizes N. We determine z as the value at which the data for different system sizes collapse onto a single curve. For the 2D-QCP, we find that rapidity-reversal symmetry holds because α = δ . The val- ues of the critical exponents are in good agreement with the DP values within the error bars. Therefore, the 2D-QCP for κ = 0 belongs to the DP class. We conclude that the quantum coherent effect is irrelevant in two dimensions.

IV. SUMMARY AND DISCUSSION
We investigated the 1D-QCP and 2D-QCP as prototypical examples of nonequilibrium absorbing phase transitions in dissipative quantum systems. A phase diagram was obtained (Fig. 1) in the parameter space (κ, ω), where these variables represent the contributions of the classical and quantum effects, respectively. When the 1D-QCP starts from a homogeneous state, the transition curve between the absorbing and In the quantum region, the critical exponent α, which is associated with the density of active sites n(t), decreases continuously as κ is increased from the quantum value to the classical DP value, as presented in Table II. Thus, the crossover from quantum to classical behavior proceeds gradually. When the 1D-QCP starts from a heterogeneous state, this soft-landing crossover does not occur. For the 2D-QCP, we find a continuous transition at κ = 0 in the DP class, which is inconsistent with the prediction by the functional renormalization group approach [34]. Interestingly, in the mean-field solution [33][34][35], the transition in the region near κ = 0 is discontinuous. This discontinuous transition changes to a continuous transition at a tricritical point, which belongs to the MF-TDP class. In the MF-TDP class, the critical exponent α MF−TDP has a different value from the corresponding DP value. However, there is no crossover region. Thus, there is a discrepancy in the crossover behavior between the 1D-QDP and the mean-field QDP. It would be interesting to investigate underlying mechanism for this discrepancy, which is left for the next problem.
The NN approach we used was applied to a dataset obtained by QJMC simulations. We believe that this approach can be applied directly to snapshots obtained in cold Rydberg atom experiments. Furthermore, our approach will be used in the near future for other dissipative systems where a transition point is hardly determined. Examples include the dissipative transverse-field Ising model, dissipative XYZ model, and dissipative anisotropic Heisenberg model.  Here, we concern experimental realization of QCP, specifically on the antiblockade mechanism matching to the branching and coagulation processes. It is known theoretically that the quantum limit κ = 0 can be reached by controlling the dephasing strength. This is realized by controlling the laser linewidths and residual Doppler broadening. The classical DP limit (strong dephasing limit) of QCP is experimentally verified in Ref. [49]; however, it is intrinsically unreachable to the quantum limit with zero dephasing, because setting the laser linewidth to zero is challenging. Thus, achieving the experimental realization of the quantum limit has yet to be resolved. It is analogous to the situation that zero temperature is unreachable precisely experimentally in the transverse Ising model. Here, we propose that using our result, a quantum effect can be found at κ 0 but in the interval [0, κ * ]. Below we derive the formula to reach this interval in terms of experimental quantities.
An essential factor in this experiment is the antiblockade effect, in which an inactive (active) spin is activated (inactivated) by detuning the excitation energy so that it is comparable to the interaction energy with the active spin of the nearest neighbor. Moreover, when the dephasing noise is sufficiently strong, the quantum coherences are suppressed, and the dynamics may be reduced to those of the CCP process [75]. Otherwise, quantum coherence is effective, and a combination of coherent and incoherent CPs can be realized [33,34]. Thus, by controlling the ratio of the Rabi frequency and the dephasing rate, one may observe crossover behavior from a quantum to a classical nonequilibrium phase transition. It has been shown that the classical parameter is related to κ = 4Ω 2 /Γ, where Ω is the Rabi frequency, and Γ is the dephasing rate [75]. Thus, for Ω < √ κ * Γ/2, QDP behavior may appear. This region is a quantum critical region where quantum fluctuations play a role in the universal behavior.
Let us consider the Rydberg atom represented by two spin states, |↓ and |↑ , which are the eigenstates of the Pauli matrix in the z direction, where the up spin state indicates Rydberg excitation. We start with the Hamiltonian in the rotating-wave approximation of N coupled Rydberg atoms on a lattice, as follows:Ĥ where Ω and ∆ are the Rabi frequency induced by the external laser and the detuning strength, respectively. The third term on the r.h.s. of Eq. (A1) describes the interaction between the up spins, where V m is a power-law decay function of distance: where C p is the dispersion coefficient [76], and x l is the position of the lth site. Note that p characterizes the interaction [76]; p = 3 for the dipole interaction associated with dorbital excitation, and p = 6 for the van der Waals interaction with s-orbital excitation. For p = 6 in Eq. (A2), the interaction decays rapidly with distance; thus, the dynamics is effectively dominated by the nearest-neighbor interaction [77].
In the limit ∆ Ω, |V m |, the Rabi oscillation is suppressed because of the large energy gap between the up and down spin states, which suggests that a spin configuration, for example, |a = | · · · ↓↓↑ · · · , may be approximately an eigenstate of H R . Therefore, given an initial spin configuration, there is no fluctuation in time. However, if V m = −∆ for the nearestneighbor pair and m, and the th spin is up, the mth spin can fluctuate with Ω. This is the so-called antiblockade mechanism, in which an excited atom facilitates Rabi oscillation at the nearest-neighbor atom. The effective Hamiltonian can be described in one dimension: where the projection operator is given byP =n −1 +n +1 − 2n −1n +1 . We are interested in the region near the transition point, where the low-density limit leads toP ≈n −1 +n +1 . Thus, the quantum Hamiltonian in the main text is obtained. Additionally, the decay in the main text is implemented by radiative decay from |↑ to |↓ using a zero-temperature heat bath.
On the other hand, the classical limit is obtained via the strong dephasing rate, which is denoted as Γ, as extensively discussed in Ref. [75]. In the strong dephasing limit, it is known that the coherent dynamics can be neglected; thus, the effective Hamiltonian in Eq. (A3) effectively becomes classical and is given by the Lindblad operators in the main text. The strong dephasing limit is derived using superoperator formalism in the next section. The resulting relation is given as κ = 4Ω 2 /Γ.
When we combine the quantum and classical settings, the relationship between the experimental realization and QCP is given by ω = Ω 1 and κ = 4Ω 2 2 /Γ, where Ω 1 and Ω 2 are the Rabi frequencies of the independent laser source [33,34]. Thus, for Ω 1 = ω c and Ω 2 < √ κ * Γ/2, QDP behavior appears. This region is a quantum critical region where quantum fluctuations play a role in the universal behavior. More details are provided in the next section. The critical behavior of the CCP can be realized using Rydberg atom experiments. Here we show that the transition rate κ of the branching and coagulation processes in the classical limit can be obtained in terms of the experimental parameters used in the Rydberg atom system.

Lindblad equation for Rydberg gases
To describe the Rydberg gases in open quantum systems, we employ the Lindblad equation as follows: We have split H R into two parts, and the Lindblad dissipator D (i) is given by Here, the Lindblad operators defined on each site,L (1) and L (2) , represent decay from |↑ to |↓ by a zero-temperature heat bath and the dephasing processes, respectively, where the detailed forms are given bŷ By using the representation of the operator with the spin configurations as the basis, for example, ρ ab ≡ a|ρ|b , where a and b denote one of the 2 N configurations of N spins, that is, |a =|↓↑ · · · ↓ , each term on the r.h.s. of Eq. (B1) can be evaluated as follows. BecauseĤ 0 is the diagonal matrix for the basis, the first term reads a|[Ĥ 0 ,ρ]|b = ∆ ( a|n |a − b|n |b ) + m V m 2 ( a|n n m |a − b|n n m |b ) ρ ab , which becomes zero when a = b. By contrast,Ĥ Ω yields the transition between the states by flipping a single spin, and the second term in Eq. (B1) can be written as where the diagonal term a = b is given by Note that the diagonal part of Eq. (B6) contains only the off-diagonal contributions of the density matrix because a|σ ± |a = 0, which will be used to derive the rate equation.
Similarly, we obtain the representations of the Lindblad operators; for the decay operator, where the diagonal term a = b is given by For the dephasing operator, which becomes zero when a = b. If the Lindblad equation consists of only this dephasing dissipator, the density operator ρ ab in the long time limit becomes the diagonal matrix; therefore, D (2) is called the dephasing operator. It is known that in the limit Γ Ω, γ, the coherent dynamics can be neglected, and thus the Lindblad equation [Eq. (B1)] effectively reduces to the classical rate equation [75]. We will briefly show the procedure below.

Derivation of transition rate in classical limit
For convenience, we introduce the superoperator [78] by mapping the density operator to the density vector, ρ ab → ρ α , where a and b are those of the spin configurations as defined above, and α is the corresponding vector index. Then we rewrite the Lindblad equation [Eq. (B1)] in terms of the vectorized density operator ρ as where the superoperatorĤ defined in a 4 N ×4 N complex space is given by because Eq. (B1) is the linear equation of ρ α . We decompose the density vector to two parts as ρ = µ ⊕ ν, where µ (which belongs to the 2 N -dimensional space M) and ν (which belongs to the 4 N − 2 N -dimensional space N) are defined by arranging ρ ab as follows: µ = (ρ 0 , ρ 1 , · · · , ρ α , · · · , ρ 2 N −1 ρ α =ρ aa ) T ∈ M , ν = (ρ 2 N , ρ 2 N +1 , · · · , ρ α , · · · ρ 4 N −1 Here, the first 2 N components of the vectors correspond to the diagonal components of the density matrix, and the others are the off-diagonal terms of the density matrix. From these decomposed vectors and Eqs. (B4)-(B9),Ĥ can be decomposed using the block matrix = Ω a a|σ x |a δ β,a a − a |σ x |a δ β,aa , (B16) where each (·) α denotes the matrix element defined in Eqs. (B4)-(B9) corresponding to the vector index α. Consequently, one can rewrite Eq. (B10) in terms of µ and ν: In the strong dephasing limit, whereĤ (2) may dominate the off-diagonal dynamics in Eq. (B20), the solution of Eq. (B20) becomes an approximately exponentially decaying function in time with a time scale of 1/Γ (∆ is also a large parameter, butĤ 0 only induces oscillation). Therefore, the full dynamics can be effectively reduced to the diagonal dynamics of µ at a slower time scale than 1/Γ. Inserting the solution of Eq. (B20), which is given by into Eq. (B19) with ν(0) = 0, we obtain Because Γ, ∆ γ, Ω, we can expand the exponential function about the small parameter using the Zassenhaus formula [79]; it then becomes We remark again that the time scale for the dynamics of µ is much larger than 1/Γ, which results in the replacement µ(t ) → µ(t) in the integrand of Eq. (B23), as in the Markov approximation [57]. In addition, we let the lower limit of the integral go to negative infinity because of the rapid decay of the exponential function. Then, rearranging the time integral in Eq. (B23) as (t − t ) → t , we obtain Now, one can see that the slow dynamics in Eq. (B24) is given only in terms of µ α , and thus it can be rewritten using the diagonal elements ρ aa of the density operator. The first term on the r.h.s. of Eq. (B24) is given by the diagonal element of Eq. (B18), which is composed only of diagonal elements µ α or ρ aa as follows: Next, becauseĤ 0 andĤ (2) are diagonal matrices, the αth component of the second term is given by Using Eqs. (B16) and (B17), one can rewrite Eq. (B25) by changing α → aa: ab,ab ,+H (2) ab,ab t + e −i H 0 where we have used , a|σ x |b b|σ x |a = a|σ x |b 2 , which becomes zero when and are different. Note that the transition rate Λ a,b = Λ b,a induced by the interaction is nonzero only when the spin configuration b is given by |b = |a 1 , where {|a 1 } is generated from |a by flipping a single spin. For example, if |a = |↑↓↓ is given for N = 3, a set {|a 1 } can be defined as {|↓↓↓ , |↑↑↓ , |↑↓↑ }. Now, using H (2) aa 1 ,aa 1 = −iΓ/2 for any a 1 and H 0 aa 1 ,aa 1 = a|Ĥ 0 |a − a 1 |Ĥ 0 |a 1 , we obtain the nonzero transition rate Λ a,a 1 in Eq. (B26), as follows: From Eqs. (B7) and (B26), therefore, one can see that the rate equation Eq. (B24) up to the second order is written as where the transition rate W a,b reads First, we consider only the nearest-neighbor interaction, V m = V 0 , for the nearest-neighbor pair ( , m); otherwise, V m = 0. Defining a as |a 1 =σ x |a , we obtain the transition rate Λ a,a in Eq. (B27) as whereP denotes the number of nearest neighbors for which the site is in the up state, that is,P ≡ m∈n.n.( )nm . Setting V 0 = −∆ and ∆ Γ, one can see that Eq. (B30) is approximately zero, except for a|P l |a = 1. In the low-density limit, where the number of up spins per site is vanishingly small (n = tr ρ n /N 1), the configurations having a small number of up spins make a major contribution toρ. Then one can assume that a|P |a for the major configuration is generally 0 or 1; thus, we obtain the following approximation: Using |a =σ x |a 1 and considering general configurations b, the transition rate Λ a,b becomes which is equivalent to the branching and coagulation processes in the ordinary CP model. Here we used a|n mσ + |b + a|n mσ − |b = | a|n mσ + |b | 2 + | a|n mσ − |b | 2 in Eq. (B32). We expect that this approximation in the low-density limit may be valid near the absorbing transition point, where the order parameter n is small.
Finally, we briefly present the diagonal component of the Lindblad equation with the Lindblad operatorL corresponding to Eq. (B28): where the transition rate is given by W a,b = | a|L |b | 2 . Thus, from Eqs. (B29) and (B32), the three Lindblad operators are obtained with κ = 4Ω 2 /Γ. In this section, we consider ω → 0. When κ γ, inactive particles become more abundant with time, and ultimately the system is fully occupied by inactive particles. Thus, the system is no longer dynamic and falls into an absorbing state. When κ γ, the system remains in an active state with a finite density of active particles. Thus, the CCP exhibits a phase transition from an active to an absorbing state as κ is decreased.
The critical exponents of the 1D-QCP were obtained using FSS of the data obtained using the QJMC method as described in the main text. To check the validity of the FSS for a small system, we consider the 1D-CCP, where κ is finite, and ω = 0 [see Eqs. (2)(3)(4)(5) in the main text]. At the critical point, we perform FSS of the 1D-CCP using the QJMC method. The observables are those defined in the main text.
First, we obtain the exponents δ + δ , η, δ , z, δ, and α directly by measuring the slopes of the double-logarithmic plots shown in Figs. 9 and 10. Then we collapse the data by using the obtained exponents to compute the dynamic exponent z. Specifically, we plot ρ d t δ+δ versus tN −z in Fig. 9(a), N a t −η versus tN −z in Fig. 9(b), and P(t)t −δ versus tN −z in Fig. 9(c) for different system sizes N. We measure the exponent z directly using the plot of R 2 (t) versus t in Fig. 9(d). In the CCP, we can classify the surviving runs, and thus we measure the exponent −δ directly using the plot of ρ d,s (t) versus t in Fig. 9(e). Next, we plot n(t)t −α versus tN −z in Fig. 10(a) for different system sizes N. The exponent ν is obtained from the rescaling plot of n(t)t α versus t(ω c − ω) ν for different ω in Fig. 10(b).
The critical exponents are thus obtained as δ + δ = 0.32 ± 0.01, η = 0.31 ± 0.02, δ = 0.16 ± 0.01, δ = 0.16 ± 0.02, z = 1.58 ± 0.03, and α = 0.16 ± 0.01. Note that δ = α. In addition, α = δ , indicating that rapidity-reversal symmetry holds. All critical exponents are in good agreement with the DP values within the error bars. Thus, we verified that the critical exponents in the CCP can be successfully obtained using the QJMC method with the same system size as in the main text.  For ω = 4.00, an exponentially decaying curve is observed. By contrast, for ω = 8.00, a stationary state converges to a finite density. At the critical point ω = 6.04, it exhibits power-law behavior. The system size is N = 20.
Appendix D: Critical behavior obtained by neural network approach with different training regions For supervised learning, it is advantageous to take a narrower test region [white region in Fig. 11(a)] because more information can be obtained in the training region. However, if the test region is too narrow to include the crossing point, the crossing point of the outputs would not be the critical point. Consequently, it is desirable to consider a test region of an appropriate size.
We took the left boundary ω = 4 in the main text because this is the value at which the order parameter n(t) decays exponentially, that is, at which the system is in the subcritical region, as shown in Fig. 12. This result was obtained using the QJMC method. However, the boundary ω = 8 was considered because n(t) behaves as it does in the supercritical state.
To boundaries, we also considered a test region of (3 ≤ ω ≤ 9) and then estimated the transition point ω c in the thermodynamic limit and the value of the exponent ν ⊥ . As shown in Fig. 11, we obtained the same values of ω c and ν ⊥ .