# Unraveling Dicke Superradiant Decay with Separable Coherent Spin States

Pedro Rosario, Luiz O. R. Solak, André Cidrim, Romain Bachelard, Johannes Schachenmayer

► To cite this version:

Pedro Rosario, Luiz O. R. Solak, André Cidrim, Romain Bachelard, Johannes Schachenmayer. Unraveling Dicke Superradiant Decay with Separable Coherent Spin States. Physical Review Letters, 2025, 135 (13), pp.133602. 10.1103/xcxr-sm9c. hal-05367170

HAL Id: hal-05367170

https://cnrs.hal.science/hal-05367170v1

Submitted on 15 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Unraveling Dicke Superradiant Decay with Separable Coherent Spin States

P. Rosario,¹,² L. O. R. Solak,¹,² A. Cidrim,¹ R. Bachelard,¹ and J. Schachenmayer²

¹Departamento de Física, Universidade Federal de São Carlos,
Rodovia Washington Luís, km 235—SP-310, 13565-905 São Carlos, SP, Brazil

²CESQ/ISIS (UMR 7006), CNRS and Université de Strasbourg, 67000 Strasbourg, France

(Dated: September 3, 2025)

We show that idealized Dicke superradiant decay from a fully inverted state can at all times be described by a positive statistical mixture of coherent spin states (CSS). Since CSS are separable, this implies that no entanglement is involved in Dicke decay. Based on this result, we introduce a new numerical quantum trajectory approach leading to low-entanglement unravelings. This opens up new possibilities for employing matrix product state (MPS) techniques for large-scale numerical simulations with collective decay processes.

**Introduction**—The Dicke model [1], a cornerstone of quantum optics, describes the dynamics of a large, all-to-all-coupled ensemble of N two-level atoms resulting from their coupling to the same single mode of the electromagnetic field. The model predicts a phenomenon coined as pure superradiance: Atoms spontaneously emit light cooperatively, leading to an intense, short burst of radiation whose intensity scales with $N^2$ [2]. Although semiclassical approaches have successfully described aspects of superradiance (e.g., burst times $t_B \propto \ln N$) [3, 4], hinting at the classicality of the phenomenon, understanding whether quantum correlations are generated throughout Dicke dynamics has long remained a challenge [5–9]. Today, a common picture is to consider the many-body state evolving in time by following a cascade over states symmetric under atom exchange. Although these Dicke states are individually entangled, a general solution for the dynamics is a mixture of them. Currently, extensive works have identified explicit new solutions of the dynamics [10–12], yet asserting their separability is considered NP-hard [5, 6]. Thus, given the inherited complexity of the task, it has been conjectured, and numerically shown only for a few emitters, that no entanglement should be generated [7, 8].

In this work, we introduce a decomposition of the time-dependent state into a mixture of coherent spin states (CSS). CSS are product states, implying that, as long as the mixture is positive, no entanglement is present. We develop a numerically efficient technique to find positive mixtures — which we test up to N = 50. This insight is used to devise a new entanglement-optimized quantum trajectory algorithm that leads to unravelings of dynamics into almost perfect CSS trajectories. These results pave the way for efficient large-scale simulations of open systems with collective processes.

**Superradiant Dicke cascade**—In the idealized Dicke problem, N two-level emitters collectively decay from an excited state $|e\rangle$ to a ground state $|g\rangle$. The decay dynamics is described by the Lindblad master equation

$$ \frac{d}{dt}\hat{\rho} = \Gamma \left[ 2\hat{L}\hat{\rho}\hat{L}^{\dagger} - \hat{\rho}\hat{L}^{\dagger}\hat{L} - \hat{L}^{\dagger}\hat{L}\hat{\rho} \right], \quad (1) $$

<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Superradiant decay represented as a fall through the collective Bloch sphere. This leads to mixed states (grey points) on the $S_z$-axis. Our ansatz consists in representing those states as mixtures of coherent spin states (CSS), Eq. (3), here sketched for N = 2 (12 arrows pointing to the surface of the sphere). (b) Evolution of the coefficients $\rho_m$ in the entangled Dicke state expansion (2) (N = 5). (c) Equivalent positive evolution of expansion coefficients $P_a$ in the CSS ansatz (3), which demonstrates the absence of entanglement.</figcaption></figure>

with the collective jump operator $\hat{L} = \sum_i \hat{\sigma}_i / \sqrt{2N}$, $\hat{\sigma}_i = |g\rangle\langle e|_i$ the single-emitter lowering operators, and $\Gamma$ the collectively enhanced decay rate. We consider a fully inverted initial state $\hat{\rho}_0 = \bigotimes_{i=1}^N |e\rangle\langle e|_i$ [13]. In our idealized setup of Eq. (1), no Hamiltonian is present. Instead, all emitters couple identically to the same mode. This model was originally envisioned as arising from a densely packed set of two-level emitters with distances very small compared to the wavelength of the transition, neglecting coherent Hamiltonian contributions [2, 14]. In more practical setups, such an idealized model can be realized using atoms in optical cavity systems [15, 16], near 1D waveguides [17], or even with atomic nuclei excited by free-electron lasers [18], to name just a few examples.

Permutation symmetry allows for numerical simulations of very large systems [19] and enables analytical
insight [10, 11]. A common way to solve the dynamics is to restrict evolution to the symmetric Dicke basis $|m\rangle$,

$$ \hat{\rho}(t) = \sum_{m=0}^{N} \rho_m(t) |m\rangle\langle m|. \quad (2) $$

For $N$ emitters, $m = 0, \dots, N$ is the number of emitters in the excited state $|e\rangle$, and $|m\rangle$ is the symmetrized wave-function of all such states [7, 11]. Dicke states are collective spin states of a spin-$S$ model with $S = N/2$, and the spin-lowering operator acts on the Dicke basis as $\hat{S}^- |m\rangle = \sqrt{S(S+1) - (S-m)(S+1-m)} |m-1\rangle$. In our convention, $\hat{L} = \hat{S}^- / \sqrt{2N}$. Hence, the dynamics can be visualized in a collective Bloch sphere, where the decay from the fully inverted state evolves through the center of the sphere [Fig. 1 (a)]. An evolution of Dicke state probabilities $\rho_m(t)$ is shown in Fig. 1(b) for $N = 5$.

**Entanglement**—Dicke states with $m \in [1, N-1]$ are entangled [20–24], and for any bipartition of emitters into two blocks $A$ and $B$, they feature a finite entanglement entropy. The von Neumann entanglement entropy is defined as $S_{VN}(|\psi\rangle) = -\sum_{\beta}^{D_B} s_{\beta} \log_2(s_{\beta})$ for a pure state $|\psi\rangle$. Here, $s_{\beta}$ are the eigenvalues of the reduced density matrix in block $B$, obtained from the partial trace $\hat{\rho}_B = \text{tr}_A(|\psi\rangle\langle\psi|)$, and $D_B$ is the Hilbert space dimension of block $B$. For a single Dicke state and a splitting of emitters into two blocks with $N_B$ (block $B$) and $N-N_B$ (block $A$), one finds $s_{\beta} = \binom{N_B}{\beta} \binom{N-N_B}{m-\beta} / \binom{N}{m}$ and $D_B = N_B + 1$ [20, 25]. For a symmetric state $S_{VN}$ is the largest for $N_B = N/2$, we thus focus below on this splitting. The entanglement in the half-excited Dicke state is proportional to $S_{VN}(|m=N/2\rangle) \propto \log(N)$ [20, 21], growing slowly with $N$ due to the symmetry-restricted sub-Hilbert space size $D_B = N/2 + 1 \ll 2^N$ for large $N$.

For a statistical mixture of pure states $|\psi^{[\chi]}\rangle$, $\hat{\rho} = \sum_{\chi} q_{\chi} |\psi^{[\chi]}\rangle\langle\psi^{[\chi]}|$ with $q_{\chi} > 0$ and $\sum_{\chi} q_{\chi} = 1$, we define an averaged trajectory entanglement (TE) as $\bar{S}_{QT}(\hat{\rho}) = \sum_{\chi} q_{\chi} S_{VN}(|\psi^{[\chi]}\rangle)$. For example, when computing super-radiant decay in the Dicke expansion [Eq. (2)] one finds that TE is finite, $\bar{S}_{QT}(\hat{\rho}(t)) = \sum_m \rho_m(t) S_{VN}(|m\rangle) > 0$, except for $t = 0$ and $t \to \infty$, see Fig. 3. Nevertheless, TE is not connected to genuine separability [26, 27] since the decomposition into pure states is highly non-unique. A proper definition of entanglement is given by the entanglement of formation [28], which is defined as the minimum TE of all pure-state ensembles representing the density matrix. Finding such a minimum is generally considered NP-hard [29, 30]. We will now introduce an alternative expansion into well-chosen CSS with vanishing TE, thus zero entanglement of formation at all times.

**Separable CSS ansatz**—Let us consider the expansion

$$ \hat{\rho}(t) = \sum_{a=0}^{N} P_a(t) \left[ \frac{1}{N_{\phi}} \sum_{b=1}^{N_{\phi}} |\theta_a(t), \phi_b\rangle\langle\theta_a(t), \phi_b| \right], \quad (3) $$

where $|\theta_a(t), \phi_b\rangle$ denotes a CSS parametrized by a (time-varying) polar angle $\theta_a(t)$ and an azimuthal angle $\phi_b$. A CSS is a product state, $|\theta_a(t), \phi_b\rangle = \bigotimes_j [\cos\theta_a(t) |g\rangle_j + \exp(i\phi_b) \sin(\theta_a) |e\rangle_j]$ [4, 31], represented by points on the Bloch sphere surface. In Eq. (3) the angles $\phi_b$ are chosen uniformly such that on average the $S_x$ and $S_y$ components vanish [see Fig. 1(a)]. We choose $N_{\phi} = 2N$ angles $\phi_b = b(\pi/N) = b\delta\phi$ ($b = 0, \dots, N_{\phi}$), but note that larger $N_{\phi} = xN$ with integer $x > 2$ are also possible. The vector $P_t$ with elements $P_a(t)$ represents the time-dependent state as a distribution of polar angles. Note that our ansatz is reminiscent of a CSS phase-space description [32]. There, states in the Dicke cascade are described by a time-dependent Wigner function, similar to $P_a(t)$. A Wigner function, however, can become negative, while here we aim at finding $P_a(t) \ge 0$.

**Positive solutions with the CSS ansatz**—In the Dicke state expansion (2), a solution to Eq. (1) at time $t$ is described by a $(N+1)$-dimensional vector $\rho(t)$ with elements $\rho_m(t)$. It is given by $\rho_t = A_t[;, 1]$, where $A_t[;, 1]$ denotes the first column of the evolution matrix $A_t = \exp(Dt)$. The $(N+1) \times (N+1)$ matrix $D$ has a bidagonal form [10], $(D)_{mn} = -\delta_{mn}\beta_{N-m}^2 + \delta_{m(n-1)}\beta_{N-m}^2$, for $m, n = 0, \dots, N$, with $\beta_m^2 = m(N - m + 1)/N$ and $\delta_{mn}$ the Kronecker delta. Note that $A_t$ is a stochastic matrix in the sense of Markov chains, satisfying the conditions $0 \le (A_t)_{mn} \le 1$ and $\sum_{m=0}^N (A_t)_{mn} = 1, \forall n$. A key feature of a mapping with a stochastic matrix is that $\rho_t$ remains positive and normalized, $\sum_m \rho_m(t) = 1$. Thus, $\rho_t$ remains a valid probability distribution during the dynamics.

We introduce a mapping between the representations,

$$ \rho_t = M_t \cdot P_t \qquad P_t = (M_t)^{-1} \cdot \rho_t, \quad (4) $$

$$ (M_t)_{mn} = \binom{N}{m} z_n(t)^{N-m} [1 - z_n(t)]^m. \quad (5) $$

Here, we define a time-dependent mapping matrix $M_t$ with $N + 1$ time-dependent functions $z_a(t) = \cos^2(\theta_a(t)/2) \in [0, 1]$ [25]. Note that $M_t$ is a stochastic matrix, guaranteeing that $\rho_t$ is positive if $P_t > 0$. However, the inverse matrix $(M_t)^{-1}$ is not necessarily stochastic. To prove that $P_t$ remains positive, one can show that the $P$-evolution matrix, $(M_t)^{-1} \cdot A_t$, is stochastic. It is possible to achieve this iteratively: Defining a discrete time index $\tau$, such that $t = \tau\Delta t$, and $A^{[\Delta t]} = \exp(D\Delta t)$, such that $A_t = (A^{[\Delta t]})^\tau$, in first order it can be shown that $(M_t)^{-1} \cdot A^{[\Delta t]} \cdot M_t$ preserves the positivity of $P$ for $\Delta t \to 0$ (see companion article [33]).

Let us now construct an explicit positive solution for $P_t$, which requires solving a system of linear equations

$$ M_t \cdot P_t = A_t[;, 1]. \quad (6) $$

Our objective is to solve this linear system with a positive vector $P_t$. Then, since $A_t[;, 1]$ is a solution to Eq. (1) in
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) Logarithm of negativity $\log_{10}(\mathcal{N})$ in the range $[-10, 0]$ as function of time and $\eta$, for $N = 30$. Before the burst time $t_B = \ln(N)/\Gamma$ (dashed lines) many positive solutions exist. Insets: Zooms in the ranges $0 \le t\Gamma \le 0.5$, $10^{-3} \le \eta \le 0.7$ and $9 \le t\Gamma \le 9.02$, $0.97332 \le \eta \le 0.97344$. After $t_B$, two passages to late times with positive solutions are highlighted: star (lower, at all time remaining $\eta \le 1$), triangle (upper). (b) $\log_{10}(\mathcal{N})$ as function of $\eta$ at $t\Gamma = 8$. For different $N$, $\mathcal{N} \to 0$ in the areas of the two passages. (c) Evolution with $P_a > 0$ for a numerically extracted $\eta(t)$ along the lower passage ($N = 30, \mathcal{N} < 10^{-6}$). (d) $\{P_a\}_{a=0}^N$ distribution as function of $\theta$ at different times ($N = 30$). At $t_B$, the peak centers around $\theta = \pi/2$.</figcaption></figure>

the Dicke state representation at time $t$, we have shown that there exists a positive solution to the Dicke dynamics also in the CSS basis (3). The transformation matrix $M_t$ depends on $N+1$ parameters, finding a positive $P_t$ is thus a multi-parameter optimization problem. We now introduce a key step that can turn it into a simpler single-parameter problem: We restrict $M_t$ to matrices with equally spaced polar angles $\theta_a(t) = \eta(t)a\pi/N$, using a single parameter $\eta(t)$. $M_t$ then belongs to the family of centrosymmetric matrices. For such matrices the inverses always exist [34], and $\eta(t)$ can in principle be an arbitrary positive number, corresponding to a set of CSS angles $\theta$ (interpreted in the range $[0, \pi]$). The challenge is to find an $\eta(t)$ that leads to positive $P$ at all times.

Note that $\eta(t)$ does not need to be a smooth function, and multiple positive-$P$ solutions can exist. However, for $N=2$, we identify the following analytical one [25]:

$$ \eta(t) = \frac{2}{\pi} \arccos \sqrt{\frac{t\Gamma e^{-t\Gamma}}{2(1 - e^{-t\Gamma}) - t\Gamma e^{-t\Gamma}}} \quad (7) $$

Obtaining $\eta(t)$ analytically for larger $N$ requires solving a set of transcendental inequalities [25].

To find positive $P_t$ numerically, we introduce the negativity $\mathcal{N} = -\sum_a \min(P_a, 0)$. In Fig. 2(a) we plot $\log_{10}(\mathcal{N})$ as a function of time $t$ and $\eta$ for $N = 30$. Note that we cut off values $\mathcal{N} < 10^{-10}$ and $\mathcal{N} > 1$ in order to avoid showing effects of machine precision noise and to ignore regions with diverging negativity. Below the burst time, $t_B = \ln(N)/\gamma$, $0 < t \le t_B$ we find a large parameter range (blue area) with many positive solutions, extending to $t \to 0$ (left inset). Note that at $t = 0$, $P_a = \delta_{a0}$ is a trivial solution. For $t > t_B$, the number of solutions reduces, but we still find several positive “passages” to late times. We highlight two passages with a star (lower) and a triangle (upper). In general, passages are very narrow, but zooming into the respective regions (right inset) reveals negativities that can be made arbitrarily small. This point can be better appreciated in

Fig. 2(b), where a slice of $\mathcal{N}$ is shown, as a function of $\eta$ at time $t\Gamma = 8$. We observe $\mathcal{N} \to 0$ at the locations of the passages (up to $N = 50$, where noisy behavior due to machine precision limitations starts to appear). Although the upper passage has some internal structure, it appears that the lower passage corresponds to a single solution with $\eta(t) \le 1$ at all times [25]. Note that for $N = 2$ the lower passage corresponds exactly to Eq. (7). Due to the narrow nature of the passage, it is difficult to find a smooth fitting function with low $\mathcal{N}$. We can however extract solutions by minimizing $\mathcal{N}$ numerically along the lower passage [25]. We note that finding a smooth solution of $\eta(t)$ numerically that fully follows the lower passage becomes more difficult with $N$ due to machine precision limitations in the linear solver. However, solutions that are allowed to jump between the passages are generally much easier to find. In addition, we have observed the existence of positive solutions for values of $\eta(t) > 1$. For a solution along the lower passage, in Fig. 2(c) we show a corresponding positive evolution of $P_a$ as a function of time (for $N = 30$). We always observe that at $t_B$, $\max(P_{t_B})$ is lowest in the evolution [25]. In Fig. 2(d) we plot the distribution of $P_a$ at various $t$, showing that $\max(P_{t_B})$ corresponds to $\theta = \pi/2$. This asserts that CSS populations are distributed around the equator of the Bloch sphere at the time of the burst. Having established that a positive CSS decomposition exists, we now introduce an algorithm finding trajectory unravelings close to CSS.

Quantum trajectory unraveling—Our quantum trajectory (QT) method emerges from a Kraus operator formalism [26, 35]. Markovian dynamics allows for a discretized evolution with timesteps $\Delta t \to 0$. Each step can be formulated as a quantum operation $\mathcal{E}(\bullet)$ with operator sum representation $\hat{\rho}_{\text{out}} = \mathcal{E}(\hat{\rho}_{\text{in}}) = \sum_k \hat{E}_k \hat{\rho}_{\text{in}} \hat{E}_k^\dagger$ [36]. For Eq. (1) there are two Kraus operators

$$ \hat{E}_0 = 1 - \Delta t \Gamma \hat{L}^\dagger \hat{L}, \quad \hat{E}_1 = \sqrt{2 \Delta t \Gamma} \hat{L}. \quad (8) $$

In QT algorithms, the time-dependent density ma-
trix $\hat{\rho}_t$ is approximated by a large number $n_T$ of pure-state samples, $\hat{\rho}_t \approx \sum_{\chi=1}^{n_T} |\psi_t^{[\chi]}\rangle\langle\psi_t^{[\chi]}|/n_T$, with the initial condition $|\psi_0^{[\chi]}\rangle = |m=N\rangle$ for all $\chi = 1, \dots, n_T$. A single timestep, which evolves a trajectory from $|\psi_t^{[\chi]}\rangle$ to $|\psi_{t+\Delta t}^{[\chi]}\rangle$, is stochastically simulated by applying $\hat{E}_k$ with probability $p_k = \langle\psi_t^{[\chi]}|\hat{E}_k^\dagger\hat{E}_k|\psi_t^{[\chi]}\rangle$, followed by a renormalization, $|\psi_{t+\Delta t}^{[\chi]}\rangle = \hat{E}_k|\psi_t^{[\chi]}\rangle/\sqrt{p_k}$. As a result, the statistical approximation of $\hat{\rho}_{t+\Delta t}$ is exact for $n_T \to \infty$, since $\mathcal{E}(|\psi_t^{[\chi]}\rangle\langle\psi_t^{[\chi]}|) = \sum_k p_k \hat{\rho}_k$ with $\hat{\rho}_k = \hat{E}_k|\psi_t^{[\chi]}\rangle\langle\psi_t^{[\chi]}|\hat{E}_k^\dagger/\langle\psi_t^{[\chi]}|\hat{E}_k^\dagger\hat{E}_k|\psi_t^{[\chi]}\rangle$ and $p_k$ a probability distribution. A quantum operation acting on $|\psi_t^{[\chi]}\rangle\langle\psi_t^{[\chi]}|$ thus amounts to randomly replacing the state $|\psi_t^{[\chi]}\rangle\langle\psi_t^{[\chi]}|$ with $\hat{\rho}_k$. In the limit $\Delta t \to 0$, this approach is equivalent to the standard QT method developed in quantum optics [37–40], which has been key to efficient simulations of open quantum many-body systems [41].

Naively performing QT with Kraus operators from Eq. (8) restricts all trajectories to the Dicke basis and thus to the $S_z$-axis of the Bloch sphere. The inset in Fig. 3(a) confirms this, showing an evolution of the normalized Bloch vector length, $\xi_\chi = 2\langle\psi_t^{[\chi]}|\hat{S}^2|\psi_t^{[\chi]}\rangle^{1/2}/N$, from $1 \to 0 \to 1$ (we denote a trajectory average by $\bullet$). Using the naive unraveling also implies that trajectories become entangled. Around $t_B$, with most of the population in states $|m \approx N/2\rangle$, TE is thus of the order $\bar{S}_{QT} \sim \log N$ [20, 21], as confirmed in Fig. 3(b). From a classical computational point of view, this is inefficient, as above we have shown that TE can be zero. Let us now present a method to drastically reduce TE.

**Low-entanglement QT unraveling**—We exploit the unitary degree of freedom of Kraus operators [26, 35],

$$ \hat{F}_n = \sum_k u_{nk} \hat{E}_k. \qquad (9) $$

Here, $\hat{F}_n$ ($n, k = 0, 1$) are Kraus operators defining the same quantum operation as $\hat{E}_k$, as long as the matrix with elements $(\mathbf{u})_{nk} = u_{nk}$ is unitary [36]. Any observable evolution can be computed in QT with Kraus operators $\hat{F}_n$ instead of $\hat{E}_k$. Entanglement entropy, however, is not an observable but a nonlinear function of the state. As such, it can substantially depend on the choice of $\mathbf{u}$ [26, 35]. Following [26] we focus on unitary matrices parametrized with $\theta_F \in [0, \pi/2]$ and $\phi_F \in [0, 2\pi]$:

$$ \mathbf{u}(\theta_F, \phi_F) = \begin{pmatrix} \cos \theta_F & \sin \theta_F \\ -\sin \theta_F & \cos \theta_F \end{pmatrix} \begin{pmatrix} e^{i\phi_F} & 0 \\ 0 & e^{-i\phi_F} \end{pmatrix}. \quad (10) $$

In the limit of $N \to \infty$, this parametrization leads to

$$ \hat{F}_n \to e^{(-1)^n \exp(-i2\phi_F) \cot(\theta_F)^{2n-1} \sqrt{2\Delta t} \hat{L}}. \quad (11) $$

A Bloch vector rotation can be performed with $\hat{R}(\theta, \phi) = e^{\exp(i\phi) \tan(\theta/2)} \hat{S}^{-}$ [31]. Therefore, we observe that for

<figure><img src="image_4.png" /><figcaption>FIG. 3. Trajectory averaged evolution ($N = 50$) of (a) normalized Bloch vector length $\bar{\xi}$ and (b) TE. We compare naive [Kraus operators (8)], $\phi$-randomized, and $\phi$-optimized unravelings [using randomized/optimized choices of Kraus operators (9), see text]. The scaling with $N$ is analyzed for averaged: (c) minimum Bloch vector length, $\xi_\chi$; and (d) maximum TE, $\bar{S}_{max}$. Dashed lines highlight the burst time.</figcaption></figure>

large $N$, applying a Kraus operator $\hat{F}_n$ to a CSS (without renormalization) becomes equivalent to mapping it to another CSS. Note that the polar rotation induced by the Kraus operator depends on $\sqrt{\Delta t}$ and $\theta_F$. We find that the precise choice of $\theta_F$ does not have any major impact on TE [25] (except close to the trivial $\theta_F = 0, \pi$). We therefore fix $\theta_F = \pi/4$ while varying $\phi_F$ only. Note that for $N \to \infty$, $\phi_F \to \phi/2$. In ansatz (3), the $\phi_b$ are homogeneously distributed. Therefore, to lower TE we propose to randomly select $\phi_F$ uniformly at each quantum operation. In Fig. 3 we demonstrate that this $\phi$-randomized algorithm leads to a strong reduction of TE. In addition, the trajectory Bloch vectors remain now almost on the sphere as $\bar{\xi} \to 1$. With increasing $N$, spin-lengths further increase while TE decreases, as shown in Fig. 3(c/d). There we plot the averaged maximum (minimum) von Neumann entanglement entropy $\bar{S}_{max}$ (spin length $\bar{\xi}_{max}$), taken over time, [$S_{max} = \max_t(S_{VN}(|\psi_t^{[\chi]}|))$, $\xi_{min} = \min_t(\xi_\chi)$]. This is in line with the expectation from Eq. (11).

We can further reduce TE by resorting to an adaptive minimization [26, 35]. Defining $|\psi_{po}^{[\chi]}\rangle_{0,1} = \hat{F}_{0,1}|\psi_t^{[\chi]}\rangle$, the post-operation entanglement is $S_{po}(\phi_F) = [p_0S_{VN}(|\psi_{po}^{[\chi]}\rangle_0) + p_1S_{VN}(|\psi_{po}^{[\chi]}\rangle_1]/2$. Using expressions from [25], we use $S_{po}(\phi_F)$ as a cost function for numerical minimization at each step. After finding $\phi_F^o$ minimizing $S_{po}(\phi_F)$, we use a QT step with the corresponding $\hat{F}_{0,1}(\phi_F^o)$. We find that this $\phi$-optimized algorithm leads to TE below $10^{-5}$ for $N = 50$ before $t_B$, and always two orders of magnitude smaller than the naive unraveling [Fig. 3(b)]. Now, $1 - \bar{\xi}_{min} \lesssim 10^{-3}$ for $N = 50$ and $\bar{\xi}_{min} \to 1$ with increasing $N$, indicating an convergence
towards CSS trajectories for large $N$. In Fig. 3(d) we observe no $N$ dependence on the residual $S_{max}$. Note that for $N > 300$ also in the $\phi$-randomized algorithm the scaling with $N$ becomes sub-linear, e.g. for $N = 800$ we find $\bar{S}_{max} = 0.0612 \pm 0.003$ (not plotted). We attribute this to the renormalization step in QT. We also note that the $\phi$-optimization is not guaranteed to find zero TE due to the “greedy” nature of the algorithm, only considering knowledge of a single trajectory at a time. Nevertheless, it is important to highlight that such entanglement entropies are extremely small. In practical scenarios, these low entropies will allow us to describe states almost as products, enabling the use of MPS with extremely small bond dimensions $\chi$ [42, 43]. While the bond dimensions usually depend exponentially on the entanglement entropies ($\chi \sim 2^{S_{max}}$), here $\chi \approx 2$ would already provide a close-to-perfect state approximation. Our unraveling scheme will therefore grant the removal of nearly all trajectory entanglement contributions from symmetric Dicke decay process, and can thus drastically improve the efficiency of simulations of quantum dynamics in more general systems where such processes are present [44–47].

**Conclusion & Outlook**—We have shown that states in Dicke superradiant decay can be written as positive mixtures of CSS [Eq. (3)]. This implies that there is no entanglement in the Dicke decay dynamics, i.e., the entanglement of formation is vanishing. The positive expansion coefficients $P_a$ ($a = 0, \dots, N$) correspond to weights of CSS with polar angles $\theta_a$ and uniformly distributed azimuthal angles. Crucially, only the spacing between different $\theta_a$ needs to evolve time-dependently [single parameter, $\eta(t)$] to find positive $P_a > 0$. We then introduced an entanglement-optimized numerical QT unraveling algorithm, based on a simple random variation of a single parametrization angle in the Kraus representation of a discretized timestep. In the large $N$ limit we have found that this reproduces almost perfect CSS trajectories.

Identifying a closed analytical form for $\eta(t)$ for arbitrary $N$ would be particularly interesting. Indeed, the insight that TE can be transformed away paves a way to new efficient algorithms, also in cases where, besides symmetric superradiant decay, Hamiltonian or dissipative processes without permutation symmetry are present. Here, combining $\phi$-randomization with matrix product states [42, 43] can lead to entanglement-optimized methods [26, 35, 48] for simulating collective dissipation on much larger scales than currently possible.

For future research, considering scenarios beyond the fully-symmetric Dicke limit has now become more tractable with our proposed method. Collective asymmetric jump operators can drive a separable quantum state into an entangled one in atomic systems with dipole-dipole interactions in free space or close to waveguides, for instance. These systems can generate multipartite entanglement, even without a coherent contribution of the master equation [49]. For some particular configurations, this entanglement has been shown to be irrelevant, allowing, for example, the correct computation of two-point local observables in a large array of particles [50]. With our method, such questions can now be explored on a full-state level: How much entanglement can be imprinted by a coherent contribution of the master equation or by asymmetric super- and sub-radiant decay channels?

**Acknowledgments**—We thank Claudiu Genes and Ruben Daraban for valuable discussions. We especially thank Nico Bassler for constructing the mathematical proof of separability in the related companion article [33]. This work was supported by the ERC Consolidator project MATHLOCCA (Grant nr. 101170485), by the CNRS through the EMERGENCE@INC2024 project DINOPARC, and by the French National Research Agency under the Investments of the Future Program project ANR-21-ESRE-0032 (aQCess). Computations were carried out using resources of the High Performance Computing Center of the University of Strasbourg, funded by Equip@Meso (as part of the Investments for the Future Program) and CPER Alsacacul/Big Data. P.R., A.C. and R.B. acknowledge the financial support of the São Paulo Research Foundation (FAPESP) (Grants No. 2022/12382-4, 2024/05564-4, 2022/06449-9, 2023/07463-8, 2023/14832-0, 2022/00209-6 and 2023/03300-7). L.O.R.S. acknowledges the financial support of CAPES-COFECUB (CAPES, Grants Nos. 88887.711967/2022-00).

[1] R. H. Dicke, Coherence in Spontaneous Radiation Processes, Phys. Rev. 93, 99 (1954).

[2] M. Gross and S. Haroche, Superradiance: An essay on the theory of collective spontaneous emission, Phys. Rep. 93, 301 (1982).

[3] J. C. MacGillivray and M. S. Feld, Theory of superradiance in an extended, optically thick medium, Phys. Rev. A 14, 1169 (1976).

[4] F. T. Arecchi, E. Courtens, R. Gilmore, and H. Thomas, Atomic coherent states in quantum optics, Phys. Rev. A 6, 2211 (1972).

[5] J. Tura, A. Aloy, R. Quesada, M. Lewenstein, and A. Sanpera, Separability of diagonal symmetric states: a quadratic conic optimization problem, Quantum 2, 45 (2018).

[6] G. Tóth and O. Gühne, Separability criteria and entanglement witnesses for symmetric quantum states, Applied Physics B 98, 617–622 (2009).

[7] E. Wolfe and S. F. Yelin, Certifying separability in symmetric mixed states of $n$ qubits, and superradiance, Phys. Rev. Lett. 112, 140402 (2014).

[8] N. Yu, Separability of a mixture of Dicke states, Phys. Rev. A 94, 060101 (2016).

[9] L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Entanglement
lement in many-body systems, Rev. Mod. Phys. 80, 517 (2008).

[10] R. Holzinger, N. S. Bassler, J. Lyne, F. G. Jimenez, J. T. Gohsrich, and C. Genes, Solving Dicke superradiance analytically: A compendium of methods, arXiv:2503.10463 10.48550/arXiv.2503.10463 (2025).

[11] R. Holzinger and C. Genes, Exact solution for Dicke superradiance, arXiv:2409.19040 10.48550/arXiv.2409.19040 (2024).

[12] D. Malz, R. Trivedi, and J. I. Cirac, Large-$n$ limit of dicke superradiance, Phys. Rev. A 106, 013716 (2022).

[13] C. T. Lee, Exact solution of the superradiance master equation. i. complete initial excitation, Phys. Rev. A 15, 2019 (1977).

[14] R. H. Lehmberg, Radiation from an $N$-Atom System. I. General Formalism, Phys. Rev. A 2, 883 (1970).

[15] M. A. Norcia, M. N. Winchester, J. R. K. Cline, and J. K. Thompson, Superradiance on the millihertz linewidth strontium clock transition, Science Advances 2, e1601231 (2016), https://www.science.org/doi/pdf/10.1126/sciadv.1601231.

[16] A. Angerer, K. Streltsov, T. Astner, S. Putz, H. Sumiya, S. Onoda, J. Isoya, W. Munro, K. Nemoto, and J. Majer, Superradiant emission from colour centres in diamond, Nature Physics 14 (2018).

[17] C. Liedl, F. Tebbenjohanns, C. Bach, S. Pucher, A. Rauschenbeutel, and P. Schneeweiss, Observation of superradiant bursts in a cascaded quantum system, Phys. Rev. X 14, 011020 (2024).

[18] A. Chumakov, A. Baron, I. Sergueev, C. Strohm, O. Leupold, Y. Shvyd'ko, G. Smirnov, R. Rüffer, Y. Inubushi, M. Yabashi, K. Tono, T. Kudo, and T. Ishikawa, Superradiance of an ensemble of nuclei excited by a free electron laser, Nature Physics 14 (2018).

[19] N. Shammah, S. Ahmed, N. Lambert, S. De Liberato, and F. Nori, Open quantum systems with local and collective incoherent processes: Efficient numerical simulations using permutational invariance, Phys. Rev. A 98, 063815 (2018).

[20] J. I. Latorre, R. Orús, E. Rico, and J. Vidal, Entanglement entropy in the Lipkin-Meshkov-Glick model, Phys. Rev. A 71, 064101 (2005).

[21] J. Schachenmayer, B. P. Lanyon, C. F. Roos, and A. J. Daley, Entanglement Growth in Quench Dynamics with Variable Range Interactions, Phys. Rev. X 3, 031015 (2013).

[22] G. Tóth, Detection of multipartite entanglement in the vicinity of symmetric Dicke states, J. Opt. Soc. Am. B, JOSAB 24, 275 (2007).

[23] L.-M. Duan, Entanglement Detection in the Vicinity of Arbitrary Dicke States, Phys. Rev. Lett. 107, 180502 (2011).

[24] B. Lücke, J. Peise, G. Vitagliano, J. Arlt, L. Santos, G. Tóth, and C. Klempt, Detecting Multiparticle Entanglement of Dicke States, Phys. Rev. Lett. 112, 155304 (2014).

[25] See Supplemental Material at [URL will be inserted by publisher] for details on: entanglement entropies of a general superposition of Dicke states; mapping matrix between Dicke and CSS representation; exact expression of $\eta(t)$ for $N = 2$; additional data.

[26] R. Daraban, F. Salas-Ramírez, and J. Schachenmayer, Non-unitarity maximizing unraveling of open quantum dynamics, SciPost Phys. 18, 048 (2025).

[27] G. Preisser, D. Wellnitz, T. Botzung, and J. Schachenmayer, Comparing bipartite entropy growth in open-system matrix-product simulation methods, Physical Review A 108, 10.1103/physreva.108.012616 (2023).

[28] C. H. Bennett, D. P. DiVincenzo, J. A. Smolin, and W. K. Wootters, Mixed-state entanglement and quantum error correction, Physical Review A 54, 3824 (2024).

[29] L. Gurvits, Classical deterministic complexity of edmonds' problem and quantum entanglement, in Proc. 35th Annual ACM Symp. on Theory of Computing, STOC '03 (Association for Computing Machinery, New York, NY, USA, 2003) p. 10–19.

[30] S. Gharibian, Strong np-hardness of the quantum separability problem, Quantum Info. Comput. 10, 343–360 (2010).

[31] A. M. Perelomov, Generalized coherent states and some of their applications, Soviet Physics Uspekhi 20, 703 (1977).

[32] L. M. Narducci, C. A. Coulter, and C. M. Bowden, Exact diffusion equation for a model for superradiant emission, Phys. Rev. A 9, 829 (1974).

[33] N. S. Bassler, Absence of entanglement growth in Dicke superradiance, arXiv (2025), 2504.XXXXX.

[34] I. J. Good, The inverse of a centrosymmetric matrix, Technometrics 12, 925 (1970).

[35] T. Vovk and H. Pichler, Quantum trajectory entanglement in various unravelings of Markovian dynamics, arXiv:2404.12167 10.48550/arXiv:2404.12167 (2024).

[36] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cambridge, England, UK, 2010).

[37] J. Dalibard, Y. Castin, and K. Mølmer, Wave-function approach to dissipative processes in quantum optics, Physical Review Letters 68, 580 (1992).

[38] R. Dum, P. Zoller, and H. Ritsch, Monte Carlo simulation of the atomic master equation for spontaneous emission, Physical Review A 45, 4879 (1992).

[39] H. Carmichael, An Open Systems Approach to Quantum Optics (Springer, Berlin, Germany, 1993).

[40] C. W. Gardiner, A. S. Parkins, and P. Zoller, Wave-function quantum stochastic differential equations and quantum-jump simulation methods, Physical Review A 46, 4363 (1992).

[41] A. J. Daley, Quantum trajectories and open many-body quantum systems, Adv. Phys. (2014).

[42] N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, Entropy Scaling and Simulability by Matrix Product States, Phys. Rev. Lett. 100, 030504 (2008).

[43] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011).

[44] A. González-Tudela and D. Porras, Mesoscopic entanglement induced by spontaneous emission in solid-state quantum optics, Phys. Rev. Lett. 110, 080502 (2013).

[45] S. Masson and A. Asenjo-Garcia, Universality of dicke superradiance in arrays of quantum emitters, Nature Communications 13 (2022).

[46] S. Cardenas-Lopez, S. J. Masson, Z. Zager, and A. Asenjo-Garcia, Many-body superradiance and dynamical mirror symmetry breaking in waveguide qed, Phys. Rev. Lett. 131, 033605 (2023).

[47] R. Holzinger and S. F. Yelin, Beyond dicke superradiance: Universal scaling of the peak emission rate (2025),
arXiv:2506.12649 [quant-ph].

[48] T. Vovk and H. Pichler, Entanglement-Optimal Trajec-
tories of Many-Body Quantum Markov Processes, Phys.
Rev. Lett. **128**, 243601 (2022).

[49] J. T. Lee, S. Cardenas-Lopez, S. J. Masson, R. Trivedi, and A. Asenjo-Garcia, Exact many-body quantum dynamics in one-dimensional baths via "superspins" (2025), arXiv:2505.00588 [quant-ph].

[50] X. H. H. Zhang, D. Malz, and P. Rabl, Unraveling super-radiance: entanglement and mutual information in collective decay (2025), arXiv:2505.13401 [quant-ph].
## Supplemental Material for:
Unraveling Dicke Superradiant Decay with Separable Coherent Spin States

P. Rosario¹,², L. O. R. Solak¹,², A. Cidrim¹, R. Bachelard¹, and J. Schachenmayer²

¹Departamento de Física, Universidade Federal de São Carlos,
Rodovia Washington Luís, km 235—SP-310, 13565-905 São Carlos, SP, Brazil

²CESQ/ISIS (UMR 7006), CNRS and Université de Strasbourg, 67000 Strasbourg, France

### I: Bipartite entanglement entropies in Dicke-state superpositions

Following [20, 21] a Dicke state for $N$ emitters $|N, m\rangle$ (with $m$ denoting the number of excited states out of $N$ emitters) can be written as the tensor product of two smaller sub-spaces, we define block A and B with sizes $N - N_B$ and $N_B$ respectively such that

$$|N, m\rangle = \sum_{l=0}^{N_B} \sqrt{p_{lm}} |N_B, l\rangle_A |N - N_B, m - l\rangle_B \quad (S1)$$

the amplitudes are given by

$$p_{lm} = \frac{\binom{N_B}{l} \binom{N-N_B}{m-l}}{\binom{N}{m}}. \quad (S2)$$

Here the states $|N_B, k\rangle_A$ and $|N - N_B, k\rangle_B$ are Dicke states of the emitters living in the ensembles A and B, respectively. Note that the sum in Eq. (S1) includes invalid Dicke states, e.g. when $m - l < 0$ and $m - l > N - N_B$, however for those cases the $p_{lm}$ are zero. Taking now an arbitrary Dicke state superposition, we can write

$$|\psi\rangle = \sum_m c_m |N, m\rangle \quad (S3)$$

$$= \sum_{m=0}^{N} c_m \sum_{l=0}^{N_B} \sqrt{p_{lm}} |N_B, l\rangle_A |N - N_B, m - l\rangle_B \quad (S4)$$

now taking the partial trace with respect to block A, we get

$$\hat{\rho}_B = \text{tr}_A (|\psi\rangle\langle\psi|) = \sum_{k=0}^{N_B} \sum_{m=0}^{N} \sum_{n=0}^{N} c_m c_n^* \sqrt{p_{kn}} \sqrt{p_{km}} |N - N_B, m - k\rangle_B \langle N - N_B, n - k|_B \quad (S5)$$

From a diagonalization of this matrix, the von Neumann entropy can be computed in the usual way.

### II: Mapping matrix between Dicke and CSS representation

Spin coherent states are characterized by a collection of individual particles pointing in the same direction, no entanglement is present. We formally write

$$|\theta(t), \phi\rangle = \bigotimes_{j=1}^{N} \left[ \cos\left(\frac{\theta(t)}{2}\right) |e\rangle + \sin\left(\frac{\theta(t)}{2}\right) e^{i\phi} |g\rangle \right]_j \quad (S1)$$

with $0 \le \theta(t) \le \pi$ and $0 \le \phi \le 2\pi$. In terms of Dicke states $\{|m\rangle\}_{m=0}^N$ with $m$ being the number of excitations, $|\theta(t), \phi\rangle$ can also be written as

$$|\theta(t), \phi\rangle = \sum_{k=0}^{N} \sqrt{\binom{N}{k}} \left(\cos\left(\frac{\theta(t)}{2}\right)\right)^{N-k} \left(\sin\left(\frac{\theta(t)}{2}\right)\right)^k e^{ik\phi} |N-k\rangle \quad (S2)$$

those states live on the surface of the generalized Bloch sphere. On the other hand, a $N$-qubit statistical mixture of Dicke states, is defined as

$$\hat{\rho}(t) = \sum_{m=0}^{N} \rho_m(t) |m\rangle\langle m| \quad (S3)$$
where for some specific populations distributions $\rho_m(t)$, the quantum state in Eq. (S3) can be full separable [8]. In the following, we show that a statistical mixture of Dicke states similar to Eq. (S3) can be obtained by combining multiple spin coherent states. For that, we propose the time-dependent ansatz

$$ \hat{\rho}(t) = \sum_{a=0}^{N_{\theta}} P_a(t) \left[ \frac{1}{N_{\phi}} \sum_{b=1}^{N_{\phi}} |\theta_a(t), \phi_b\rangle\langle\theta_a(t), \phi_b| \right] \quad (S4) $$

with $N_{\theta}$ and $N_{\phi}$ defining how many angles $\theta$ and $\phi$ are needed to reproduce the target state in Eq. (S3). Therefore, choosing $N_{\phi} = 2N$, $N_{\theta} = N$ and $\phi_b = \frac{2\pi}{N_{\phi}}b = \frac{\pi}{N}b$, we arrive to the statistical mixture of Dicke states

$$ \begin{aligned} \hat{\rho} &= \sum_{m=0}^{N} \left[ \sum_{a=0}^{N_{\theta}} P_a(t) \binom{N}{m} \left( \cos\left(\frac{\theta_a(t)}{2}\right) \right)^{2N-2m} \left( \sin\left(\frac{\theta_a(t)}{2}\right) \right)^{2m} \right] |N-m\rangle\langle N-m|, \\ &= \sum_{m=0}^{N} \left[ \sum_{a=0}^{N} P_a(t) \binom{N}{m} z_a(t)^{N-m} (1-z_a(t))^m \right] |N-m\rangle\langle N-m| \quad \text{with } z_a(t) = \cos^2\left(\frac{\theta_a(t)}{2}\right). \end{aligned} \quad (S5) $$

note that, the quantity inside the brackets can be written as the product of a stochastic matrix times a probability vector as $M \cdot P$, with

$$ P = [P_0(t), P_1(t), \dots, P_N(t)]^T \quad (M)_{mn} = \binom{N}{m} z_n(t)^{N-m} [1 - z_n(t)]^m, \quad m, n = 0, 1, \dots, N. \quad (S6) $$

Therefore, any probability distribution obtained from the product $M \cdot P$ corresponds to a full separable state iff $P \ge 0$.

### III: Exact expression of $\eta(t)$ for $N=2$

Defining the general coherent state representation

$$ P^{[\tau]} = (M^{[\tau]})^{-1} \cdot \rho^{[\tau]}, \quad (S1) $$

for a given $\rho^{[\tau]}$, we need to find $\eta(t)$ such that $\forall \tau \quad P^{[\tau]} \ge 0$. For the case $N=2$, taking $\theta_0 = 0$, $\theta_1 = \eta(t)\pi/2$ and $\theta_2 = \eta(t)\pi$ we get the set of inequalities (here, $\Gamma \equiv 1$)

$$ P_1 = \frac{1}{8} \left( 8e^t + (8 - 8e^t + 3t) \csc^2\left(\frac{\eta(t)\pi}{4}\right) + 2(-1 + e^t - t) \csc^4\left(\frac{\eta(t)\pi}{4}\right) - t \sec^2\left(\frac{\eta(t)\pi}{4}\right) \right) \ge 0 \quad (S2) $$

$$ P_2 = - \frac{e^{-t} \left(-1 + e^t - \frac{3}{2}t + \left(-1 + e^t - \frac{t}{2}\right) \cos(\eta(t)\pi)\right) \csc^4\left(\frac{\eta(t)\pi}{4}\right)}{2 \left(1 + 2 \cos\left(\frac{\eta(t)\pi}{2}\right)\right)} \ge 0 \quad (S3) $$

$$ P_3 = \frac{e^{-t} \csc^4\left(\frac{\eta(t)\pi}{4}\right) \left(-2 + 2e^t - t - t \sec^2\left(\frac{\eta(t)\pi}{4}\right)\right)}{8 \left(1 + 2 \cos\left(\frac{\eta(t)\pi}{2}\right)\right)} \ge 0 \quad (S4) $$

taking the mapping $y = \cos^2(\eta(t)\pi/2)$ and using Chebyshev polynomials properties, we arrive to the expression presented in the main text for $\eta(t)$. When $N > 2$, the number of inequalities increases like $N+1$ and obtaining $\eta(t)$ becomes analytically unfeasible. However, this case is numerically addressed in the main text.

### IV: Numerical considerations and additional data

Lower passage $\eta(t)$ solutions— At each time, $P$ can be obtained by solving the linear system Eq. (6) while varying the parameter $\eta(t)$. We find that in this way negativites along the evolution can be made arbitrary small up to numerical precision limits. We search for solutions of $\eta(t)$ following the lower passage. For the plots in this work we show corresponding evolutions of $P$ using values with negativities $N < 10^{-6}$. Note that in all cases for numerical stability we truncate probabilities in the vector $\rho_t$ below relative machine precision ($10^{-16}$). We show some selected
<figure><img src="image_5.png" /><figcaption>FIG. S1. (a) Numerically extracted positive-$P$ solutions close to the lower passage for different values of $N$. All solutions remain close to the analytical result for $N = 2$ (red line). The smooth behavior may suggest the existence of a general smooth solution for $N > 2$. (b) $P_a$ evolution for $N = 20$ and (c) $N = 40$ ($N < 10^{-6}$). The vertical dashed lines indicate the burst time.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. S2. QT simulations for $N = 50$ with randomized $\phi_F$ taking the average over 100 trajectories: (a) TE evolution using various choices of $\theta_F$. No major variation in TE is observed when comparing $\theta_F = \pi/8, \pi/4, 3\pi/4$; (b) TE evolution for different timesteps, showing invariance with respect to $\Delta t$</figcaption></figure>

evolution of $\eta(t)$ extracted numerically in Fig. S1(a). It is important to point out that there is only a little variation in the solutions for different $N$ and that the solutions remain rather close to the exact analytical curve for $N = 2$. In Fig. S1(b/c) we also show the evolution of $P$ for $N = 20$ and $N = 40$, respectively, together with the burst time $t_B$ as vertical dashed line. Note that not in all cases we found smooth $\eta(t)$ and $P_a(t)$, since multiple solutions exist and our numerical optimizer can jump between them.

*Independence of TE on the choice $\theta_F$ and $\Delta t$—* For the simulations related to transformed Kraus operators according to Eq. (9), in Fig. S2(a) we verify that selecting different values of $\theta_F$ have no major impact on the TE evolution. In addition we verify in Fig.S2(b) that our QT simulations are converged in the timestep $\Delta t$ and do not depend on $\Delta t$.