Bath-induced collective phenomena on superconducting qubits: synchronization, subradiance, and entanglement generation

A common environment acting on a pair of qubits gives rise to a plethora of different phenomena, such as the generation of qubit-qubit entanglement, quantum synchronization and subradiance. Here we define time-independent figures of merit for entanglement generation, quantum synchronization and subradiance, and perform an extensive analytical and numerical study of their dependence on model parameters. We also address a recently proposed measure of the collectiveness of the dynamics driven by the bath, and find that it almost perfectly witnesses the behavior of entanglement generation. Our results show that synchronization and subradiance can be employed as reliable local signatures of an entangling common-bath in a general scenario. Finally, we propose an experimental implementation of the model based on two transmon qubits capacitively coupled to a common resistor, which provides a versatile quantum simulation platform of the open system in any regime.

A common environment acting on a pair of qubits gives rise to a plethora of different phenomena, such as the generation of qubit-qubit entanglement, quantum synchronization and subradiance. Here we define timeindependent figures of merit for entanglement generation, quantum synchronization and subradiance, and perform an extensive analytical and numerical study of their dependence on model parameters. We also address a recently proposed measure of the collectiveness of the dynamics driven by the bath, and find that it almost perfectly witnesses the behavior of entanglement generation. Our results show that synchronization and subradiance can be employed as reliable local signatures of an entangling common-bath in a general scenario. Finally, we propose an experimental implementation of the model based on two transmon qubits capacitively coupled to a common resistor, which provides a versatile quantum simulation platform of the open system in any regime.

I. INTRODUCTION
Common baths inducing dissipation on more than one qubit at a time are known to give rise to different collective phenomena [1]. Much studied has been superradiance, i.e. the enhancement of the emission power of atoms jointly dissipating into a common environment [2]. However, the physics of this system is even richer, with the emergence of additional effects such as subradiance -a complementary effect of superradiance [3,4], qubit-qubit entanglement generation [5], and quantum synchronization [6]. The goal of this paper is to map out the parameter regimes of these collective phenomena and to characterize them by timeindependent figures of merit. Specifically, the questions we want to address are the following: is there a one-to-one relation between entanglement generation, quantum synchronization and subradiance between a pair of qubits? Is it possible to find a regime where the bath induces one of the above collective phenomenon, but is not able to generate a different one? If we know that, for instance, by tuning a certain model parameter we enhance the generation of quantum synchronization and/or subradiance, can we say the same for entanglement? If the answer to the latter question is affirmative, then we may employ quantum synchronization and subradiance as local signatures for the presence of an entangling bath, without the resourceintensive complete state tomography. Mathematically speaking, we aim to investigate the collective features of the quantum map (dynamical semigroup) describing the open dynamics, or, equivalently, of the corresponding Markovian master equation [1]. We provide specific experimental prescriptions for measuring each of these figures of merit: one would need to start from different initial states, e.g. Bell singlet state for subradiance, maximally overlapped state for entanglement, etc.
Note that in addition to quantum synchronization and entanglement generation, we have chosen to study subradiance instead of superradiance between a pair of qubits. Despite its fragility against local noise, subradiance, which manifests itself through the emergence of a slowly-decaying collective mode in the atomic emission, has been reported in different experiments, see for instance Refs. [7][8][9][10]. Furthermore, while the presence of a slowly-decaying collective mode is a patent manifestation of subradiance, in the case of only two qubits superradiance does not exhibit a clear enhancement of the emission power [2]. Therefore, a figure of merit for subradiance is more suitable for the quantitative and comparative study we will carry on in this work. Finally, we will also compare the results with a measure of the "collectiveness" of the dynamics [11] in order to understand whether or not there is a clear, monotonic relation between the strength of the above-mentioned collective phenomena and the spatial correlations induced by the bath. This measure has been tested in a recent experiment [12].
In order to apply our findings to a physical platform of interest for quantum technologies, throughout the paper we will consider a model of two superconducting qubits and, in particular, we will propose to test our predictions through an experimental implementation involving transmon qubits and resistors. Superconducting qubits such as transmons are among the best candidates in the race for the creation of a quantum computer [13][14][15][16] and configurations such as the ones discussed here represent a versatile simulation platform for quantum thermodynamics [17][18][19]. Within this challenge, one of the major issues is the unavoidable presence of noise [15,[20][21][22][23], which induces dissipation and decoherence and breaks all the most fundamental quantum features, such as superposition or entanglement. Here the study of correlated dissipative noise induced by a common bath, which is the topic of our paper, is of particular importance. This type of noise appears inevitably due to the coupling to the electromagnetic modes of the nearby waveguides, resonators, and other circuit elements [24][25][26]. Its investigation has being receiving more and more interest in the recent past; for example, the two-qubit spectroscopy of spatiotemporally correlated noise between two superconducting qubits has been recently proposed [27]. The interest in correlated noise is motivated e.g. by its detrimental effects on the performance of quantum error correction protocols [28], which must be correspondingly modified [11,[29][30][31][32][33][34][35][36][37][38][39][40][41][42]. A key factor that for this purpose needs to be taken into account is the generation of qubit-qubit entanglement due to the action of the common bath during the evolution [35,43,44]. The use of synchronization and/or subradiance as local signatures of entanglement generation may represent a useful tool for the analysis of these scenarios.
Following this proposal, our findings suggest a novel application of quantum synchronization for quantum tasks. Recently, transient synchronization was reported to allow for probing the spectral density of a single qubit immersed in a dissipative environment [45]. Different forms of synchronization have also been reported as beneficial for atomic clocks operation [46,47]. Finally, the superconducting platform we will discuss could also be implemented as one of the first experimental realizations of quantum synchronization [48][49][50].
The paper is structured as follows: Sec. II presents the model and discusses the system dynamics starting from the master equation of the two qubits, while we introduce and analyse the measures associated to each phenomenon in Sec. III. We discuss the results in Sec. IV, addressing some analytical limits in Sec. IV A, and drawing a numerical comparison between the phenomena for several scenarios in Sec. IV B. We propose the superconducting experimental implementation of our model in Sec. V, where we also provide the simulation of the system dynamics in a concrete case. Finally, we draw some concluding remarks in Sec. VI.

II. MODEL
In order to characterize the collective dissipative effects, we consider two superconducting transmon qubits [24] embedded in a common bath. Specifically, in the experimental proposal of Sec. V the bath will be a resistor jointly coupled to both transmon qubits. Such a resistor can be modeled as an infinite collection of bosonic modes, each of them corresponding to an independent fictitious LC circuit, which induce dissipation as described by the Caldeira-Leggett model [51]. Remarkably, the circuit Hamiltonian of such a system depends only on the circuit topology, i.e. on how circuit elements are connected between each other, and, for instance, the qubit-qubit distance does not play any role. In other schemes, the role of a common bath may be played by an external electromagnetic environment whose wavelength is larger than the qubit-qubit distance [24,52], by a structured phononic environment [53] or by a resonator to which both qubits are dispersively coupled [27]. The qubits frequencies can be different, reflecting either the inherent variability due to fabrication imperfections or the intentional detuning achieved by applying bias magnetic fields. In order to identify the effects induced by a collective environment, we focus on a set-up in which the qubits are not directly interacting, as this could mask the collective origin of dissipative effects. The Hamiltonian of the overall system reads: where we have introduced the system Hamiltonian H S = ω1 2 σ z 1 + ω2 2 σ z 2 , with eigenvectors |gg , |ge , |eg , |ee . Here, |e and |g are respectively the excited and the ground state of a qubit, ω 1 and ω 2 are the frequencies of respectively the first and second qubit, while the dimensionless coefficients g 1 and g 2 express the weights of the dissipative coupling of each qubit. The parameter µ is the coupling constant in the units of energy, while each f k is a real dimensionless number representing the strength of the coupling of the k−th mode of the bath to the qubits. The numbers {f k } k define the spectral density of the common bath [1] (see Appendix A), which we will take as Ohmic. This choice is well-justified if the bath is a common resistor, as in the experimental proposal of Sec. V: the properties of resistor-induced noise in circuit QED are well-understood and, for instance, the spectral density is known to be Ohmic [17,51], leading to the correct Johnson-Nyquist spectrum. In other words, the noise generated by a common resistor is the standard thermal noise in the quantum regime [51]. Finally, the weak-coupling condition reads ω 1 , ω 2 µ/ . Motivated by the experimental decay observed in transmon qubits (see e.g. Ref. [23]), we assume that the system follows a Born-Markov master equation [1], and also allow for an additional phenomenological local bath on each qubit, inducing local dissipation characterized by the time T 1 . More details can be found in Appendix A. With these prescriptions, the state of the two-qubit system ρ S obeys the following master equation under the partial secular approximation [54] d dt where H LS is the Lamb-shift Hamiltonian which reads while D is the dissipator defined as: The master equation (2) is valid under the weak-coupling assumption ω 1 , ω 2 µ/ and under the assumption of Markovian bath, which is in general fulfilled when the spectral density is Ohmic [1,54]. Note that we have employed a partial secular approximation instead of a full one [54], i.e. we keep slowly rotating terms driven by the frequency detuning ω 1 − ω 2 in the derivation of Eq. (2). This allows us to treat scenarios with very small detuning, where the standard master equation with full secular approximation fails (for further details we refer the reader to Ref. [54]). The coefficients of the master equation (2) contain both a contribution coming from the common bath with Ohmic spectral density, introduced in Eq. (1), and a contribution proportional to 1/T 1 , i.e. due to the phenomenological local bath acting on each qubit. Their form can be found in Appendix A.
The solution of the master equation is obtained by finding the eigenvalues together with the right and the left eigenvectors of the Liouvillian L. As discussed in Appendix B, the master equation is in partial secular approximation and therefore displays the phase-covariance symmetry on the superoperator level [55]. This means that the Liouvillian superoperator commutes with the number superoperator, defined as N = [N, · ], where N counts the total number of qubit excitations. Thanks to this symmetry, we can divide the Liouvillian superoperator into five blocks [55], which we label using the subscript d such that d = −2, −1, 0, 1, 2: L = 2 d=−2 L d . The blocks have respectively dimension 1, 4 × 4, 6 × 6, 4 × 4, 1, and two of them can be trivially obtained from the others: We name the eigenvalues of the block L d as {λ , the right eigenvectors as {τ and the left eigenvectors as {τ . The latter are normalized such that Tr (τ is the Hilbert-Schmidt inner product of the space of operators acting on the system (see e.g. Ref. [56]). Using this notation and the master equation (2), the state of the system at time t is readily obtained 1 : where p defines the initial conditions. Eq. (5) represents the normal modes decomposition of the open dynamics of the system.

III. FIGURES OF MERIT FOR COLLECTIVE PHENOMENA
Next, we introduce the figures of merit of synchronization, subradiance, entanglement generation and correlations induced by the presence of a common bath during the relaxation dynamics of the two qubits, and comment on their experimental characterization. As anticipated, our aim is not to detect the simultaneous appearance of these four phenomena for some common initial conditions, but to characterize the capability of the common bath to induce some (in general different) open system dynamics displaying them. Therefore, the proposed figures of merit will not depend on time, but only on the parameters of the master equation (2). If we can establish that, for instance, the degree of synchronization witnesses that of bath-induced entanglement, then, in an experiment, the former (locally measured) will be instrumental for the choice of the proper conditions to optimize the entanglement production.
Note that the measures of entanglement and collectiveness will be taken as the supremum over time of some time-local indicator (namely, entanglement negativity and measure of correlations in the dynamics [11]). In contrast, the figures of merit for synchronization and subradiance will depend on the timescales during which these phenomena occur, and they will not be based on a time-local indicator. Indeed, the four measures aim to characterize the overall presence of each phenomenon for a given bath-driven quantum dynamics (i.e. quantum map), and not at a specific time of the evolution, as relevant for our purposes.

A. Quantum synchronization
The emergence of transient synchronized dynamics of local observables has been predicted in different quantum systems in recent years [47,[57][58][59][60][61][62][63][64][65][66][67][68][69][70]. This phenomenon is induced by some forms of dissipation into an environment (see for instance Refs. [6,71] for a review), and in some instance it can persist asymptotycally [58,64,70]. Other forms of synchronization in the quantum regime have also been explored, such as entrainment in presence of an external driving or synchronization among self-sustained oscillators [72][73][74]. Also the latter phenomenon is favored by collective dissipation as recently shown for optomechanical systems [63]. For the purpose of this work we focus on transient spontaneous synchronization mediated by the bath. In particular, we aim to investigate the synchronization of some peculiar oscillating coherences in the qubit dynamics, namely the ones captured by the mean values of the local observables σ x 1 (t) and σ x 2 (t). Without any external environment action, each qubit's coherences will oscillate as σ x k (t) = cos(ω k t + φ k ), where φ k is a fixed phase that depends on the initial conditions. Therefore, in the non-dissipative scenario each qubit oscillates at a different frequency. In contrast, the presence of dissipation can induce a regime where both qubits oscillate at a synchronized frequency [60]. In particular, under the dynamics driven by Eq. (2), the mean values of these local observables can be completely recovered by analyzing the Liouvillian block L 1 only [55,56]: with c Here, Re(λ (1) j ) and Im(λ (1) j ) denote the real and imaginary parts of the eigenvalue λ (1) j . Note that we have chosen to investigate synchronization of the mean values of σ x k . Any other choice of observable detecting oscillating coherences, such as σ y k , would be perfectly possible as well, and would lead to analogous results [56]. For the sake of simplicity, we focus on σ x k only. A well-known figure of merit to detect the synchronized dynamics of two observables is the Pearson coefficient [6], whose definition can be found in Eq. (C1) of Appendix C. Being a temporal correlation of the measured coherences dynamics of each qubit, this quantity can be easily experimentally measured. For our theoretical analysis in different parameter regimes we introduce a time-independent figure of merit. This will only depend on the coefficients of the master equation (2), having fixed some specific initial conditions. In fact, the emergence of synchronization is due to a mode of the block L 1 (let us say the one with eigenvalue λ (1) 4 ) decaying much slower than any other mode [57,58,60]. According to Eq. (6), at a certain time t S , σ x 1 (t) and σ x 2 (t) will synchronize at frequency Im(λ 4 ). A measure of synchronization can be introduced as follows: Syn describes the ratio between the decay time of the slowest-decaying mode λ 4 and the synchronization time t S , defined as the time at which the contribution to σ x 1 (t) and σ x 2 (t) of this mode is 100-times bigger than the one of any other mode (see the definition in Eq. (C2) of Appendix C). Choosing the number 100 in the definition of synchronization time (and correspondingly, in the factor log 100 in Eq. (7)) is heuristic. We do so in order to characterize the "disappearance" of all modes but the synchronized one, however, other choices are possible. Eventually, what we are characterizing is the trade-off between the dissipation timescale and synchronization timescale, and the qualitative behavior of the figure of merit Syn should not depend on this numerical choice. It is important to choose the same number in Eq. (7) and in Eq. (C2) (i.e., if we set log 50 in Eq. (7) then we will define the synchronization time as the time at which the synchronization mode is 50-times larger than any other one). By doing so, Syn > 1 indicates that synchronization appears before the "relaxation" of the slowest decaying mode. The bigger Syn, the more detectable the synchronization in the system. We choose ρ S (0) = |ψ Syn ψ Syn | , as initial state of the evolution, where |ψ Syn = (cos π/4 |e + sin π/4 |g ) ⊗ (cos π/3 |e + i sin π/3 |g ), as such a state has a large amount of initial coherence and at the same time is not symmetric under the exchange of the two qubits, which avoids pathological behaviors for some choice of parameters due to symmetry. In the case of transmon qubits, this separable state can be prepared by standard techniques using microwave Rabi pulses with frequencies ω 1 and ω 2 applied to to each qubit respectively.

B. Subradiance
In order to introduce a measure able to capture the appearance of subradiance in our model, we focus on the observables describing the population of the excited state of each qubit, defined as: The dynamics of P e 1 and P e 2 can be completely recovered by analyzing the Liouvillian block L 0 only [55,56]: where h As said in the introduction, subradiance entails the presence of a slowly-decaying collective mode in the dynamics of both the observables P e 1 (t) and P e 2 (t) defined in Eq. (8), which remains alive even when all the other modes have disappeared. L 0 always has at least one zero eigenvalue, corresponding to the steady state of the dynamics [55]. We do not consider the latter, and we focus on the remaining five. Analogously to the problem of quantum synchronization, we introduce the measure of subradiance as: where λ (0) 5 is the slowest-decaying mode of L 0 (apart from the non-decaying one) and t B is the subradiance time at which the slowest-decaying component in the mean value of P e 1 (t B ) and P e 2 (t B ) is 100-times bigger than that of any other mode (apart from the steady state one), whose expression can be found in Eq. (C4) of Appendix C. As for the case of synchronization, choosing the number 100 in the above definition is heuristic and other choices are possible (see the discussion in the previous section). We assume to start the evolution in the subradiant state (|eg − |ge )/ √ 2, which in a perfectly symmetric scenario (ω 1 = ω 2 , g 1 = g 2 and without local dissipation) lives in a decoherence-free subspace. This state is a Bell state that can be prepared by standard methods (see also Section 5). For example, starting with the qubits in the state |g |g , we can apply on the first qubit an X gate (σ x with our notations) followed by a Hadamard gate, then use this qubit as a control for a CNOT gate with the second qubit as target. Finally, we apply another X gate on the first qubit. As Sub describes the ratio between and the decay time of the slowest-decaying mode λ (0) 5 and the subradiance time, Sub > 1 implies that subradiance appears before the "relaxation" of the slowest decaying mode. If we consider the symmetrical scenario leading to a decoherence-free subspace, we have Sub = ∞, since the initial state would be a steady state of the evolution. In general terms, the bigger Sub the more detectable the subradiance in the system.
As a measure of entanglement of the state of the system we choose the negativity, which is a well-defined, easily computable entanglement monotone [100]. We expect any other measure of entanglement, such as entanglement of formation, to produce very similar qualitative results [101]. Negativity is given by where r j (t) are the eingevalues of the partial transpose of the state of the system ρ S (t) with respect to the second qubit. ρ S (t) evolves according to the master equation (2), therefore the negativity depends on time as well, and we choose as measure of the entangling power of the bath the maximum value over the whole evolution: This definition clearly depends on the initial conditions. As separable initial state of the evolution, we set ρ S (0) = ρ C = 1/4 j,k,l,m=e,g |jk lm|, which is the maximally coherent state. We choose so by considering that the majority of quantum algorithms start from ρ C , given that the latter is obtained by applying a Hadamard gate on the state with all the qubits in 0 (here |gg ) [102], which on a superconducting platform can be realized by applying resonant microwave Rabi pulses to each qubit. Some different choices of initial state will be discussed in Appendix D and Sec. IV B 3. The computation of N M requires the complete knowledge of the state of the system ρ S (t) at any time t. This can be obtained experimentally by means of the full-state tomography [95,103], based on joint measurements on both the qubits at the same time, which is highly resource-intensive.

D. Collectiveness of the dynamics
To quantify the collectiveness of the dynamics of two qubits driven by a common bath, we employ a measure recently introduced by Rivas and Müller [11]. It quantifies how much a given quantum evolution E(t) of two subparties 1 and 2 differs from a totally uncorrelated dynamics written as The measure is based on the Choi-Jamiołkowski isomorphism which makes use of the state |Ψ SS , living in H S ⊗ H S , where H S is the Hilbert space of the system and H S a copy of it. |Ψ SS is defined as: The subscripts indicate that the generic single-qubit state |j refers to the qubit 1 or 2 of H S or to the qubit 1 or 2 of H S . Let us name E(t) = exp(Lt) the quantum map describing the evolution until time t. Then, we can introduce a 16 Since E(t) is completely positive, Φ E (t) is positive-semidefinite and has trace 1, i.e. is a state in the Hilbert space H S ⊗ H S . The measure of correlations in the dynamicsĪ is defined as the normalized quantum mutual information of Φ E (t): where S is the von Neumann entropy and Tr 11 (22 ) is the partial trace on the qubit 1 and 1 (2 and 2 ). 0 ≤Ī(E(t)) ≤ 1, it is null if and only if the dynamics is uncorrelated, and it satisfies a fundamental principle which allows one to consider the correlations in the dynamics as a resource (see Ref. [11] for details), therefore it is a well-defined measure. It can be shown that it reaches the valueĪ( is a measure for the correlation in the quantum map that evolves the state until time t. Since we would like to have a figure of merit independent of time, following the original reference [11] we choose: For convenience, let us term the measure in Eq. (16) as collectiveness. The experimental computation ofĪ M requires quantum process tomography which is highly non-trivial [12], and in the present work we consider this figure of merit only as an abstract indicator to be compared with the other measures.

IV. RESULTS
Although synchronization, subradiance, entanglement generation and correlations in the dynamics of decoupled qubits share a common origin in this system, namely the presence of a collective bath, they are essentially four different physical phenomena. From a mathematical perspective, this is displayed by the block structure of the Liouvillian superoperator discussed in Appendix B. We observe that quantum synchronization is described by the block L 1 (Eq. (6)), subradiance by the block L 0 (Eq. (9)), while both the negativity Eq. (11) and the collectiveness Eq. (15) may in general depend on all the blocks of the Liouvillian. We know that, only in the case with zero temperature, a parametric dependence induces the same separation between the decay rates characterizing synchronization and subradiance [56], but in more general scenarios there is no clear relation between the phenomena, given that each Liouvillian block has an independent behavior. Therefore, we first discuss the behaviour of the four figures of merit for some simple limits of the model parameters, which can be computed analytically, and then we perform a more general, deep numerical investigation in order to compare the strength of each phenomenon in different scenarios and their possible concomitance. Note that the discussion on the analytical limits is valid not only for the case of transmon qubits immersed in an Ohmic thermal bath, but for any system of two qubits dissipatively coupled to an environment.
A. Analytical limits In the zero temperature case, the form of the eigenvalues of each block of the Liouvillian superoperator is known, as discussed in Appendix E. In particular, we find that the separation between the eigenvalues corresponding to the two slowest modes of both L 0 and L 1 has the same expression, and can be written as (the principal square root is taken): with ∆ω = ω 1 − ω 2 . Re(∆λ) provides the separation between the intensity of the decay rates of the slowest modes, that plays an important role in the definition of Syn and Sub respectively in Eqs. (7) and (10). The case of high detuning is trivial, since it leads to the validity of the full secular approximation [54], i.e. the qubit-qubit correlations in the master equation become negligible. As for the time-scales separation, it does not arise since ∆λ is almost purely imaginary. More interesting for our purpose is the case of small detuning. We start assuming balanced weights, g 1 = g 2 . Under these Looking at ∆λ as a function of the detuning, we easily observe that d d∆ω Re(∆λ) < 0, i.e. the separation between the real part of the eigenvalues of the two slowest eigenmodes decreases as a function of the detuning. According to the discussion in Appendix E, for g 1 = g 2 and ω 1 = ω 2 the eigenvectors associated to the slowest eigenvalues have the maximum projection over σ x 1 and σ x 2 for synchronization, and P e 1 and P e 2 for subradiance. Then, we can conclude that, in the balanced case g 1 = g 2 and for T = 0, increasing the detuning hinders Syn and Sub. Moreover, for g 1 = g 2 and finite detuning, decreasing the strength of the qubit-bath coupling also hinders Syn and Sub, since lim µ→0 + ∆λ ≈ √ −∆ω 2 , which has zero real part. As for entanglement and collectiveness, these phenomena depend on the amount of correlations that dynamically build up among the qubits. Therefore, they will not arise in the limits in which such correlations disappear. In particular, N M = 0 and I M = 0 (i) for µ 2 1/T 1 (the collective action of the bath is negligible with respect to the local incoherent dissipation), (ii) for g 2 → 2, g 1 → 0 and viceversa (the bath acts only on a single qubit and not on the other), and (iii) for ∆ω = O(ω 1 ) (the detuning is of the order of the qubit frequencies, and therefore the full secular approximation applies). We observe that entanglement and collectiveness display the same behavior as synchronization and subradiance in these limits.

T → ∞
In the infinite temperature limit, exact analytical results can be derived for g 1 = g 2 , ∆ω = 0, so that there is a single dissipative coefficient γ = γ ↓ jk = γ ↑ jk . Under these assumptions, first of all we observe that the measure of subradiance has an infinite value, Sub = ∞, since the subradiant state lives in a decoherence-free subspace, i.e. λ , and therefore, in general terms we bring forward the subradiance time. For continuity, we expect a similar behavior also when the detuning is not zero anymore (and therefore the subradiance measure has a finite value), that is to say, subradiance is improved for higher temperatures.
In the case of synchronization, we can compute the eigenvalues of the block L 1 , whose expressions can be found in Appendix E. We observe that the difference between the two slowest eigenvalues is purely imaginary: with s + = s ↓ 12 + s ↑ 21 . Therefore, the real part of the difference is zero and synchronization can never appear, since there will always be at least two modes with different frequency in the dynamics of the observables in Eq. (6). Hence, Syn → 0 for T → ∞.
Finally, let us address entanglement and collectiveness. For these figures of merit we are not able to derive an analytical solution at any time. However, previous works have shown that the entanglement generated by the action of a common bath vanishes at infinite temperature [83,104]. Let us understand why. First, note that as we increase the temperature, we strongly increase the dissipative coefficient γ, proportional to coth(β ω/2), while we do not observe such a strong enhancement of the collective Lamb-shift term s + (see e.g. Eq. (A1) in Appendix A). As extensively discussed in the literature, the capability of a common bath to entangle a couple of non-interacting qubits strongly depends on the collective Lamb-shift term [105]. This is also suggested by our analysis of the sufficient conditions for an entangling bath presented in Appendix D: if s + → 0 and T → ∞, the sufficient condition is not satisfied anymore. Therefore, we can expect that, in the limit γ s + , this inequality for the sufficient condition in Appendix D is only weakly satisfied, an indication of small entanglement. If this is true, we obtain that N M → 0 + for T → ∞. In Appendix E we show that this hypothesis is correct through a careful study of the trade-off between γ and s + . Furthermore, we can employ a similar argument for the measure of collectiveness, showing that it decreases (although does not vanish) for T → ∞: if γ s + , the collective Lamb-shift term becomes negligible in the master equation and the collective action of the quantum map derives only from the dissipative part. On the contrary, when the Lamb-shift term is relevant, its action brings an additional collective term to the master equation, increasing the value ofĪ M . The analysis in Appendix E confirms this claim.

B. Numerical comparison in different scenarios
In the following, we will explore several situations by studying how the figures of merit vary for different parameters in the two qubits platform introduced in Sec. II. We will analyze their behavior against the detuning, the bath temperature, the systembath coupling, and the unbalancing of the coupling, also exploring the dependence on initial conditions. For the sake of clarity, we will first discuss the case g 1 = g 2 and then move to the unbalanced case, corresponding to a situation in which collective dissipation acts with different strength on the two qubits. The results are displayed in Fig. 1 for the balanced case g 1 = g 2 , where the four indicators are plotted against the detuning (1(a-d)) and against the temperature (1(e-h)) and in Fig. 2 for the unbalanced case g 1 = g 2 . Further results are presented in Fig. 3. In all the scenarios, anticipating the relevant case for the platform we are going to introduce in Sec. V, we set an Ohmic spectral density of the bath (see Eq. (A4)).

Balanced couplings
Let us start analyzing the role played by the system-bath coupling: the higher the coupling strength, the stronger the effects of the common bath, and therefore the larger the four measures. This can be observed, for instance, in Figs. 1(a)-(d) and in Figs. 3(e)-(h). If µ is too weak, the unitary evolution driven by the system Hamiltonian will play the only relevant role in the dynamics, and no collectiveness, entanglement, synchronization or subradiance will appear. For instance, with the parameters employed in Figs. 1(a)-(d) (i.e. g 1 = g 2 = 1, β = 10/ ω 1 , T 1 = 3 × 10 5 /ω 1 ), we observe that for µ 10 −2.5 ω 1 and ∆ω ≈ 0.01ω 1 neither synchronization nor subradiance emerge before thermalization, while both negativity and collectiveness are negligible.
As for the detuning, in all scenarios, N M andĪ M decrease as soon as ∆ω increases, consistently with a vanishing dissipative coupling. Indeed, in spite of the presence of a common bath, in this regime the dynamics can be well approximated by a full We observe that synchronization shows a slight asymmetry between the scenarios with ∆ω > 0 and ∆ω < 0. This is due to the fact that Syn is particularly sensible to the change of the absolute value of the frequencies ω 1 and ω 2 , and not only of the detuning ∆ω.
More complex is the dependence with the temperature, as depicted in Figs. 1(e)-(h) and Figs. 3(e)-(h). In this case, the measure of subradiance displays notable differences with respect to the other ones: it is greatly enhanced by higher temperatures, as found analytically, unlike the other phenomena. We indeed observe that, even for the case with finite detuning, the slowest decay rate is not significantly affected by the value of the temperature, while all the remaining ones are speeded up when temperature increases, and therefore Sub is enhanced for higher temperatures. As for synchronization, in agreement with the analytical limit and what was partially discussed in the literature [106,107], we observe that while a higher temperature speeds up the relevant decay rates (formally the eigenvalues of the block L 1 , see Eq. (6)), it does not increase the gap between the two slowest eigenvalues. As a consequence, the synchronization time does not change considerably, while the relaxation time occurs way before than at lower temperatures. The lack of a significant transient where synchronization can be observed is indeed captured by a decrease of the measure of Eq. (7). Furthermore, entanglement smoothly disappears as temperature increases, vanishing for β ω 1 1, in agreement with the analytical limit of Sec. IV A. Finally, we observe that the collectiveness decreases as well when the temperature increases, although, contrary to the negativity and synchronization, its value never vanishes for T → ∞. Indeed, even at infinite temperature the master equation describing the dynamics driven by a common bath is different from the one driven by two independent local baths.

Unbalanced couplings
The maximal collective effect is achieved when the coupling between the qubits and the bath is balanced, i.e. g 1 = g 2 . A less relevant scenario appears when this is not true anymore, and the collectiveness of the evolution inevitably decreases. Let us analyze this case: the bigger the imbalance between the weights of the dissipative coupling g 1 and g 2 , the smaller the collective effects in the dynamics. Indeed, both collectiveness and negativity decrease as the imbalance increases (see Figs. 3(k)-(l) and 2(c)-(d)). If we weaken the coupling of one of the two qubits by varying g 1 and g 2 , we monotonically decrease the values of N M andĪ M , reaching zero in the trivial limit g 1 → 0, g 2 → 2. On the other hand, a significative imbalance between the weights of the dissipative interaction g 1 and g 2 breaks the possibility of using synchronization and subradiance as signatures of entanglement generation. Indeed, the appearance of synchronization and subradiance depends not only on the collectiveness of the evolution, but also on the gap between the eigenvalue of the slowest eigenmode and all the remaining ones. By unbalancing the coupling, we may enhance this gap, and therefore increase the values of Syn and Sub (as for instance already discussed in Ref. [45]). This behavior is observed in Figs. 2 and 3(i)-(j), which shows that, by unbalancing g 1 and g 2 , the separation between the eigenvalues of both synchronization and subradiance is enhanced in a stronger way than the collectiveness of the relevant eigenmodes is decreased, and therefore we obtain higher values of Syn and Sub. Finally, let us comment that, in the unbalanced scenario, the behavior of synchronization and subradiance when varying the temperature is now completely different from the balanced case (see Figs. 3(i)-(j)). Indeed, the arguments we used in the previous section to understand the behavior of the figures of merit as a function of the temperature are now not valid anymore, and, for example, for some values of g 1 and g 2 we obtain a higher synchronization for higher temperatures.

Dependence on the initial conditions
Let us finally try to understand whether our results depend on the specific initial states we have chosen. In order to do so, for the case of synchronization and subradiance, we have computed a modified form of their figure of merit in which we do not consider the initial conditions. In particular, we have set each p Focusing on the entanglement generation, we can explore different initial conditions by considering the states discussed in Appendix D. We have found that the sufficient conditions for having entanglement at time t → 0 + starting from a given initial FIG. 4. Circuit diagram representing two transmon qubits, characterized by the capacitances C1 and C2, capacitively coupled to the same resistor R, respectively through the capacitors CL and CR. The resistor plays the role of a common bath with Ohmic spectral density, that can be heated up by an external bias current I b .
state, actually can be associated to the behavior of the entanglement generation throughout the whole evolution. For instance, we observe that no entanglement is generated by the bath when starting the dynamics in the pure states |ee or |gg , which as proven in Appendix D can never display entanglement at t → 0 + . We have then considered the initial states |eg and |ge . We observe that, in both cases, the maximum value of negativity generated during the evolution is higher than the one observed in Figs. 1, 2 and 3, where the choice was ρ S (0) = ρ C . Curiously, the sufficient conditions of Appendix D also suggest that these states have, in general, a stronger power of generating entanglement at t → 0 + . Moreover, we observe that the asymmetry between qubit 1 and qubit 2 is reflected in the behavior of N M as a function of the detuning: for instance, if we start the dynamics in |eg we find stronger entanglement for ∆ω = +x than for ∆ω = −x, and viceversa. Still, the dependence as a function of ∆ω does not change, and all the results of the previous discussion hold. In all the other scenarios the behavior of N M is analogous to the one depicted in the figures of the main text. We therefore conclude that the dependence on the initial conditions does not affect the findings of our work, and the same conclusions can be drawn for general scenarios.

C. Extensions to multi-qubit systems
In Sec. IV A we have discussed how the four figures of merit behave in some particular limits, which are valid independently of the chosen model or physical scenario, and hold for any system of two qubits dissipatively coupled to a thermal bath. In Sec. IV B we have performed an extensive numerical investigation for different scenario inspired by two transmon qubits immersed in an Ohmic thermal bath. The extensions of these results to the multi-qubit case would deserve a separate study, but we can anticipate some issues in this more complex scenario. Adding a single qubit, one may take inspiration from the discussion in Ref. [108], focused on three harmonic oscillators. Here, the different strengths of each coupling interaction with the bath play a relevant role, and different conditions for decoherence-free subspaces may appear in some parameter manyfolds. Furthermore, it is crucial to distinguish between synchronization of only two of the three qubits immersed in a common bath and synchronization of the whole system. Correspondingly, one may decide to compare it with multipartite entanglement [109] of the system, or just bipartite entanglement in a qubit pair (similar issues have been investigated in Ref. [110]). A recent work also studies different synchronization features in a three-qubit system [107]. The case of multiple qubits is even more complex, and one may also address the coexistence of common and local baths acting on a subsystem of the qubits only, as studied for harmonic oscillators in Refs. [58,64].

V. AN EXPERIMENTAL PROPOSAL TO TEST OUR FINDINGS
We propose here an experimental platform that would allow one to test our predictions in a controllable environment and beyond mathematical approximations, eventually also allowing for the exploration of further regimes. The phenomena we have discussed can be observed by using two superconducting qubits, in a simple experiment as depicted in Fig. 4: two transmon qubits of frequency ω 1 and ω 2 are coupled to the same resistor through the capacitors C L and C R . Let us assume that the qubits have total capacitances (shunt plus intrinsic junction capacitance) C 1 and respectively C 2 . In the weak coupling limit C L , C R C 1 , C 2 , using standard quantum network analysis in the presence of dissipation [111], we obtain the Hamiltonian Eq. (1). Here the resistor plays the role of a thermal bath with an Ohmic spectral density proportional to the resistance R. The parameters µ, g 1 and g 2 can be tuned by varying the values of the capacitors in the circuit. Indeed, we obtain: and µ is a constant with the units of energy, proportional to the resistance and to (C L + C R )/(C 1 + C 2 ), that with the above assumption satisfies the weak coupling limit. The temperature of the thermal bath is the effective temperature at which the resistor is dissipating energy, typically some tens of millikelvin. Let us consider a specific set of parameters: according to the discussion in the previous section, synchronization and subradiance are good signatures of entanglement generation in the balanced case, i.e. g 1 = g 2 . To obtain so, we choose C L and C R such that C L /(C L + C 1 ) = C R /(C R + C 2 ). As reasonable qubit frequencies, we choose ω 1 = 2π × 5 GHz and ω 2 = 2π × 4.95 GHz, so ∆ω = 0.01ω 1 = 2π × 0.05 GHz. In this configuration, state preparation can be achieved by standard methods in circuit quantum electrodynamics. Single qubit operations can be realized by applying microwaves tones in resonance with the corresponding qubit frequency. Since the coupling is fixed by the choice of the capacitors, two-qubit operations such as the CNOT needed to prepare the subradiant state can be realized naturally in this configuration by cross-resonance. In this method, two microwave pulses in resonance with the other qubit's frequency are applied simultaneously, which results in an effective tunable qubit-qubit coupling [112][113][114]. Note that in order to achieve high-fidelity state preparation, additional flux bias can be added to the transmons, enabling fast tunability of the frequencies. By suitably tuning the magnitude of the coupling capacitances and resistance, we set the coupling constant as µ = 10 −1.5 ω 1 ≈ h × 0.16 GHz. Taking into account the values that are usually measured in real experiments with transmon qubits [23], we choose the local relaxation time T 1 = 3 × 10 5 /ω 1 ≈ 10 µs. Finally, we set the inverse temperature β such that ω 1 β = 10, i.e. β ≈ 3 × 10 24 J −1 . This is reasonable, considering that it corresponds to a temperature T ≈ 24 mK. This choice allows us to neglect the leakages to the third level of the transmon qubit, that would be proportional, according to Maxwell-Boltzmann distribution, to a very small factor of the order of ≈ 1/e 20 . Higher temperatures (as in some parameter regions of Figs. 1(e)-(h) and 3(e)-(l)), however, may induce some non-negligible leakages to higher energy levels. In this case, cavity filters at the input of the two qubits can be installed to reduce the effect of spurious excitations to frequencies other than the qubit frequencies [26] -see for instance the scheme in Ref. [17], where a couple of resonant LC circuits are employed to suitably filter the radiation coming from separate thermal baths.
We can now investigate the open dynamics of the transmon qubits under the action of the resistor, with the parameters set as above, for the two different initial states associated to the figures of merit of synchronization and subradiance, so as to be able to track the evolution of respectively σ x 1,2 (t) and P e 1,2 (t). The complete information about these observables is acquired by local measurements on each qubit using the standard dispersive scheme. In this scheme, P e 1,2 (t) are obtained directly as populations of the excited state, while for measuring σ x 1,2 (t) we need to apply an X-pulse before the measurement. After performing these measurements, we can compute the synchronization and subradiance measures Syn and Sub as discussed in Sec. III. The results of our simulation are plotted in Fig. 5, where we also compare them with the negativity and collectiveness as a function of time (subfig. (c)). The Pearson coefficient depicted in Fig. 5(a) shows the appearance of an antisynchronized regime (upper inset) after an incoherent transient (lower inset). Still, it does not correctly detect the synchronization time t S defined in Eq. (C2), since it provides a value higher than 99% for t ≈ 4.2 × 10 2 /ω 1 . On the contrary, a more accurate spectral analysis shows that the time at which the slowest decaying mode in the dynamics of σ x 1,2 (t) is 100−times bigger than all the other ones is t S ≈ 8.3 × 10 2 /ω 1 , highlighting how the Pearson coefficient alone is not enough accurate to extract the value of Syn. Fig. 5(b) depicts the decay of P e 1,2 (t) as a function of time. We observe the emergence of a subradiant mode which remains alive also at late times, and a spectral analysis reveals the subradiance time Eq. (C4) t B ≈ 3.4 × 10 2 /ω 1 . From this analysis, we extract the value of the figures of merit Syn ≈ 428, Sub ≈ 529, according to Eqs. (7) and (10). Finally, let us focus on the negativity and collectiveness in Fig. 5(c): N (t) increases at t → 0 + , showing that entanglement is generated as soon as the dynamics begins, as predicted in Appendix D. After reaching the maximum value at t ≈ 21/ω 1 , it decreases towards zero. Then, a less intense re-birth of entanglement is observed around t ≈ 50/ω 1 , as already predicted in similar models [54]. The collectiveness displays a similar behavior, reaching the top at later times t ≈ 25/ω 1 , and then decreasing with an oscillating behavior. Consistently with the original paper [11], it decays toward 0 at infinite time. From the evolution we obtain the values N M ≈ 0.37 andĪ M ≈ 0.81. It is important to note that all these characteristic times are smaller than the relaxation time T 1 ; this gives a realistic temporal window for performing measurements and observing these effects.

VI. CONCLUSIONS
In this work we have investigated the behavior of quantum synchronization, subradiance, entanglement, and collectiveness of the dynamics, in a broad scenario of two superconducting transmon qubits dissipating into a common thermal bath. Although the four phenomena are essentially different, and a general one-to-one correspondence between them cannot be established, they share a common cause, namely the action of the collective bath. Inspired by this, we have addressed in which scenarios synchronization and subradiance (equipped with well-defined, time-independent quantifiers) can act as signatures of entanglement generation. We also have compared synchronization and subradiance with a measure of the collectiveness of the evolution, i.e. of how much the dynamics differs from one described by two separable quantum maps acting locally on each qubit. Finally, we have provided an experimental proposal for a system of two transmon qubits coupled to the same resistor, which mimics a thermal bath.
We have found that quantum synchronization acts as a good signature of entanglement generation, i.e. the behavior of its measure varies in the same way of the measure of entanglement production as a function of the model parameters, except in the scenario with the qubits coupled to the bath in an unbalanced way, which favors the emergence of a slowly-decaying eigenmode, while it decreases the collectiveness of the dynamics. The measure of subradiance follows the same behavior, but it fails to describe entanglement generation for high temperatures, where the slowly-decaying subradiant eigenmode is not affected by temperature, while the decay rate of all the other modes are enhanced. It is interesting to notice that, in the scenarios where synchronization and subradiance cannot be considered as reliable signatures of entanglement generation, the common bath can easily induce a certain collective phenomenon while struggles to produce a different one. Even more curiously, there are cases for which we may experimental detect intense subradiance and quantum synchronization (e.g. for very strongly unbalanced weights, see Figs. 2 and 3(i)-(l)), which are clearly collective phenomena, while the measure of collectiveness is very weak.
Our results also shed new light on the relation between entanglement generation and collectiveness. We have observed that, at least if the initial state of the system is symmetric with respect to the switching of the two qubits, their values follow exactly the same behavior, with the only exception that the collectiveness decreases but does not reach a null value for T → ∞, while entanglement does. Our findings seem to suggest that the two measures are closely related, and a tighter connection between them may be found. Finally, our analysis establishes quantitatively the correlated dissipation strength needed for a pair of transmon qubits to display effects such as entanglement and spontaneous synchronization.
Remarkably, the advantage of employing synchronization or subradiance as signatures of entanglement generation consists in the fact that their measures can be computed by monitoring the mean value as a function of time of local observables only, without relying on joint measurements, which are on the contrary necessary to calculate any measure of entanglement, such as the negativity. We expect this to be especially useful when the system tomography is not available. In particular, it would be interesting to extend our results to systems of multiple qubits that realize a quantum computing platform, for which tomography is not feasible anymore. Furthermore, the experimental proposal we discuss may be implemented with the currently available technology and know-how about superconducting qubits. Its realization would be the first demonstration of spontaneous quantum synchronization, which does not require the presence of an external driving field, in the deep quantum regime.

Acknowledgements
The numerical analyses have been performed with the assistance of QuTiP, the Quantum Toolbox in Python [115]. This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agree- We assume that the bath is in a thermal state at temperature T . The coefficients of the master equation (2) read: where Γ β (ω) is the one-side Fourier transform of the autocorrelation functions of the bath, which depends on the spectral density Eq. (A3) and on the temperature of the thermal bath through β = 1/k B T [1,54]: where N β (ω) = 1 e β ω −1 = 1 2 coth β ω 2 − 1 , and J(ω) is the spectral density of the bath, defined as: where f k are introduced in Eq. (1). In the analysis of Sec. IV we have employed an Ohmic spectral density, i.e.: where ω C is a cut-off frequency that we have set as ω C = 20ω 1 , and we have chosen to renormalize the spectral density by the frequency of the first qubit.
In the decay coefficients of Eq. (A1) we have added a local dissipative decay rate 1/T 1 on each qubit, which represents the action of a phenomenological local bath at zero temperature. Its effects are therefore not relevant for the absorption coefficients, and we have also neglected its contribution to the Lamb-shift. These choices are motivated by the experimental decay observed in a single superconducting qubit, which reaches to a good approximation the ground state at infinite time, corresponding to a thermal bath at zero temperature. On the contrary, we allow for the common bath to have a finite temperature, which experimentally can be realized by local heating of the resistor connecting the qubits.

Appendix B: Blocks of the Liouvillian superoperator
The master equation (2) is in partial secular approximation and therefore a fundamental symmetry on the superoperator level emerges [55]: [L, N ] = 0, where N = [P e 1 + P e 2 , · ] is the number superoperator (P e 1 and P e 2 are defined in Eq. (8)). For this reason, the Liouvillian superoperator can be block-diagonalized with each block labeled by a different eigenvalue of N . Let us work in the basis of the space of operators {|jk lm|} j,k,l,m=e,g . N is diagonal in this basis, and its eigenvalues correspond to the number of excited qubit states in the ket minus the number of excited qubit states in the bra. Therefore, a single eigenvalue d can assume the values d = −2, −1, 0, 1, 2, and the Liouvillian superoperator is divided into five independent blocks labeled by d: L = 2 d=−2 L d . In particular, the block L 0 will act on the basis vectors |ee ee|, |eg eg| , |eg ge| , |ge eg| , |eg eg|, |gg gg|, L 1 on |ee eg| , |ee ge| , |eg gg| , |ge gg|, L −1 on |eg ee| , |ge ee| , |gg eg| , |gg ge|, L 2 on |ee gg| and L −2 on |gg ee|. In these bases, we have L −1 = L * 1 and L −2 = L * 2 (more details in Refs. [55,56]). The spectral analysis of each block of the Liouvillian superoperator leads to the solution of the dynamics as given by Eq. (5).
Appendix C: Derivation of the figures of merit

Quantum synchronization
The Pearson coefficient is defined in a temporal window ∆t as: In a scenario in which the dynamics of the observables reaches the perfect in-phase synchronization at a given time t S , we observe that the value of the Pearson coefficient stabilizes to 1 at the same time, for a suitable time window ∆t [57,60,107]. A slightly different definition of the Pearson coefficient allows us to catch the emergence of synchronization with a given phase shift [6].
Let us now focus on a time-independent measure of quantum synchronization. It will depend only on the separation between the real part of the eigenvalues of the modes of L 1 , on each corresponding inner product c (1) jk and on the initial conditions, i.e. on the state at which we initialize the system at t = 0. We start the dynamics from the state ρ S (0) = |ψ Syn ψ Syn | defined in Sec. III A, and we monitor the observables σ x 1 (t) and σ x 2 (t). The free evolution of each qubit would give σ x k (t) = cos(ω k t + φ k ) where φ k is a phase depending on the initial conditions, while the presence of a common bath can lead to the synchronization of their dynamics, i.e. to a situation in which they oscillate at the same frequency. To do so, their evolution must be monochromatic, that is to say, the main contribution to the mean value of σ x k (t) in Eq. (6) must be given by a single eigenmode of the Liouvillian block L 1 only. This happens when all the other modes have disappeared due to their faster decay. Let us say that synchronization emerges at the time t S , when the contribution of the slowest-decay mode with eigenvalue λ (1) 4 to σ x 1 (t) and σ x 2 (t) is 100-times bigger than the one of any other mode. Intuitively, at t S the Pearson coefficient will have reached a value C ∆t (t S ) > 0.99. Using the notation introduced in Eq. (6), the synchronization time t S is estimated by inverting and then maximazing over j = 2, 3, 4 and k = 1, 2. Finally, we get: We observe that the synchronization time depends on the separation between the slowest decay rate and all the other ones, and on the inner product between the eigenvector of the slowest-decaying mode and the initial state, σ x 1 (t) and σ x 2 (t). If the two slowest decay rates do not coincide, there will always be a finite synchronization time at which the mean values of the qubit observables will oscillate at the same frequency. However, decoherence appears along the evolution and the synchronization time may be way longer than the time at which the system has lost all of its coherences, and the detection of synchronization would be almost impossible. To take into account both phenomena, we define the figure of merit of synchronization as the ratio between the synchronization time and the time at which the slowest-decaying mode is reduced to 1/100 of its initial value, and we obtain Eq. (7): .

(C3)
Experimentally, Syn can be estimated by monitoring the mean value of σ x 1 (t) and σ x 2 (t) until the time at which the spectral density of the signal reveals the presence of a single mode only, with an error of the 1%. Then, the decay rate of the remaining mode can be found through an exponential fit of the carrier wave. The higher Syn, the easier will be to detect quantum synchronization, since the signals of the mean value of σ x 1 (t) and σ x 2 (t) will be more intense at the synchronization time. If Syn = 1, then the synchronization appears exactly when the only remaining mode is reduced to 1/100 of its initial value.

Subradiance
The derivation of the figure of merit for subradiance is analogous to the one of quantum synchronization: we want to identify a collective slowly-decaying mode of the system decay which remains alive when all the other modes have vanished. The only differences consist in the observables we want to monitor, that now are P e 1 (t) and P e 2 (t) instead of σ x 1 (t) and σ x 2 (t), and in the initial conditions, that read ρ S (0) = (|eg −|ge )( eg|− ge|) 2 .
In a scenario without a collective bath acting on both the qubits, each of them would decay with an independent decay rate. On the contrary, the presence of a common bath creates a slowly-decaying collective mode whose component in both P e 1 (t) and P e 2 (t) survives after all the other ones have vanished, apart from the mode pertaining to the steady state which does not decay at all (say the one with eigenvalue λ (0) 6 = 0). We want to identify the long-surviving component, let us say corresponding to the eigenvalue λ (0) 5 , and the time at which it emerges. Following the discussion about quantum synchronization, we choose the subradiance time t B as the time at which the slowest-decaying component in the mean value of P e 1 (t B ) and P e 2 (t B ) is 100-times bigger than that of any other mode (apart from the steady state one). Using Eq.(9) and analogously to Eq.(C2), we have: As for the case of quantum synchronization, we define the figure of merit of subradiance as the ratio between the subradiance time t B and the time at which the component of the slowest-decaying mode λ (0) 5 is reduced to 1/100 of its initial value, and we obtain Eq. (10).
Analogously to the estimation of the figure of merit of quantum synchronization, in order to compute the value of the subradiance measure Eq. (10) in a real experiment one has to track the signals of the mean values P e 1 (t) and P e 2 (t) . In this way, we are able to detect the time t B at which one single decaying component of the signals is 100-times bigger than any other one (after having removed the steady-state one), and to estimate the corresponding decay rate.
• If ρ S (0) = |ge ge|, the sufficient condition for generating entanglement at time t → 0 + reads γ ↓ 22 γ ↑ 11 < γ ↓ 12 +γ ↑ 21 2 2 + s 2 + . In particular, we see that the bath always generates entanglement if T = 0. In any case, we observe that the quantity s 2 + is relevant to provide a general sufficient condition for the bath to be entangling. An analogous result is found for ρ S (0) = |eg eg|.
• If ρ S (0) = ρ C = 1/4 j,k,l,m=e,g |jk lm|, i.e. the state we choose as initial state of the evolution in our analysis, as discussed in Sec. III C, the condition is Once again, we recognize the important role the Lamb-shift terms s ± play in assuring a sufficient condition to generate entanglement.
Note that, a priori, if for an initial state we find that entanglement is not generated at time t → 0 + , this does not mean that entanglement will never appear during its evolution. Indeed, it may be the case that, after a certain time t * , the state is represented by a density matrix ρ for which the condition M 22 M 33 < |M 23 | 2 is verified. In other words, the method we have followed [99] is useful to provide sufficient conditions, while in order to find a necessary one we would need to check that M 22 M 33 > |M 23 | In the zero temperature case, we can derive analytically the eigenvalues and some of the eigenvectors of the Liouvillian superoperator L, as already discussed in Ref. [56]. In particular, the eigenvalues of the block L 0 read: where taking the principal square root. Note that we have sorted the eigenvalues in descending order of absolute value of their real part. The subradiant mode we are interested in has eigenvalue λ (0) 5 . The associated right eigenvector is τ .
Note that, according to Eq. (C2), we need to consider the difference between the real parts of λ 5 and the rest of eigenvalues of L 0 , excluding λ (0) 6 because it indicates the steady state eigenmode. Therefore, the smallest separation between the eigenvalues is given by Re(λ Let us now focus on the block L 1 . Its eigenvalues are: where once again we have sorted them in desceding order of absolute value of the real part. The right eigenvector associated to the slowest eigenvalue λ .
Once again, we find Re(λ 3 ) = Re(V ). Let us now calculate this quantity: As discussed in the main text, by assuming g 1 = g 2 and consequently γ ↓ 11 − γ ↓ 22 ≈ 0, clearly we have d d∆ω Re(V ) < 0. (1) 4 = 1 2 . This means that for the unbalanced scenario with g 1 = g 2 , the zero detuning configuration provides maximal symmetric weights for the definition of subradiance (h (0) 5k ) and synchronization (c (1) 4k ), therefore the smaller the detuning the higher Syn and Sub.

T → ∞
Let us consider the limit of T → ∞, and for simplicity g 1 = g 2 , ∆ω = 0. In this simplified scenario, the master equation depends on only three parameters: ω = ω 1 + s 11 = ω 2 + s 22 , γ = γ ↓ jk = γ ↑ jk , s + = s − = s ↓ 12 + s ↑ 21 , where s jj = s ↓ jj − s ↑ jj is a local Lamb-shift correction. Clearly, the measure of subradiance Sub always takes an infinite value, since, as discussed in the main text, the symmetric scenario makes a decoherence-free subspace emerge, associated to the subradiant mode. Let us now focus on synchronization. The eigenvalues of the block L 1 are the following: We observe that the two slowest eigenvalues have the same real part, while different imaginary part. This means that synchronization can never appear (Syn = 0), and after a transient time we will observe the presence of only two coexisting modes (corresponding to λ 3 and λ 4 ) oscillating with different frequencies. Therefore, we observe that the limit T → ∞ in the balanced scenario increases the value of the measure of subradiance Sub, while it decreases the value of Syn, not allowing synchronization to emerge. This result is not valid anymore in the presence of unbalanced weights, i.e. g 1 = g 2 .
Let us now consider entanglement and collectiveness in the same scenario. Unfortunately, a readable analytical solution for the negativity N (t) in Eq. (11) and collectiveness of the dynamics up to time t,Ī(E(t)) in Eq. (15), is not available. However, we can understand their behavior by looking at the ratio between the Lamb-shift collective term s + and the dissipative collective term γ. To do so, we analyze the behavior of the following master equation: with Liouvillian corresponding to a scenario with a thermal bath at infinite temperature, but keeping the possibility to freely set the parameters ω, γ and s + , which now do not depend on the spectral density of a specific bath. To work with dimensionsless units, we set = 1 and ω = 1, so that we only need to focus on γ and s + . In Fig. 6 we depict the entanglement and collectiveness as a function of γ and s + . In Figs. 6(a),(c) we fix s + and make the time vary: we observe that, for fixed s + , higher values of γ lead to lower values of both entanglement and collectiveness. In Figs. 6(b),(d) we fix a early time of the evolution, so as to compute the figures of merit in the first rising oscillation (compare with the simulated curve in Fig. 5(c)), and vary γ and s + . We clearly identify a trade-off between these two parameters: the