Laure Moiroux-Arvis
email: <laure.moiroux-arvis@inrae.fr>

François Christophe Cariou

François Pinet

Jean-Pierre Chanet

CIDEA : Robot behavior adaptation from interactive communication with buried sensor nodes. Application to smart agriculture

Keywords: Wireless Underground Sensor Networks, Soil moisture, Mobile robotics, LoRa, Agriculture

come

I. INTRODUCTION

In recent years, the development of mobile robotics has increased significantly in agriculture, whether it be to increase farm productivity, perform environmentally friendly operations or relieve human operators from tedious and unhealthy work [START_REF] Oliveira | Advances in agricultural robotics: a state-of-the-art review and challenges ahead[END_REF]. Many challenges remain, however, to be faced to obtain fully operational robots. One of them is to provide for the robot the capacity to take into account the ground conditions in the field to adapt its behavior accordingly. In particular, when the soil moisture varies significantly, the speed of the robot and the tuning of its agricultural tool have to be adapted, see the illustration in Figure 1.

Figure 1. Variability of the soil moisture in a field requiring to adapt the working speed of the robot and its tool accordingly Nevertheless, the access to the information of soil moisture, locally and in real time, is a real issue. Microclimate weather stations and connected devices can be positioned at the edge of the field, but they are not sufficient to obtain local and accurate information of the soil moisture within the areas [START_REF] Tenzin | Low cost weather station for climate-smart agriculture[END_REF]. In addition, such systems can obviously not be positioned and multiplied inside the field as they represent obstacles for the farming activities.

The recent advances in Wireless Underground Sensor Networks (WUSNs) and Internet of Underground Things (IoUT) have however opened new opportunities. These approaches are based on the development of sensor nodes able to operate underground at a few dozens of centimeters deep; see Figure 2 and (Moiroux-Arvis et al., 2022), [START_REF] Saeed | Towards the Internet of Underground Things: a systematic survey[END_REF], [START_REF] Salam | Current advances in Internet of Underground Things[END_REF]. Different measurements can be carried out, as the soil moisture and temperature. The data are emitted by the radio modules, also buried underground, using Low-Power Wide Area Networks (LPWAN) as LoRa technology [START_REF] Augustin | A study of LoRa: long range & low power networks for the Internet of Things[END_REF]. A device located aboveground, e.g. a gateway or a data collector embedded on an Unmanned Aerial Vehicle (UAV), can communicate with the buried sensor nodes and retrieve the measurements, see [START_REF] Cariou | Data collection from buried sensor nodes by means of an Unmanned aerial vehicle[END_REF], [START_REF] Popescu | Advanced UAV-WSN System for Intelligent Monitoring in Precision Agriculture[END_REF]. The advantage to bury and dissimulate all the components underground, including the radio modules and the antennas, is to protect the system against potential damages due to farming activities but also adverse weather conditions, animals, theft and vandalism. The main current applications of WUSNs are smart irrigation, landslide detection, and monitoring of underground infrastructure, see [START_REF] Akyildiz | Wireless underground sensor networks: Research challenges[END_REF], [START_REF] Ferreira | UnderApp: a system for remote monitoring of landslides based on wireless underground sensor networks[END_REF] and [START_REF] Silva | CPS)2 Integration of center pivot systems with wireless underground sensor networks for autonomous precision agriculture[END_REF].

The development of buried sensor nodes is, however, a challenging task as the electromagnetic radio waves are highly attenuated in the soil, about 20-300 times worse than in the air [START_REF] Da Silva | The future of wireless underground sensing networks considering physical layer aspects[END_REF]. This attenuation is, moreover, dependent on numerous environmental factors (e.g. volumetric water content (VMC), burial depth, soil composition). The communication ranges are therefore considerably reduced and highly variable. In practice, the Underground to Underground (UG2UG) communication link between two buried sensor nodes is limited to a few meters. However, as highlighted in our last work [START_REF] Moiroux-Arvis | Evaluation of LoRa technology in 433-MHz and 868-MHz for underground to aboveground data transmission[END_REF], the Underground to Aboveground (UG2AG) communication link can reach more than 250 meters under certain conditions, i.e. 868-MHz radio modules with LoRa communication, transmit power of +14 dBm/25 mW and a burial depth about 15 cm.

Based on these first results, the objective of this paper is to demonstrate that a data collector node embedded on a mobile robot is able, in real time, to retrieve the soil moisture of a buried sensor node when passing at proximity and adapt the working speed of the robot accordingly. This is the objective of the project named "CIDEA" which is funded by the CAP 20-25 International Research Centre -Innovative Transportation and Production Systems (CIR-ITPS). The principle scheme is presented in Figure 3. II. METHODOLOGY The methodology consists first to bury underground a sensor node at a depth d 1 below the surface, see Figure 3. This sensor node regularly measures the soil moisture and emits a data frame using its LoRa radio module.

A collector node is embedded on the mobile robot and connected to its CAN bus. It has reading access to various information on the CAN bus (e.g. GPS position, status of the implement, robot's speed) and can write a CAN message to change the speed of the robot if required.

When the robot is inside the coverage range of the buried sensor node (a 1 in Figure 4), the collector node receives and decodes the data frames which contain the soil moisture information. It will take the control of the speed of the robot with respect to the RSSI (Received Signal Strength Indication) signal level: if the strength of the RSSI signal is above a threshold (-65 dBm), the soil moisture information is taken into account and the speed of the robot is adapted according to the curve presented in Figure 5. When the robot leaves the delock area (-80 dBm), the collector node gives back the control of the speed to the controller of the robot. A Truebner SMT100 probe is connected to measure the soil moisture. The sensor node is powered with a rechargeable lithium battery 3.7 V/8.8 Ah. A temporal reference DS1378 is implemented to switch between the deep-sleep and active modes. The node is activated during the measurement of soil moisture (the measurements are stored on a 8 bit memory EEPROM with I2C-bus interface) and during predetermined time-windows to communicate the data from the radio module. When a time-window is active, the sensor node puts its radio module in listening mode and waits for a message from the collector node. When the communication is established, the sensor node transfers the data to the collector node.

B. The collector node

The collector node is built around a microcontroller Atmega328-AU (Microchip) running at 8 MHz, see Figure 7. It is powered through a rechargeable lithium battery 3.7V / 8.8Ah. It is equipped with a LoRa radio module RFM95W (HopeRF) to communicate with the buried sensor node, and a CAN bus interface to exchange information with the robot (read and write access). All the data of the experiments are stored in a SD card in the collector node enabling to analyze the results afterwards.

C. The mobile robot

The mobile robot used for the tests is presented in Figure 9. This robot is electrically propelled (Lithium-ion batteries 12.8 kWh). The wheelbase is 1.40 m for a total length of 2.50 m. The track is adjustable from 0.70 m to 1.30 m. It weights about 550 kg with a maximal payload of 150 kg. It has four independent wheeled-motors and four independent steering actuators. It is able to reach speeds up to 10 m/s in the wide track configuration. The implement at the rear is a horizontal sprayer. Four nozzles are placed at different working distances, they can be activated independently.

The architecture of this robot is built around three main controllers (steering, speed, implement) joined by two CAN bus. The collector node is connected to the CAN bus enabling to take the control of the speed of the robot.

A. The experimental setup

The sensor node is buried at the depth d 1 =15 cm. Four soil moistures were tested: 0%, 5%, 10% and 20%. At the beginning of each experiment, the robot is manually controlled with its joystick (speed and steering). When the collector node detects the buried sensor node with a RSSI value higher than -65 dBm (i.e. very good reception, that enables to limit the size of the experimental area), it controls the speed of the robot in place of the joystick (from 0 to 1 m/s, see Figure 10). When the RSSI signal becomes lower than -80 dBm, the collector node gives back the control of the speed to the joystick.

Figure 10. Speed of the robot with respect to the soil moisture

B. Results with 0% soil moisture

For the first experiment, the probe is not connected leading to a measured soil moisture of 0%. According to Figure 10, that leads to a setpoint for the speed of 1 m/s. Figure 11 presents the trajectory of the robot and Figure 12 the RSSI level measured by the collector node. The time during which the collector node controls the speed of the robot is depicted in red. The speeds of the robot (setpoints and measurement) are presented in Figure 13. At the beginning (t=0s), the buried sensor node is detected but with a RSSI value equals to -80 dBm (i.e. not taken into account), see Figure 12. The robot is then moved and controlled with the joystick. At t=40s, the RSSI signal becomes upper than -65 dBm. The robot is then controlled by the collector node until the time t=65s (i.e. the RSSI value becomes lower than -80 dBm). This is also the case between t=110s and t=122s. On the top of Figure 13, we clearly observe the commands of the collector node (v=1 m/s), the commands of the speed of the robot from the joystick and the collector node (middle of Figure 13), and the effective speed of the robot (bottom of Figure 13).

C. Results with 5%, 10% and 20% soil moisture The probe is then connected to the buried sensor node. Successively, three experiments are carried out with a dry soil (5% humidity), a slightly wet soil (10% humidity) and a wet soil (20% humidity). The results are presented in the Figures hereafter. For the three experiments, it can be clearly observed that when the buried sensor node is detected with a RSSI value higher than -65 dBm, the collector controls the speed of the robot. In Figure 16, it can be observed for example that the speed of the robot was about 1.5 m/s at the beginning of the experiment. When the robot enters in the speed control area with the buried sensor node, its speed is constrained to be reduced at 0.7 m/s. In Figure 19, the speed is reduced to 0.5 m/s, and to 0.3 m/s in Figure 22. V. CONCLUSION The objective of this paper was to demonstrate the possibility to adapt the behavior of a mobile robot to the variations of the soil moisture encountered in a field. To access information on soil status as close as possible to the intervention area of the mobile robot, we proposed to develop a system based on sensor nodes buried at a few dozens of centimeters deep. The main challenge is however to develop these underground objects with relevant communication ranges. Constraints on energy consumption and transmit power also have to be considered as the nodes have to work ideally several years without battery replacement.

We developed a set of buried sensor nodes able to measure the soil moisture and transmit that information to an aboveground collector node. The radio communication is based on the technology LoRa at 868-MHz. The collector node is embedded in a mobile robot and connected to its bus CAN. It has read and write access, in particular to control the speed of the robot when required. The first experimental results are presented. When the robot is in the communication range of a buried sensor node with a RSSI value upper than a threshold, the moisture information is collected and used to adapt the speed of the robot accordingly.

VI. DISCUSSION

The first results obtained in this paper have considered a single sensor node located underground. In the future work, we intend to study several sensor nodes positioned in a field. That requires to carefully study and manage the overlapping areas. In particular from the knowledge of the position of the robot and the locations of the different moisture sensors in the field, the setpoint for the robot's speed could be determined with precision. Obviously, the speed curve with respect to the soil moisture could also be adapted to the situation. The exponential shape of Figure 8 was only used to rapidly reduce the speed of the robot with respect to the soil moisture. A more in-depth study on the evolution of the speed of the robot according to the humidity of the soil could be carried out. Moreover, the system can take advantage of the circulation of the mobile robot in the field to collect the soil moisture data. These measurements are local and accurate, and are particularly useful to improve the management of the crops (e.g. irrigation planning, spatial distribution of the water).

Funding:

The work presented in this paper is part of the project called "CIDEA". This research is financed by the French government IDEX-ISITE initiative 16-IDEX-0001 (CAP 20-25) -Innovative Transportation and Production Systems (CIR-ITPS). It aims to demonstrate a proof-ofconcept applied to an agricultural context. The opportunity to bury and dissimulate a panel of sensors underground also opens new perspectives in other activity sectors.

Figure 2 .

 2 Figure 2. Buried sensor node measuring the moisture of the soil and emitting the data through radio communication

Figure 3 .

 3 Figure 3. The robot adapts its working speed with respect to the soil moisture information collected from a buried sensor node

Figure 4 .

 4 Figure 4. Speed control area of the collector node

Figure 6 .

 6 Figure 6. The buried sensor node

Figure 7 .

 7 Figure 7. The collector node on the robot

Figure 8 .

 8 Figure 8. Diagram of the program in the collector node

Figure 9 .

 9 Figure 9. Experimental robot used for the tests

Figure 11 .

 11 Figure 11. Trajectory of the robot -Test 0%

Figure 14 .

 14 Figure 14. Trajectory of the robot -Test 5%

Figure 17 .

 17 Figure 17. Trajectory of the robot -Test 10%


David Wellnitz

Guillermo Preisser

Vincenzo Alba

Jerome Dubail

Johannes Schachenmayer
email: schachenmayer@unistra.fr

Rise and fall, and slow rise again, of operator entanglement under dephasing

The operator space entanglement entropy, or simply 'operator entanglement' (OE), is an indicator of the complexity of quantum operators and of their approximability by Matrix Product Operators (MPO). We study the OE of the density matrix of 1D many-body models undergoing dissipative evolution. It is expected that, after an initial linear growth reminiscent of unitary quench dynamics, the OE should be suppressed by dissipative processes as the system evolves to a simple stationary state. Surprisingly, we find that this scenario breaks down for one of the most fundamental dissipative mechanisms: dephasing. Under dephasing, after the initial 'rise and fall' the OE can rise again, increasing logarithmically at long times. Using a combination of MPO simulations for chains of infinite length and analytical arguments valid for strong dephasing, we demonstrate that this growth is inherent to a U (1) conservation law. We argue that in an XXZ spin-model and a Bose-Hubbard model the OE grows universally as 1 4 log 2 t at long times, and as 1 2 log 2 t for a Fermi-Hubbard model. We trace this behavior back to anomalous classical diffusion processes.

The study of quantum many-body systems through the prism of their quantum entanglement continues to prove extremely fruitful [START_REF] Amico | Entanglement in many-body systems[END_REF][START_REF] Eisert | Colloquium: Area laws for the entanglement entropy[END_REF]. In particular, the growth of entanglement in time-evolving quantum many-body systems is of fundamental interest [START_REF] Calabrese | Evolution of entanglement entropy in one-dimensional systems[END_REF][START_REF] Fagotti | Evolution of entanglement entropy following a quantum quench: Analytic results for the xy chain in a transverse magnetic field[END_REF][START_REF] Žnidarič | Many-body localization in the Heisenberg XXZ magnet in a random field[END_REF][START_REF] Alba | Entanglement and thermodynamics after a quantum quench in integrable systems[END_REF][START_REF] Jonay | Coarsegrained dynamics of operator and state entanglement[END_REF][START_REF] Lukin | Probing entanglement in a many-body-localized system[END_REF]: not only is it useful to characterize the dynamics, but the amount of entanglement also indicates whether a quantum evolution can be efficiently simulated on a classical computer. In one dimension (1D) the connection can be made via the concept of matrix product states (MPS) [START_REF] Vidal | Efficient Simulation of One-Dimensional Quantum Many-Body Systems[END_REF][START_REF] Verstraete | Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems[END_REF][START_REF] Schollwöck | The density-matrix renormalization group in the age of matrix product states[END_REF][START_REF] Paeckel | Time-evolution methods for matrix-product states[END_REF]. An MPS is a decomposition of a many-body state vector into a product of χ × χ matrices (where the entries of the matrices are local kets). In such a representation, the bipartite von Neumann entanglement entropy S is bounded by max[S] = log 2 (χ). Consequently, to represent a physical state |ψ(t) with entanglement entropy S(t) as an MPS, the matrix size (or 'bond dimension') has to grow at least as χ ∝ 2 S(t) with time. For example, an evolution where S increases linear in time can therefore be considered computationally hard [START_REF] Schuch | Entropy Scaling and Simulability by Matrix Product States[END_REF].

The past few years have seen the arrival of novel experiments capable of synthetically engineering quantum many-body models in controllable and clean environments, e.g. using optically trapped ultracold atoms [START_REF] Bloch | Quantum simulations with ultracold quantum gases[END_REF][START_REF] Adams | Rydberg atom quantum technologies[END_REF][START_REF] Browaeys | Many-body physics with individually controlled Rydberg atoms[END_REF][START_REF] Morgado | Quantum simulation and computing with Rydberg-interacting qubits[END_REF], molecules [START_REF] Gadway | Strongly interacting ultracold polar molecules[END_REF], or ions [START_REF] Blatt | Quantum simulations with trapped ions[END_REF]. Since such experiments are currently bringing the goal of analog quantum simulation into sight [START_REF] Cirac | Goals and opportunities in quantum simulation[END_REF][START_REF] Georgescu | Quantum simulation[END_REF], the question of entanglement growth, and thus classical simulability, has become very important.

Every experiment has small couplings to its environment, and should therefore be considered as an open quantum system described by a density matrix ρ. Analogously to MPSs for pure states, also a matrix product operator (MPO) form of the density matrix ρ can be defined. An MPO form allows to easily express the density matrix as Schmidt decomposition between a left and a right block:

ρ = a λ a τ [L] a τ [R] a , (1)
with Tr(τ

[L/R] a τ [L/R] b
) = δ ab , and λ a the Schmidt coefficients [schematically this is depicted in Fig. 1(b)]. The bipartite entropy of this decomposition is given by the 'operator space entanglement entropy' or simply 'operator entanglement' (OE) defined as [START_REF] Jonay | Coarsegrained dynamics of operator and state entanglement[END_REF][START_REF] Zanardi | Entangling power of quantum evolutions[END_REF][START_REF] Zanardi | Entanglement of quantum evolutions[END_REF][START_REF] Wang | Quantum entanglement of unitary operators on bipartite systems[END_REF][START_REF] Prosen | Operator space entanglement entropy in a transverse Ising chain[END_REF][START_REF] Dubail | Entanglement scaling of operators: a conformal field theory approach, with a glimpse of simulability of long-time dynamics in 1+1d[END_REF][START_REF] Zhou | Operator entanglement entropy of the time evolution operator in chaotic systems[END_REF][START_REF] Alba | Operator entanglement in interacting integrable quantum systems: The case of the rule 54 chain[END_REF][START_REF] Wang | Barrier from chaos: operator entanglement dynamics of the reduced density matrix[END_REF][START_REF] Styliaris | Information Scrambling over Bipartitions: Equilibration, Entropy Production, and Typicality[END_REF]]

S OP = - a λ 2 a log 2 λ 2 a . (2)
The OE quantifies how many Schmidt values are at least needed for faithfully approximating decomposition (1), and thus indicates the efficiency of an MPO representation [31]. It can be easily computed numerically [START_REF] Prosen | Operator space entanglement entropy in a transverse Ising chain[END_REF][START_REF] Zhou | Operator entanglement entropy of the time evolution operator in chaotic systems[END_REF][START_REF] Alba | Operator entanglement in interacting integrable quantum systems: The case of the rule 54 chain[END_REF][START_REF] Noh | Efficient classical simulation of noisy random quantum circuits in one dimension[END_REF][START_REF] Rakovszky | Dissipation-assisted operator evolution method for capturing hydrodynamic transport[END_REF] and is amenable to analytical treatment [START_REF] Dubail | Entanglement scaling of operators: a conformal field theory approach, with a glimpse of simulability of long-time dynamics in 1+1d[END_REF][START_REF] Bertini | Operator Entanglement in Local Quantum Circuits I: Chaotic Dual-Unitary Circuits[END_REF][START_REF] Bertini | Operator Entanglement in Local Quantum Circuits II: Solitons in Chains of Qubits[END_REF]. We stress that OE is not necessarily connected to genuine quantum entanglement between distinct blocks of spins when ρ is a mixed state. Still, it is a crucial quantity as it puts severe restrictions on the possibility to approximate ρ by an MPO. Furthermore, OE can give insights into quantum many body effects such as quantum chaos and information scrambling [START_REF] Zhou | Operator entanglement entropy of the time evolution operator in chaotic systems[END_REF][START_REF] Wang | Barrier from chaos: operator entanglement dynamics of the reduced density matrix[END_REF][START_REF] Styliaris | Information Scrambling over Bipartitions: Equilibration, Entropy Production, and Typicality[END_REF].

Here, we analyze the far-from-equilibrium dynamics of the OE, S OP , in open many-body quantum systems. Our models include coherent nearest-neighbor Hamiltonian couplings that compete with incoherent single-particle dephasing [START_REF] Rossini | Coherent and dissipative dynamics at quantum phase transitions[END_REF][START_REF] Cai | Algebraic versus Exponential Decoherence in Dissipative Many-Particle Systems[END_REF][START_REF] Medvedyeva | Exact Bethe Ansatz Spectrum of a Tight-Binding Chain with Dephasing Noise[END_REF][START_REF] Foss-Feig | Solvable Family of Driven-Dissipative Many-Body Systems[END_REF][START_REF] Žnidarič | Relaxation times of dissipative many-body quantum systems[END_REF] at rate γ [see sketch in Fig. 1(a)]. Under dephasing, fluctuations and coherences can decay towards equilibrium in a universal algebraic and subdiffusive manor [START_REF] Poletti | Interaction-Induced Impeding of Decoherence and Anomalous Diffusion[END_REF][START_REF] Poletti | Emergence of Glasslike Dynamics for Dissipative and Strongly Interacting Bosons[END_REF][START_REF] Ren | Noise-Driven Universal Dynamics towards an Infinite Temperature State[END_REF][START_REF] Bouganne | Anomalous decay of coherence in a dissipative many-body system[END_REF]. We treat dissipation in a Lindblad master equation. Dephasing arises due to a coupling with the environment in which local magnetization is preserved, e.g. for laser driven transitions due to laserphase fluctuations [START_REF] Plankensteiner | Laser noise imposed limitations of ensemble quantum metrology[END_REF][START_REF] Gardiner | Quantum noise[END_REF], or due to spontaneous photoabsorption of lattice photons in optical lattices [START_REF] Pichler | Nonequilibrium dynamics of bosonic atoms in optical lattices: Decoherence of many-body states due to spontaneous emission[END_REF]. We compute the evolution of S OP for an infinite MPO representation of the full density matrix using an infinite time evolving block decimation algorithm (iTEBD) with re-orthogonalization [START_REF] Orús | Infinite time-evolving block decimation algorithm beyond unitary evolution[END_REF]. Surprisingly, we find that for the magnetization conserving XXZ model and well-defined initial magnetization [see text below Eq. ( 3) for the definition of our models], the OE exhibits a universal logarithmic growth at long times [see Fig. 1(c/d)]. An identical universal behavior is also observed for particle number conserving Bose-and Fermi-Hubbard models [see discussion below]. Strikingly, as shown in Fig. 1(e), this logarithmic growth breaks down if the symmetry is broken by initial states, or for Liouvillians without magnetization conservation (see Appendix for more examples). In the latter scenarios S OP saturates, or even vanishes at long times. In this paper we explain this behavior by considering symmetryresolved OE in combination with known results for classical models of interacting particles [START_REF] Harris | Diffusion with "Collisions" between Particles[END_REF][START_REF] Levitt | Dynamics of a Single-File Pore: Non-Fickian Behavior[END_REF][START_REF] Arratia | The Motion of a Tagged Particle in the Simple Symmetric Exclusion System on Z[END_REF].

Model -We focus on infinite spin-1/2 chains, which evolve under the general Hamiltonian ( ≡ 1):

Ĥ = 1 4 i J x σx i σx i+1 +J y σy i σy i+1 +J z σz i σz i+1 + h z 2 i σz i . (3)
Here, σx,y,z i denote standard Pauli matrices defined in a local basis |↓, ↑ i , J x,y,z are the respective nearestneighbor spin couplings, and h z is a field strength along the z direction. Our Hamiltonian (3) includes: i) the XXZ model, with J x = J y ≡ J, h z = 0; ii) an XYZ model, with J x ≡ J, J x = J y , h z = 0; and iii) a transverse Ising model, with J x ≡ J, J y = J z = 0, h z = 0. We are interested in the dynamics of highly-excited states. Here, we choose pure Néel product states polarized along the z direction, ρ0 = |ψ 0 ψ 0 | with |ψ 0 = i |↑ 2i-1 |↓ 2i , or a tilted Néel state along the x direction,

|ψ 0 = i |→ 2i-1 |← 2i with | i = (|↑ i ∓ |↓ i )/ √ 2.
Dynamics is governed by a Lindblad master equation,

d dt ρ = -i[ Ĥ, ρ] + i D [i] ρ ≡ Lρ , (4)
with the local dephasing super-operators

D [i] ρ = γ/2(σ z i ρσ z i -ρ
) and the Liouvillian super-operator L. MPO decomposition and OE -For L spins, the full many-body density matrix ρ of a spin-1/2 system is a 2 L × 2 L hermitian matrix with unit trace. The amount of information encoded in ρ can be effectively compressed using matrix product decompositions [START_REF] Schollwöck | The density-matrix renormalization group in the age of matrix product states[END_REF][START_REF] Paeckel | Time-evolution methods for matrix-product states[END_REF]. This can be done in different ways [START_REF] Weimer | Simulation methods for open quantum many-body systems[END_REF]: For instance, by decomposing ρ into a particular (not unique) statistical mixture of pure states, while using an MPS for the latter. Then, the Lindblad dynamics can be computed using quantum trajectories [START_REF] Daley | Quantum trajectories and open many-body quantum systems[END_REF]. Alternatively, one can use a direct MPO representation for ρ, e.g. simply by effectively vectorizing local density matrices [START_REF] Zwolak | Mixed-State Dynamics in One-Dimensional Quantum Lattice Systems: A Time-Dependent Superoperator Renormalization Algorithm[END_REF], or by constructing MPOs in a locally purified form which preserves positivity [START_REF] Verstraete | Matrix Product Density Operators: Simulation of Finite-Temperature and Dissipative Systems[END_REF][START_REF] Werner | Positive Tensor Network Approach for Simulating Open Quantum Many-Body Systems[END_REF].

Here, we decompose ρ into a canonical MPO form [START_REF] Zwolak | Mixed-State Dynamics in One-Dimensional Quantum Lattice Systems: A Time-Dependent Superoperator Renormalization Algorithm[END_REF], which is formally achieved by an iterative application of the Schmidt decomposition from Eq. ( 1), until each spin n is described by a matrix of unique local operators γ [n] an,an+1 , i.e. ρ = {an} n λ an γ[n] an,an+1 . By choosing local basis operators, êin , for the density matrix of spin n, we can then decompose:

ρ = {in} χ {an} n Γ [n] in anan+1 λ [n]
an n êin ,

where Γ [n] are three-dimensional tensors and λ [n] the Schmidt vectors. Tensors are truncated at a maximum MPO bond dimension χ. All results shown are numerically converged in χ (see Appendix). For êi we choose local eigenoperators for the multiplication with σz from the left and right, to take advantage of the magnetization conservation (see below). Both the initial state and the Hamiltonian are invariant under translation by two lattice sites. As a consequence, in Eq. ( 5) one has [START_REF] Orús | Infinite time-evolving block decimation algorithm beyond unitary evolution[END_REF][START_REF] Vidal | Classical Simulation of Infinite-Size Quantum Lattice Systems in One Spatial Dimension[END_REF], and only two Γ and λ tensors are needed to encode a density matrix. The time evolution is then computed with a fourth order Trotter decomposition of the matrix exponential of the Liouvillian exp(L∆t) [START_REF] Orús | Infinite time-evolving block decimation algorithm beyond unitary evolution[END_REF][START_REF] Zwolak | Mixed-State Dynamics in One-Dimensional Quantum Lattice Systems: A Time-Dependent Superoperator Renormalization Algorithm[END_REF][START_REF] Sornborger | Higher-order methods for simulations on quantum computers[END_REF]. Importantly, since the dynamics is non-unitary, a naïve implementation of this algorithm destroys the orthogonality of the decomposition in Eq. ( 5), such that with time, the λs do not correspond to orthogonal Schmidt bases anymore. We fix this by re-orthogonalizing the tensors after updates [START_REF] Orús | Infinite time-evolving block decimation algorithm beyond unitary evolution[END_REF].

Γ [n+2] = Γ [n] and λ [n+2] = λ [n]
When considering the XXZ model, the total magnetization Ŝz = n σz n is conserved. This means that ρ stays an eigenoperator of Ŝz in the sense that Ŝz ρ(t) = M ρ(t) at all times (for the Néel state, M = 0). Note that alternatively, one could also define a condition for multiplication from the right. Due to magnetization preservation, the τ [R] a matrices in Eq. ( 1) can be chosen to be eigenoperators of the 'right-half magnetization' of the chain, Ŝz R ≡ n>0 σz n (w.l.o.g. we define the right half as

n > 0), Ŝz R τ [R] a = M R τ [R] a . Similarly, one can choose τ [L]
a to be eigenoperators of Ŝz L ≡ n≤0 σz n with M L = -M R . This means that the index a in Eq. ( 1) becomes a composite index a → (M R , a ), where a distinguishes the Schmidt coefficients corresponding to the same M R :

ρ = M R √ p M R a λ M R ,a τ [L] -M R ,a τ [R] M R ,a . (6)
Here we defined

λ M R ,a ≡ λ M R ,a / √ p M R , with p M R = a λ 2 M R
,a the probability of having magnetization M R in the right half. The existence of the conservation law makes our simulations much more efficient, since the block-diagonal form of the tensors can be exploited.

Logarithmic increase of OE: Numerical results -In simulations in Fig. 1 we noticed a distinctive different behavior of OE growth at long times (log-growth) for the magnetization conserving XXZ model compared to other models breaking this conservation law. Quite generically, at times t γ -1 , the dynamics is dominated by the Hamiltonian part in Eq. ( 4). Sufficiently pure states at such short times can be approximated by the state |ψ t = e -i Ĥt |ψ 0 . In that case the OE is simply twice the entanglement entropy of |ψ t (see e.g. [START_REF] Dubail | Entanglement scaling of operators: a conformal field theory approach, with a glimpse of simulability of long-time dynamics in 1+1d[END_REF]), and it is well established that the latter generically grows linearly in time (in the absence of disorder). At times t γ -1 , the initial coherence is destroyed by dephasing, and the OE decreases (see Appendix for a more detailed discussion on the parameter dependence of the peak-heights). This 'rise and fall' is clearly visible in Fig. 1, and it is typical for OE dynamics, and also other quantities such as the mutual information [START_REF] Carollo | Emergent dissipative quasiparticle picture in noninteracting Markovian open quantum systems[END_REF][START_REF] Alba | Hydrodynamics of quantum entropies in Ising chains with linear dissipation[END_REF]. Typically, under the dynamics in Eq. (4) the system is expected to relax to a simple stationary state characterized by the conserved quantities of Eq. ( 4), or to the identity if there is no conservation law. The OE at late times converges towards the OE of that stationary state. This is visible in our simulations of the XYZ and Ising models, see Fig. 1(e). In this case only the parity Π = i σz i is preserved by the dynamics. Since the initial Néel state is an eigenstate of Π, the stationary density matrix is a projector on a fixed parity sector, 1 2 (1 ± Π), with the O(1) entropy S OP = 1 (= log 2 2). For the XXZ chain with the initial tilted Néel along the x direction, even parity conservation is broken, and the stationary density matrix becomes the identity, S OP = 0. In stark contrast, for the XXZ chain and initial Néel state, after the rise and fall dynamics, the OE increases again at long times, see Fig. 1(c), and this second increase is logarithmic in time, see Fig. 1(d). More precisely, we find the long-time behavior S OP (t → ∞) = η log 2 (tJ) + S 0 , which we will also understand analytically below. The prefactor η converges universally to η = 1/4 independent on the precise values of γ and J z , as shown in Fig. 2(a), and has also been observed with additional disorder [START_REF] Medvedyeva | Influence of dephasing on many-body localization[END_REF]. The offset S 0 depends on the characteristic time-scale of the long-time diffusive dynamics set by J, J z and γ (see Appendix). Note that we find the evolution of OE in the XXZ model to be independent of the signs of J and J z .

Mechanism for logarithmic growth: Abelian symmetry and anomalous charge diffusion -To also analytically understand this behavior, we now consider the XXZ model evolution in the strong dephasing limit γ J. The dissipators in the master equation (4) project the density matrix onto its diagonal part ρdiag = σ ρ σσ |σ σ|, where the σ denote all binary vectors of spin-z configurations. The dynamics then reduces to a classical master equation for the probability p σ = ρ σσ , dp σ /dt = σ M σσ p σ . The stochastic matrix M was determined in Ref. 37 in second-order perturbation theory starting from Eq. ( 4). It takes the form of an effective ferromagnetic Heisenberg Hamiltonian, M = -J 2 /(8γ) i σx

i σx i+1 + σy i σy i+1 + σz i σz i+1 -1 .
Importantly, M is the stochastic matrix of the Symmetric Simple Exclusion Process (SEP) [START_REF] Mallick | The exclusion process: A paradigm for nonequilibrium behaviour[END_REF][START_REF] Bernard | Transport in quantum chains under strong monitoring[END_REF], a model of classical hard-core particles that perform random walks.

Crucially, in the SEP, the mean squared displacement of a tagged particle grows as X 2 t ∝ √ t [50, 51, 64-67], as opposed to ∝ t for a usual random walk. This anomalous diffusion is universally found in problems of so-called 'single-file diffusion' [START_REF] Lin | From Random Walk to Single-File Diffusion[END_REF][START_REF] Hahn | Single-File Diffusion Observation[END_REF][START_REF] Wei | Single-File Diffusion of Colloids in One-Dimensional Channels[END_REF], when classical particles diffuse in one-dimensional channels without bypassing each other. Here, it is now also tied to anomalous scaling of the particle number fluctuations between the left and right half-systems. If ∆N (t) = M R (t)/2 is the excess number of particles (w.r.t. the initial Néel state) in the right half-system at time t, and if we tag the particle initially at the origin, then one can esti-mate ∆N (t) ρ 0 X t , where ρ 0 = 1/2 is the particle density in the Néel state. Consequently, ∆N (t) 2 ρ 2 0 X 2 t ∝ √ t. More generally, the probability distribution of M R (t) is found to obey a scaling form at long times [START_REF] Derrida | Current Fluctuations of the One Dimensional Symmetric Simple Exclusion Process with Step Initial Condition[END_REF]:

p(M R (t) = m) ∼ t→∞ exp √ t G m √ t
. Here, the large-deviation function G is non-positive, symmetric [G(u) = G(-u)], diverges when |u| → ∞, and has a single minimum at u = 0 (see [START_REF] Derrida | Current Fluctuations of the One Dimensional Symmetric Simple Exclusion Process with Step Initial Condition[END_REF]). In particular, away from the tails the distribution is Gaussian with standard deviation δ = t 1/4 / |G (0)|. The Shannon entropy associated to number fluctuations is then

S num (t) = m∈2Z -p(m) log 2 p(m) - 2 πδ 2 e -m 2 2δ 2 log 2 2 πδ 2 e -m 2 2δ 2 dm 2 = log 2 δ + log 2 πe/2 = t→∞ 1 4 log 2 t + O(1). ( 7
)
It is no coincidence that S OP grows in the same way as the number fluctuations S num at long times (see below).

Away from the γ/J 1 limit, the XXZ chain no longer reduces to the SEP. Nevertheless we find that the same type of anomalous scaling persists. This is confirmed in Fig. 2, where we show that, for times accessible numerically, p(M R = m) is approximately Gaussian [Fig. 2(b)] with width δ ∝ t 1 4 [Fig. 2(c)]. Thus, even though the exact correspondence with the SEP breaks down at finite γ/J, the scaling of S num in Eq. ( 7) remains unchanged. This result is also consistent with previous studies of transport in the XXZ and related models [START_REF] Žnidarič | Dephasing-induced diffusive transport in the anisotropic Heisenberg model[END_REF][START_REF] Žnidarič | Exact solution for a diffusive nonequilibrium steady state of an open quantum[END_REF][START_REF] Eisler | Crossover between ballistic and diffusive transport: the quantum exclusion process[END_REF][START_REF] De Nardis | Subdiffusive hydrodynamics of nearlyintegrable anisotropic spin chains[END_REF].

Symmetry-resolved OE -We now show how the relation between S OP and S num can be revealed in so-called symmetry resolved operator entanglement. From Eq. ( 6) we can derive a decomposition of the OE into the form:

S OP = M R p M R S res (M R ) + S num (p M R ), (8)
where the 'symmetry-resolved entanglement entropies' are

S res (M R ) = -a λ 2 M R ,a log 2 ( λ 2 M R ,a
), and S num is given in Eq. ( 7). Such symmetry-resolved entropies have attracted attention recently [START_REF] Lukin | Probing entanglement in a many-body-localized system[END_REF][START_REF] Goldstein | Symmetry-resolved entanglement in many-body systems[END_REF][START_REF] Xavier | Equipartition of the entanglement entropy[END_REF][START_REF] Parez | Quasiparticle dynamics of symmetry-resolved entanglement after a quench: Examples of conformal field theories and free fermions[END_REF][START_REF] Barghathi | Rényi Generalization of the Accessible Entanglement Entropy[END_REF][START_REF] Barghathi | Operationally accessible entanglement of one-dimensional spinless fermions[END_REF]. In Fig. 2(d) we display S res (M R ) for different values of M R . Also S res exhibits the 'rise and fall' phenomenon discussed above, but independent of M R they decrease to very small values at late times. This means that the logarithmic increase S OP is solely due to the growth of S num .

Fermi-and Bose-Hubbard model -To demonstrate the generality of the logarithmic OE growth, we discuss two additional paradigmatic many-body setups featuring number conservation: i) a Fermi-Hubbard (FH) model, ĤFH = -J f n,σ (ĉ † σ,n ĉσ,n+1 + h.c.) + In both cases we consider dephasing D [k]

U f n ĉ † ↑,n ĉ † ↓,
ρ = γ Lk ρ L † k -γ( L † k Lk ρ+ ρ L † k Lk )/2
, where Lσ,n = ĉ † σ,n ĉσ,n and Ln = b † n bn in the FH and BH case, respectively. As demonstrated in Fig. 3, both models also exhibit a long-time logarithmic OE growth. In the FH model, we observe S OP ∼ log 2 (tJ f )/2. For γ J f , U f , this can be understood analytically by considering the FH chain as the sum of two chains, one for each spin degree of freedom, both of which are described by the SEP and exhibit S OP ∼ log 2 (tJ f )/4. Here, interactions contribute only at higher orders (see Appendix). The BH model exhibits S OP ∼ log 2 (tJ b )/4 analogous to the XXZ model. For large U b γ J b the creation of doublons is energetically suppressed, leading back to the SEP. For finite interaction strength (γ J b , U b ), a different classical limit is reached, which also features logarithmic OE growth with a prefactor close to 1/4 stemming from a 'Symmetric Inclusion Process' (see Appendix and Ref. [START_REF] Bernard | Transport in quantum chains under strong monitoring[END_REF]).

Conclusion -

We showed that in a dissipative system possessing a U (1) conservation law the operator entanglement grows logarithmically at long times. We pinpointed the mechanism that leads to this logarithmic growth, and identified its prefactor with the [possibly anomalous] exponent characterizing the fluctuations of the charge associated with the U (1) symmetry. Our results and methods are of general interest to studies of imperfect quantum computation and quantum simulation platforms, currently pushing into a regime where they may offer a quantum advantage. The entanglement entropy dynamics we study here connects directly to similar results obtained for discrete quantum circuit models [START_REF] Zhou | Operator entanglement entropy of the time evolution operator in chaotic systems[END_REF][START_REF] Noh | Efficient classical simulation of noisy random quantum circuits in one dimension[END_REF][START_REF] Li | Quantum Zeno effect and the many-body entanglement transition[END_REF][START_REF] Chan | Unitary-projective entanglement dynamics[END_REF], or to other dissipative simulation methods such as quantum trajectories [START_REF] Coppola | Growth of entanglement entropy under local projective measurements[END_REF][START_REF] Botzung | Engineered dissipation induced entanglement transition in quantum spin chains: From logarithmic growth to area law[END_REF]. An understanding of the destructive processes of the environment on dynamics are essential, and the interplay between dissipation and coherent couplings can lead to interesting physics or state engineering (e.g. [START_REF] Mark | Interplay between coherent and dissipative dynamics of bosonic doublons in an optical lattice[END_REF][START_REF] Zhu | Suppressing the Loss of Ultracold Molecules Via the Continuous Quantum Zeno Effect[END_REF][START_REF] Shchesnovich | Control of a Bose-Einstein condensate by dissipation: Nonlinear Zeno effect[END_REF][START_REF] Müller | Measurement-Induced Dark State Phase Transitions in Long-Ranged Fermion Systems[END_REF]). In the future, it will be interesting to investigate how the presence of more complex symmetries such as SU(N) impacts entanglement dynamics. Whenever two ∆t values are given, different ∆ts were used for different parameters for historic reasons, but are clearly converged in either case. σ+ i+1 . In the limit of strong dissipation, we can then consider perturbation theory for small L 1 .

The eigenstates of L 0 are product states in the z eigenbasis. Of these eigenstates, only those without off-diagonal correlations do not decay. Perturbative coupling between the states can be computed by [see also Eq. (S6) in Sec. E]

L eff = L 0 -PL 1 (L 0 ) -1 L 1 P , (S1)
where P projects into the L 0 eigenstates with zero decay rate. According to Ref. [START_REF] Cai | Algebraic versus Exponential Decoherence in Dissipative Many-Particle Systems[END_REF], the transfer rate between sites i and (i + 1) depends on the state of the neighboring sites (i -1) and (i + 2). Thus, focussing only on the two central spins, the transfer rates can be read from the prefactors of

L i,(i+1) eff |• • • ↑↑↓↑ • • • • • • ↑↑↓↑ • • •| = J 2 4γ |• • • ↑↓↑↑ • • • • • • ↑↓↑↑ • • •| -|• • • ↑↑↓↑ • • • • • • ↑↑↓↑ • • •| , (S2) L i,(i+1) eff |• • • ↑↑↓↓ • • • • • • ↑↑↓↓ • • •| = J 2 γ 4γ 2 + J 2 z |• • • ↑↓↑↓ • • • • • • ↑↓↑↓ • • •| -|• • • ↑↑↓↓ • • • • • • ↑↑↓↓ • • •| . (S3)
Then for τ 1 we can replace e τ L0 with P. This leads to

+ τ γ 2 PL 1 P ⊥ -L ⊥ 0 -1 P ⊥ L 1 P |ρ 0 + . . . ,
where L ⊥ 0 is the restriction of L 0 to the subspace orthogonal to its kernel, and L ⊥ 0 -1 is its inverse on that subspace.

Finally, one notes that PL 1 P = 0, since L 1 acts on the density matrix ρ as -i[H, ρ]. Thus, the result of second-order perturbation theory is that, in the strong dephasing limit, the density matrix remains diagonal on long times scales, |ρ(t) P |ρ(t) . The slow evolution of the diagonal part |ρ diag (t) := P |ρ(t) is given by

d dt |ρ diag (t) = - 1 γ PL 1 P ⊥ L ⊥ 0 -1 P ⊥ L 1 P |ρ diag (t) . ( S6
)
This equation has been obtained for the XXZ chain by Cai and Barthel [START_REF] Cai | Algebraic versus Exponential Decoherence in Dissipative Many-Particle Systems[END_REF] and also for other models, including models of bosons, by Bernard, Jin and Shpielberg [START_REF] Bernard | Transport in quantum chains under strong monitoring[END_REF]. Let us now apply it to the Fermi-Hubbard model (S1).

The key observation is that the interacting term U N j=1

ĉ † ↑j ĉ↑j ĉ † ↓j ĉ↓j in the Hamiltonian (S1) acts diagonally in the computational basis, so it does not contribute to the term -1 γ PL 1 P ⊥ L ⊥ 0 -1 P ⊥ L 1 P in Eq. (S6). Therefore, in the strong dephasing limit, the spin components ↑ and ↓ simply decouple, and each spin component follows its own strong-dephasing dynamics. One can then simply set U = 0 in the Hamiltonian (S1) which splits into two independent models of non-interacting fermions, one for each spin component. Consequently, the strong-dephasing limit of the Fermi-Hubbard model consists of two decoupled Symmetric Exclusion Processes, one for the spins ↑, the other for the spins ↓, as claimed in the main text. and the Lindlblad equation with dephasing takes the form (see e.g. Ref. [START_REF] Pichler | Nonequilibrium dynamics of bosonic atoms in optical lattices: Decoherence of many-body states due to spontaneous emission[END_REF])

d dt ρ = -i[ Ĥ, ρ] + γ N j=1 D[ b † j bj ](ρ). (S8)
Like for the XXZ chain and for the Fermi-Hubbard model, the strong dephasing limit can be analyzed in second-order perturbation theory [START_REF] Cai | Algebraic versus Exponential Decoherence in Dissipative Many-Particle Systems[END_REF][START_REF] Bernard | Transport in quantum chains under strong monitoring[END_REF]. The general formula (S6), where γL 0 is the dissipative part of the Lindblad equation and L 1 is the unitary part, applies also to the Bose-Hubbard model. As in the Fermi-Hubbard case, the interaction term

U 2 N j=1 b † j bj ( b † j bj -1)
in the Hamiltonian (S7) acts diagonally in the single-site Fock basis. Consequently, it does not contribute to the term -1 γ PL 1 P ⊥ L ⊥ 0 -1 P ⊥ L 1 P. At strong dephasing the Hubbard interaction is irrelevant, and the model can be analyzed simply by setting U = 0. The resulting model of non-interacting bosons with dephasing has been analyzed in Ref. [START_REF] Bernard | Transport in quantum chains under strong monitoring[END_REF]. For the convenience of the reader, here we briefly sketch a simple derivation of that strong-dephasing limit.

For notational simplicity, we focus on N = 2 two sites only, and label the basis states as |n 1 n 2 where n 1 , n 2 ∈ N are the number of bosons on site 1 and 2. Diagonal density matrices are of the form

ρ diag = n1,n2 p(n 1 , n 2 ) |n 1 n 2 n 1 n 2 | . (S9)
Under strong dephasing, the effective dynamics within the space of diagonal density matrices must be a classical stochastic process for the probability distribution p(n 1 , n 2 ). The goal is to determine the rates that define the classical master equation of that process.

We proceed step by step to compute the r.h.s of Eq. (S6). First,

L 1 P |ρ diag = L 1 |ρ diag is the vectorized form of -i Ĥ, ρ diag = -i H, n1,n2 p(n 1 , n 2 ) |n 1 n 2 n 1 n 2 | = iJ n1,n2 p(n 1 , n 2 ) (n 1 + 1)n 2 |n 1 + 1, n 2 -1 n 1 n 2 | + n 1 (n 2 + 1) |n 1 -1, n 2 + 1 n 1 n 2 | -(n 1 + 1)n 2 |n 1 , n 2 n 1 + 1n 2 -1| + n 1 (n 2 + 1) |n 1 , n 2 n 1 -1, n 2 + 1| ,
which we write as

L 1 |ρ diag = iJ n1,n2 p(n 1 , n 2 ) (n 1 + 1)n 2 |n 1 + 1, n 2 -1 n 1 n 2 | + n 1 (n 2 + 1) |n 1 -1, n 2 + 1 |n 1 n 2 -(n 1 + 1)n 2 |n 1 , n 2 n 1 + 1n 2 -1| + n 1 (n 2 + 1) |n 1 , n 2 |n 1 -1, n 2 + 1 . Second, we observe that L[ b † 1 b1 ] |n 1 ± 1, n 2 |n 1 , n 2 = (n 1 ± 1)n 1 - (n1±1) 2 +n 2 1 2 = -1 2 . Then (L ⊥ 0 ) -1 P ⊥ L 1 P |ρ diag = -iJ 2γ n1,n2 p(n 1 , n 2 ) (n 1 + 1)n 2 |n 1 + 1, n 2 -1 |n 1 n 2 + n 1 (n 2 + 1) |n 1 -1, n 2 + 1 |n 1 n 2 -(n 1 + 1)n 2 |n 1 , n 2 |n 1 + 1n 2 -1 + n 1 (n 2 + 1) |n 1 , n 2 |n 1 -1, n 2 + 1 .
Applying again the vectorized form of -i[ Ĥ, .], and projecting onto diagonal configurations, one arrives at

-PL 1 P ⊥ (L ⊥ 0 ) -1 P ⊥ L 1 P |ρ diag = J 2 2γ n1,n2 p(n 1 , n 2 ) (n 1 + 1)n 2 2 (|n 1 + 1, n 2 -1 |n 1 + 1n 2 -1 -|n 1 , n 2 |n 1 n 2 ) + n 1 (n 2 + 1) 2 (|n 1 -1, n 2 + 1 |n 1 -1n 2 + 1 -|n 1 , n 2 |n 1 n 2 ) + (n 1 + 1)n 2 2 (|n 1 + 1, n 2 -1 |n 1 + 1n 2 -1 -|n 1 , n 2 |n 1 n 2 ) + n 1 (n 2 + 1) 2 (|n 1 -1, n 2 + 1 |n 1 -1, n 2 + 1 -|n 1 , n 2 |n 1 , n 2 ) .
In conclusion, the classical master equation obeyed by the probability distribution p(n Interestingly, the resulting classical dynamics exhibits a signature of bosonic bunching in the original quantum model. The classical stochastic process defined by Eq. (S11) belongs to a class of models known as 'Symmetric Inclusion Processes' in the statistical physics literature, see e.g. Refs. [START_REF] Giardina | Correlation inequalities for interacting particle systems with duality[END_REF][START_REF] Grosskinsky | Dynamics of condensation in the symmetric inclusion process[END_REF]. In order to determine the dynamics of number fluctuations in the model described by Eq. (S11), we numerically compute n samples = 5000 different sample trajectories for a chain of length L = 100. Fig. S4 shows the time trace of the variance of the number of particles in the left half of the chain var(N L ) [corrected by n samples /(n samples -1)]. We compute the error of the variance using the jackknife method [START_REF] Miller | The jackknife-a review[END_REF] by separating the data into M = 500 bins with 10 trajectories each. Then, we compute the variances v m for all possible subsets of M -1 different bins (i.e. deleting 1 bin for each subset). Finally, we compute the error from these variances as

error[var(N L )] = M -1 M m [v m -var(N L )] 2 . (S12)
We use LsqFit.jl to fit the numerical data with a power law of the form var(N L ) = a × (t/τ ) b , and find a ≈ 0.49 and b = 0.53, close to b = 1/2. We find that this fit describes the data well, as shown in Fig. S4. Here, τ = γ/J 2 is the characteristic time scale of the dynamics.

If we thus assume that N L is Gaussian distributed with standard deviation δ ∼ (t/τ ) b/2 , we can use Eq. ( 7) of the main text to compute S num (t) = log 2 [(t/τ ) b/2 ] + O(1), close to log 2 (t/τ )/4 + const. used in the main text.

FIG. 1 .

 1 FIG. 1. (a) We compute dynamics of spin chains with coherent nearest-neighbor couplings (double-arrows) and local dephasing at rate γ (wiggly arrows). (b) We analyze the growth of OE (SOP) for a bipartition of an infinite chain into a left/right block [from a Schmidt decomposition of the density matrix ρ, see Eq. (1)]. (c) Time evolution of SOP for a Néel state in the XXZ model for different values of γ = J/4, J/2, J, 2J, 4J in order of increasing darkness (anisotropy Jz = -J). (d) Same as (c), demonstrating logarithmic growth at long times (log-scale time axis). Grey dashed: Analytic long-time prediction: SOP = log 2 (Jt)/4 + const. (e) SOP for dynamics in models breaking magnetization conservation. Green dash-dotted: XYZ model (Jx = J, Jy = 0.8J, Jz = -J/2, γ = J/2), Red dashed: transverse field Ising model (hz = J, γ = J/2), blue solid: XXZ model with initial Néel state in x-direction (Jz = -J/2, γ = J/2). Results converged for time step ∆tJ = 0.2 [∆tJ = 0.5 for panel (d) at long times] and different values for the bond dimension χ = 256, 512, 1024 (see Appendix for convergence).

FIG. 2 .

 2 FIG. 2. (a) Numerical determination of η for long times. η at time t0J is obtained as the local tangent to SOP(t) = S0 +η log(tJ) at t0J. We find η → 1/4 (grey dashed line) for all parameters and t0 → ∞. Top (blue): Fixed Jz = -J and various γ/J = 1/4, 1/2, 1, 2, 4. Bottom (red): Fixed γ = J and various -Jz/J = 0, 1, 2, 3, 4. (b) Probabilities p(MR) for right-half magnetization MR of the infinite chain (see text) at increasingly late times (200 ≤ tJ ≤ 1200 from light to dark, Jz = -J, γ = J/2. Lines are Gaussian fits. (c) Variance δ 2 of Gaussian fits [as in (b)] as a function of time for Jz = -J and γ/J = 1/4, 1/2, 1, 2, 4 (light to dark). The grey dashed line indicates δ 2 ∼ √ tJ (double-log scale). (d) Symmetry-resolved operator entanglement entropies Sres as a function of time (see text). For short times and the larger MR = 6, 8 the probabilities p(MR) are exponentially suppressed and no sub-machine-precision data could be extracted. Same parameters as panel (b). Results converged for ∆tJ = 0.5 [(a)-(c)] and ∆tJ = 0.1 [(d)] and different values of χ = 512, 1024, 2048 (see Appendix).

FIG. 3 .

 3 FIG. 3. Logarithmic OE growth in the Fermi-Hubbard (a) and Bose-Hubbard (b) models with dephasing for different interaction strengths U f /b . The grey dashed lines indicate the analytically expected long-time growth (see Appendix). The initial states are: One fermion per site with alternating spins |• • • ↑↓↑ • • • (a) and alternating sites with 0 or 1 boson |• • • 101 • • • (b). Parameters: γ f = 8J f (a) and γ b = 2J b (b), χ = 512, 256, ∆tJ f /b = 1/2, maximum bosons per site nmax = 4.

  n ĉ↑,n ĉ↓,n with creation operators for spinfull fermions on site n, ĉ † σ,n (σ =↑, ↓); and ii) a Bose-Hubbard (BH) model, ĤBH = -J b n ( b † n bn+1 + h.c.) + U b /2 n b † n b † n bn bn , for bosons created by b † n .

1 Figure

 1 FIG. S2. Rise and fall of operator entanglement without magnetization conservation. (a) We consider the XXZ model with dephasing, as in the main text, for Jz = -J, γ = J/2. The initial state is a fully polarized state n exp[iθσ (x) n ] |↓ n , which are not and eigenstate of Ŝz (except for θ = 0). (b) We consider the XXZ model with Jz = -J for a Néel initial state as in the main text. Instead of dephasing, we consider non-conserving dissipation D ± n ρ = γ/2(2σ ± n ρσ ∓ n -σ∓ n σ± n ρ -ρσ ∓ n σ± n ), and for various rates γ. Simulations for χ = 128, ∆tJ = 0.1

2 .Ĥ

 2 Strong dephasing limit of the Bose-Hubbard model: a Symmetric 'Inclusion' ProcessWe define the Bose-Hubbard Hamiltonian as

FIG. S4 .

 S4 FIG. S4. Fit for the variance growth for the Symmetric Inclusion Process. The blue crosses are determined by a Monte Carlo simulation of Eq. (S11) with 5000 trajectories. The grey dashed line is a power law fit to the numerical data.

  (τ -τ 2 )PL 1 P ⊥ e τ2L0 P ⊥ L 1 P |ρ 0 + . . . dτ 2 e τ2L0 P ⊥ L 1 P |ρ 0 + . . .

	|ρ(τ )	P |ρ 0 + dτ 2 exp τ γ PL 1 P |ρ 0 + τ 2 2γ 2 PL 1 PL 1 P |ρ 0 + τ 1 γ 2 0 τ γ ∞ PL 1 P P |ρ 0 + PL 1 P ⊥ τ γ 2 0
		exp	τ γ	PL 1 P P |ρ 0

  1 , n 2 ) is (n 2 + 1)[p(n 1 -1, n 2 + 1) -p(n 1 , n 2 )] + 1)n 2 [p(n 1 + 1, n 2 -1) -p(n 1 , n 2 )]. (S10)This result generalizes straightforwardly to the case with more sites, N > 2, i.e. [(n i-1 + 1)p(. . . , n i-1 + 1, n i -1, . . . ) + (n i+1 + 1)p(. . . , n i -1, n i+1 + 1, . . . ) -(n i-1 + n i+1 + 2)p(. . . , n i , n i+1 , . . . )]. (S11)What is interesting about the master equation (S11) is that it is very different from the one of non-interacting particles undergoing a random walk, which would rather be of the form

		d dt (n 1 d p(n 1 , n 2 ) = J 2 γ n 1 + J 2 γ dt J 2 p(n 1 , n 2 , . . . , n N ) = γ
	d dt	p(n

i n i 1 , n 2 , . . . , n N ) ∝ i n i [p(. . . , n i -1, n i+1 + 1, . . . ) + p(. . . , n i-1 + 1, n i -1, . . . ) -2p(. . . , n i-1 , n i , n i+1 , . . . )].

Acknowledgments We thank P. Calabrese, O. Castro-Alvaredo, A. Grabsch, G. Misguich, S. Murciano, G. Pupillo and G. Schütz for helpful discussions. This work was supported by LabEx NIE under contract ANR-11-LABX0058 NIE, and the QUSTEC program, which has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sk lodowska-Curie grant agreement number 847471. This work is part of the Interdisciplinary Thematic Institute QMat, as part of the ITI 2021-2028 pro-gram of the University of Strasbourg, CNRS and Inserm, and was supported by IdEx Unistra (ANR-10-IDEX-0002), SFRI STRAT'US project (ANR-20-SFRI-0012), and EUR QMAT ANR-17-EURE-0024 under the framework of the French Investments for the Future Program.

Our MPO codes make use of the intelligent tensor library (ITensor) [87]. 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 Alsacalcul/Big Data.

APPENDICES

In Sec. A, we provide details on the numerical convergence. In Sec. B, we discuss a breaking of magnetization conservation on different levels. In Sec. C we show numerical results on the parameter dependence of the short time peaks of the OE. In Sec. D, we discuss the perturbative estimate for the offset of the operator entanglement. In Sec. E, we identify the strong dephasing limits of the Fermi-Hubbard model, and of the Bose-Hubbard model, which are instrumental for explaining the prefactors of the logarithmic growth of the OE observed in the main text.

Appendix A: Details on numerical convergence For solving the time-evolution with the iTEBD algorithm [START_REF] Orús | Infinite time-evolving block decimation algorithm beyond unitary evolution[END_REF][START_REF] Vidal | Classical Simulation of Infinite-Size Quantum Lattice Systems in One Spatial Dimension[END_REF] we make use of a Trotter decomposition of the matrix exponential of the super-operator governing the full dissipative dynamics [START_REF] Zwolak | Mixed-State Dynamics in One-Dimensional Quantum Lattice Systems: A Time-Dependent Superoperator Renormalization Algorithm[END_REF]. In order to reach long timescales, we implement 4-th order decomposition [START_REF] Sornborger | Higher-order methods for simulations on quantum computers[END_REF], which allowed for time-step converged results even up to step sizes of ∆t = 1/(2J) [see Fig. S1(a)]. ∆t signifies the full time step, which is composed of individual gates with time-steps of length ∆t/12 and ∆t/6. In particular, we use the method

Throughout the paper we only show simulations up to times until which the results are converged in the MPO bond dimension χ, defined in the truncated MPO decomposition in Eq. ( 5). In practice, we ran simulations repeatedly, doubling bond dimensions in the different runs until lines become visually indistinguishable. We varied the bond dimension in the range χ = 128, 256, 512, 1024, 2048. Note that e.g. for χ = 1024, a maximum possible entropy of S OP ≤ log 2 (χ) = 10 is theoretically supported. Naturally, the evolution producing the largest OE required the largest values of χ. The convergence plots in Fig. S1(b) and Fig. S1(c) demonstrate this procedure for "worst-case" scenarios (i.e. the data with the largest OE values and for long times). We choose the data from Fig. 1(d) with γ = J/4, simulated up to long times, and the data for Fig. 3(a) with J z /J = -1 and γ = J/8, which reached the largest values in the short time peak S OP ∼ 6.

In particular, the following parameters have been used for the figures in the main text: In this section, we give a few more examples without magnetization conservation. First, we discuss fully polarized states that are tilted at an angle θ, n exp[iθσ

n ] |↓ n . For 0 < θ < π, these states are not eigenstates of Ŝz . We find that after the initial growth, the OE decays to zero, just like for the initial tilted Néel state. For the special cases θ = 0(π), the state

) is an eigenstate of the Liouvillian, and thus the OE remains zero at all times.

If magnetization conservation is broken on the level of the Liouvillian, we find a similar rise and fall of OE, which remains zero at long times [see Fig. S2(b)]. For the case of both an incoherent pump and and decay at rate γ, after the initial rise and fall, we expect the steady state to be the (trivial) infinite temperature state with S OP = 0. Whether a state can be efficiently approximated by an MPO is determined by the OE of that state, which has a local maximum at short times. Thus, until not too long times, the approximability of the dynamics is determined by the value at this maximum S peak . This value is shown in Fig. S3(a) as a function of the dephasing rate γ and the interaction strength J z . We find that generally the peak-height grows as S peak ∼ 1/γ for small γ. Interestingly, S peak decreases when increasing the spin interaction strength |J z |. As a consequence, the dynamics under strong dephasing can always be simulated, while for small dephasing the dynamics can only be simulated in the presence of strong interactions (large |J z |).

Appendix D: OE offset scaling

Whether a state can be efficiently approximated by an MPO is determined by the OE of that state, which is at long times dominated by the logarithmic growth S OP = log 2 (tJ)/4 + S 0 . Since logarithmic growth is very slow, the offset is also important. We find that the offset S 0 decreases with increasing |J z | or γ [see Fig. S3(b)] when extrapolating to t 0 → ∞. The approximately equal spacing of lines in a regime of large γ (when doubling γ) indicates a scaling of S 0 ∼ log 2 (1/γ) in this regime. This is confirmed by the perturbative argument given below, which leads to S 0 = log 2 (J/2 4γ 2 + J 2 z )/4 + O(1). In the analytical model, the offsets are given by the time scale of the transfer, as the multiplicative transfer rate under the logarithm becomes an additive constant (see below). Cai and Barthel have computed the transfer rates perturbatively for small J γ, as discussed in the following [START_REF] Cai | Algebraic versus Exponential Decoherence in Dissipative Many-Particle Systems[END_REF]. We start by writing the evolution operator as a sum

Here, the Hamiltonian is separated into interaction Ĥ0 = J z /4 i σz

All other possibilities can be obtained by permuting bras and kets, up and down, and reading states from right to left.

Here, L i,(i+1) eff

indicates that only the effective Liouvillian acting on site i and i + 1 is given, i.e.

. Since both scenarios are combinatorically equally likely for randomly arranged spins, we compute an effective transfer rate r by taking the geometric average, which yields

From this, we can compute the parameter scaling of S 0 using S(t) = log 2 (rt)/4 + O(1) = log 2 (Jt)/4 + log 2 (r/J)/4 + O(1) as We define the Fermi-Hubbard Hamiltonian as

where ĉ † ↑j (resp. ĉ † ↓j ) is the creation operator of a fermion with spin up (resp. down) on site j. We consider the strong dephasing limit of the Lindblad equation

where D[ Â](ρ) = Âρ Â † -1 2 { Â † Â, ρ}. For our purposes, it is more convenient to write this equation in vectorized form: we vectorize the density matrix ρ → |ρ -we write the 4 N × 4 N matrix ρ as a 4 2N -vector |ρ -, and write its evolution as

where L 0 contains the dissipative part of Eq. (S2), and L 1 contains the unitary part in (S2). [Here, contrary to what is done in Ref. [START_REF] Cai | Algebraic versus Exponential Decoherence in Dissipative Many-Particle Systems[END_REF] and in Sec. D, we do not treat the diagonal and non-diagonal parts of the Hamiltonian differently.

Here Importantly, the subspace of density matrices ρ that satisfy L 0 |ρ = 0 is precisely the one of matrices that are diagonal in the computational basis. Let us call P the projector onto that subspace, and P ⊥ = 1 -P its orthogonal projector.