\section*{Method}
\noindent\textbf{Galaxy cluster sample selection.} 
We construct our sample of high-$z$ galaxy clusters using the latest {\it Planck} Sunyaev-Zel'dovich (SZ) catalog, i.e. PSZ2 \cite{planckcoll16}. This catalog provides a reliable cluster mass estimation for galaxy clusters up to $z\sim1$. Although the exact effect of selection biases is still being debated and may overestimate the observed fraction of clusters with radio halos, 
there is no clear consensus in the literature that suggests the SZ sample chosen here disproportionately favours merging over relaxed clusters \cite{eckert+11,rossetti+17,andrade-santos+2017}. We select all the objects at $\rm DEC\geq 20^\circ$, to match the part of the sky with the best LOFAR sensitivity, and $z\geq0.6$. No mass-threshold has been applied in our cluster selection. 
These selection criteria result in a sample of 30 objects (see Supplementary Table 1). The sample includes clusters where precise redshift measurements became available recently \cite{amodeo+18,barrena+18,burenin+18,sereno+18,streblyanska+18,vanderburg+16,zohren+19} and a previous-discovered {\it Planck} cluster, i.e. PLCK\,G147.3--16.6 \cite{vanweeren+14}. We used the optical images from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS \cite{panstarss16}) to refine the accuracy of the cluster center coordinates, using the brightest cluster galaxy (BCG), when detected, or approximately the center of the distribution of the galaxies. We also inspected the {\it Chandra} and/or {\it XMM-Newton} X-ray image, when available, to verify the correspondence of the cluster coordinates with the center of the thermal gas distribution. Among these 30 {\it Planck} clusters, 21 objects were already covered as part of the LOFAR Two-metre Sky Survey (LoTSS \cite{shimwell+19}). From these 21 clusters one was affected by bad ionospheric observing conditions (i.e. PSZ2\,G088.98+55.07) and one lays on the same line of sight of a cluster at $z=0.3$ (i.e. PSZ2\,G097.52+51.70). Those are therefore excluded from the final sample (see Extended Data Figure 1), consisting of 19 objects at median redshift of $z=0.7$.

\bigskip
\bigskip

\noindent\textbf{LOFAR observations, data reduction, images and flux measurements.} 
The LOFAR observations were carried out together with the LOFAR Two-metre Sky Survey (LoTSS \cite{shimwell+19}) in the 120--168 MHz frequency range. The survey consists in 8 hours of observation for each pointing, with a field of view of $\sim2.6^\circ$, full resolution of about $6''$, and median noise of about $70~\mu$Jy beam$^{-1}$. Given the large sky coverage already achieved by the survey, our targets were often observed by more than one pointing (see Supplementary Table 2), which further improved our signal to noise. For each LoTSS pointing, we performed standard data reduction \cite{shimwell+19}, which includes direction-independent and dependent calibration, and imaging of the full LOFAR field of view using \texttt{prefactor} \cite{vanweeren+16,williams+16,degasperin+19}, \texttt{killMS} \cite{tasse14,smirnov+tasse15} and \texttt{DDFacet} \cite{tasse+18}. We additionally improved the quality of the calibration, using the products of the pipeline, subtracting all the sources outside a region of $\sim15^\prime\times15^\prime$ surrounding the target, and performing extra phase and amplitude self-calibration loops in each target sub-field. To verify the quality of the calibration, each pointing is imaged independently using \texttt{WSClean} \cite{offringa+14}, with the wideband deconvolution mode (\texttt{channelsout=6}). For all our imaging, we employ automatic clean masks for the deconvolution. We utilize a $3\sigma_{\rm rms}$ masking threshold, with $\sigma_{\rm rms}$ the local map noise, and cleaning down to the $1\sigma_{\rm rms}$ level inside the mask. The images have a central frequency of 144 MHz, and are provided in the Extended Data Figure 1. 

\bigskip

\noindent All the clusters in our sample were carefully visually inspected in the full-resolution image to search for extended radio emission. 
To emphasize the presence of diffuse radio emission, we also produced low-resolution source-subtracted images. We first applied an $uv$-cut to the data, to filter out emission associated to linear sizes larger than 500 kiloparsecat the cluster redshift and 
to create a clean component model of the compact sources. For PSZ2\,G089.39+69.36 we apply an $uv$-cut corresponding to 400 kpc.
During this step, we employed multiscale deconvolution \cite{offringa+17}, using scales of $\rm [0,4,8,16]\times pixelscale$ (with the pixel size of $1.5''$) to properly include and subtract somewhat extended radio galaxies. In addition, for the automatic deconvolution we lower the mask threshold to $1\sigma_{\rm rms}$ to subtract the faintest contaminating sources (see Supplementary Figure 1). 
Finally, we subtracted the compact sources model from the visibilities 
and tapered the $uv$-data to $\sim15''$ resolution. It is important to note that in case of an extended radio galaxy, with linear size $\gtrsim 500$~kpc, this method cannot properly subtract the radio emission from the $uv$-data. For that reason, the extended double-lobed radio galaxy just north of PLCKG147.3--16.6 has been manually excluded and blanked in the low-resolution source-subtracted image (Figure 1 and Extended Data Figure 1). 
\bigskip

\noindent We used these low-resolution source-subtracted images to define the largest linear size (LLS) of the diffuse radio emission, following the $2\sigma_{\rm rms}$ contour. Systems with $\rm LLS<500~kpc$ were not included in our statistics.
Roundish, centrally-located significantly detected sources (see the next paragraph) that follow the ICM distribution were classified as radio halos; elongated, peripheral structures were classified as radio relics. We did not identify clear examples of fossil radio plasma sources (i.e. ``radio phoenices''), these typically have small sizes (i.e. 200~kpc or smaller), irregular shapes, and are not centered on the cluster center \cite{degasperin+17,mandal+20}. Similarly, we did not find clear examples of  ``mini-halos'', as these diffuse sources usually have LLS smaller than 500 kpc and are located in dynamically relaxed systems \cite{giacintucci+17}. Uncertain classifications,
due to contamination from extended radio galaxies in the cluster region and/or  ambiguous shapes,
were excluded from our final statistical analysis. 
If we measure the LLS following the $3\sigma_{\rm rms}$ contour, we would classify 8 radio halos instead of 9 (see Supplementary Table 3). 


\bigskip

\noindent To determine the integrated flux densities for the radio halos, we first measured the total cluster radio flux densities (i.e. compact sources and diffuse emission) from the low-resolution image, following the $2\sigma_{\rm rms}$ radio contour to fully cover the extension of the diffuse radio emission (see Extended Data Figure 1). In the Supplementary Table 3, we also report the integrated flux density measurements following the $3\sigma_{\rm rms}$ radio contour. %
We then mathematically subtracted the flux densities of the compact sources within the cluster region, measured in the full-resolution $uv$-cut image, to obtain the flux density associated to the diffuse radio sources \cite{cassano+19}. 
For PSZ2\,G091.83+21.16, which also hosts a bright radio relic, we visually separated it from the radio halo region (see Supplementary Figure 2). 
The uncertainties on the flux densities are estimated adding in quadrature the $\sigma_{\rm rms}$ statistical uncertainty, the 15\% systematic error associated to the LOFAR flux scale calibration, $f_{\rm cal}$ \cite{shimwell+19}, and the uncertainty due to the source subtraction in the cluster region, $\sigma_{\rm sub}$ \cite{cassano+13}, according the following Equation:  

\begin{equation}
\Delta S = \sqrt{(f_{\rm cal} S_{\rm 144MHz})^2 + \sigma_{\rm rms}^2 N_{\rm beam,h} + \sigma_{\rm sub}^2} \, .
\end{equation}
Here, $\sigma_{\rm sub}^2=\sum_i\sigma_{\rm rms}^2 N_{{\rm beam,s}_i}$, and $N_{\rm beam, h}$ and  $N_{{\rm beam, s}_i}$ are the number of beams covering the halo and the subtracted $i$ sources, respectively. Uncertainties associated to possibly missed faint compact sources are not included, but these should be smaller than the systematic error associated to the LOFAR flux scale calibration.

\bigskip

\noindent We confirmed the presence of diffuse radio emission when the measured integrated radio flux density is larger than five times the flux uncertainty, i.e. $\frac{S_{\rm 144MHz}}{\Delta S}\geq 5$ (see Column 5 in Table 1).
Flux density measurements directly obtained from the low-resolution source-subtracted images are also provided in the Supplementary.


\bigskip

\noindent Finally, we calculated the $k$-corrected radio luminosity at 1.4 GHz:

\begin{equation}
P_{\rm 1.4GHz} =  4\pi D_L^2 \dfrac{S_{\rm 1.4GHz}}{(1+z)^{\alpha+1}} \, ,
\end{equation}
with $D_L$ the luminosity distance assuming standard $\Lambda$CDM cosmology (i.e. $H_{0} = 70$ km s$^{-1}$~Mpc$^{-1}$, $\Omega_{m} = 0.3$, and $\Omega_{\Lambda} = 0.7$) and $S_{\rm 1.4GHz}=S_{\rm 144MHz}\left (\frac{\rm 1.4GHz}{\rm 144MHz} \right )^\alpha$. Since the spectral indices of the radio halos are not measured, we assume a Gaussian distribution of
spectral indices $\alpha$ with mean value $-1.5$ and scatter $0.3$ for computing the radio luminosity. These values were chosen to encompass the typical range of spectral indices found in radio halos \cite{vanweeren+19}, and to take into account that they could be steeper than local systems \cite{cassano+06}. Only for PLCK\,G147.3--16.6, we used the spectral index obtained combining our observations with the GMRT literature data \cite{vanweeren+14} (see Supplementary). The radio luminosities are listed in Table 1,
whose uncertainties are obtained by performing Monte Carlo simulations including the flux density uncertainties and the spectral index uncertainty. 

\bigskip
\bigskip


\noindent\textbf{X-ray observations and data reduction.}
Twelve clusters among our LOFAR sample have also {\it Chandra} and/or {\it XMM-Newton} observations (see Supplementary Table 4 and Extended Data Figure 2).


\bigskip

\noindent {\it Chandra} data were reduced using the \texttt{chav} software package (\url{http://hea-www.harvard.edu/\textasciitilde alexey/CHAV/}) with \texttt{CIAO v4.6} and applying the \texttt{CALDB v4.7.6} calibration files \cite{vikhlinin+05}. The processing includes the application of gain maps to calibrate photon energies, filtering out counts with ASCA grade 1, 5, or 7 and bad pixels, and a correction for the position-dependent charge transfer inefficiency (CTI). Periods with count rates with a factor of 1.2 above and 0.8 below the mean count rate in the 6--12 keV band were also removed. Standard blank-sky files were used for background subtraction. The final exposure-corrected images were made in the 0.5--2.0~keV band using a pixel binning of a factor of 4 (i.e $2^{\prime\prime}$) and combining the different ObsIDs together where relevant (i.e., PSZ2\,G160.83+81.66). 


\bigskip

\noindent {\it XMM-Newton} observations were reduced using version 17.0.0 of the dedicated Science Analysis System (SAS). After converting the observation data files to unfiltered event lists, we extracted light curves in bins of 100 seconds in the energy ranges 10--12~keV for the MOS and 10--14~keV for the pn detectors onboard. Good time intervals were created excluding periods when the mean count rates in these light curves were different from the corresponding mean by more than 2 standard deviations. Images were extracted in the 0.4--7.0 keV band, using only the highest quality (FLAG==0) single to quadruple MOS events (PATTERN$\leq$12) and single to double pn events (PATTERN$\leq$4). Images were weighted by the respective exposure maps for each detector and combined, accounting for a factor $\sim$2.5 difference between the typical expected count rates for MOS and pn at a given source flux (due to the slightly different effective areas and the fact that the light focused by two of the XMM-Newton telescopes is split evenly between the MOS and RGS detectors). Since we are interested mostly in the morphology of high-redshift, thus relatively compact, clusters positioned close to the aim point, the effects of vignetting and instrumental particle background are minimal and have not been corrected for this analysis.

\bigskip
\bigskip


\noindent\textbf{Turbulent energy flux at high redshift.}
According to re-acceleration models, relativistic electrons in radio halos are re-accelerated by turbulence. Assuming that turbulence in the ICM is mainly injected by cluster mergers, then the turbulence injection rate is $\rho v_i^3/L_{\rm inj}$ \cite{casano+brunetti05}, with $v_i$ the impact velocity between the two merging clusters, $\rho$ the cluster mean density and $L_{\rm inj}$ the injection scale of turbulence (which can be assumed to be of the same order of the cluster size, about 1 Mpc at both $z=0.2$ and $z=0.7$).
The relative impact velocity of two subclusters with virial masses $M_v$ and $\Delta M_v$ which collide (at a distance $R_{v}$ between the centers) starting from an initial distance $d_0$ with zero velocity is given by \cite{sarazin02}:

\begin{equation}
v_i\simeq \sqrt{2G\frac{(M_{v}+\Delta M_v)}{R_{v}}
\left (1-\frac{1}{\eta_{v}}\right)} \, .
\label{eq:vi}
\end{equation}
Here, $d_0=\eta_v R_{v}$, $\eta_v\simeq 4\big(\frac{ M_{v}+ \Delta M} {M_{v}}\big)^{1/3}$, and $R_{v}$ is the virial radius of the main cluster, i.e., the radius at which the ratio between the average density in the cluster and the mean cosmic density at the redshift of the cluster is given by \cite{kitayama+suto96}:

\begin{equation}
\Delta_{c}(z)=
18\pi^2(1+0.4093\,\omega(z)^{0.9052}) \, ,
\label{Dc}
\end{equation}
where $\omega(z)\equiv \Omega_{f}(z)^{-1}-1$ with:

\begin{equation}
\Omega_{f}(z)
=\frac{\Omega_{m,0}(1+z)^3}{\Omega_{m,0}(1+z)^3+\Omega_{\Lambda}} \, .
\label{Omz}
\end{equation}



\bigskip

\noindent The virial mass, $M_v$, and the virial radius are thus related by:

\begin{equation}
R_{v}=\left [\frac{3M_{v}}{4\pi\Delta_{c}(z)\rho_{m}(z)}
\right]^{1/3}
\label{Rv}
\end{equation}
where $\rho_{m}(z)=2.78\times10^{11}\,\Omega_{m,0}\,(1+z)^3~h^2$ M$_{\odot}$ Mpc$^{-3}$ is the mean mass density of the Universe at redshift $z$. The ratio of the gas density at high and low redshift can be assumed to scale as the ratio of the two virial densities, i.e., $\frac{\rho_{v,z=0.7}}{\rho_{v,z=0.2}}=\frac{M_{v,z=0.7}}{\left (R_{v,z=0.7} \right)^3}\times \frac {\left (R_{v,z=0.2} \right )^3}{M_{v,z=0.2}}$. In the case $M_{v,z=0.7}=M_{v,z=0.2}$, the density ratio becomes simply proportional to $\left (R_{v,z=0.2}/R_{v,z=0.7} \right )^3$. Specifically, considering both at low and high-$z$ $M_v\simeq 1.0\times 10^{15}$ M$_{\odot}$ and $\Delta M_v\simeq M_{v}/3$, the above equations give a ratio $\left (v_{i,z=0.7}/v_{i,z=0.2} \right)^3\sim 1.5$ and $\rho_{v,z=0.7}/\rho_{v,z=0.2}\sim 2.2$. The resulting ratio of turbulent energy flux derived at high and low redshift is then $\sim 3.3$.


\bigskip

\noindent The above calculation is approximate, because we only took kinematic effects and $\Lambda$CDM cosmology into account. We thus also use simulations in \cite{casano+brunetti05} to calculate the injection of turbulence during simulated cluster merging histories, and derived the turbulent energy injected in the volume swept by the infalling subcluster. We found that in this case the distributions of the turbulent energy fluxes derived for clusters with $M_v\simeq 0.9-1.4\times 10^{15}\,M_{\odot}$ in the two redshift ranges $0.2-0.3$ and $0.7-0.8$ differ by a factor of  $\sim 3-4$, in a good agreement with the simpler derivation reported above.


\bigskip
\bigskip



\noindent\textbf{Radiative lifetime and electrons acceleration in high-$z$ clusters.} 
In the paper we have assumed that $\eta_{\rm rel}$, i.e. the fraction of turbulent energy flux that is absorbed by relativistic electrons (Eq. 1), is constant with redshift. This depends on the energy budget of the population of re-accelerated seed electrons, and on the microphysics of the nonlinear interaction between turbulence and particles. The latter is assumed to be redshift-independent (e.g. \cite{brunetti+lazarian07}). Thus, a constant $\eta_{\rm rel}$ implies a similar energy budget of the seed electrons that are accumulated in high- and low-$z$ ICM.


\bigskip

\noindent More specifically, we used Eq. 1 to infer a lower limit to the ratio of the magnetic field in high- and low-$z$ radio halos, thus our conclusions may change in the case that clusters at higher redshift host a population of seed electrons with a much larger energy budget. In this section we show that this possibility is very unlikely, as the higher injection rate of electrons in the ICM expected from high-redshift sources is balanced by the stronger electron energy losses.


\bigskip

\noindent The budget of seed electrons that accumulate in the ICM can be estimated as $N_e \propto Q_e \times T_{\rm max}$ (with $Q_e$ the injection rate and $T_{\rm max}$ the maximum electron lifetime). The lifetime of relativistic electrons in the ICM, which is determined by the energy losses due to Coulomb interaction with the thermal plasma and to Inverse Compton scattering with the CMB photons and synchrotron losses, can be expressed as \cite{brunetti+jones14}:

\begin{equation}
T=4\left \{\frac{1}{3}\tilde{\gamma}\beta_{z}+\frac{\tilde{n}}{\tilde{\gamma}}\left[1.2+\frac{1}{75}\ln \left (\frac{\tilde{\gamma}}{\tilde{n}}\right )\right]\right\}^{-1} \, ,
\label{eq:Tlife}
\end{equation}
where we put $\tilde{\gamma}=\gamma/300$, $\tilde{n}= n/10^{-3}$, $\beta_z=\left ( \frac{B}{3.2} \right )^2+(1+z)^4$. Here, $\gamma$ is the Lorentz factor, $n$ the number density in cm$^{-3}$ and $B$ the magnetic field strength in $\mu$G. This lifetime increases with the energy $\gamma$ of relativistic electrons at lower energy, it reaches a maximum and then decreases as a function of $\gamma$.


\bigskip

\noindent The maximum of $T$ is obtained for $\tilde{\gamma}\approx\sqrt{\frac{3\tilde{n}}{\beta_z}}$. For $B< B_{\rm CMB}$, it is determined by the combination of IC and Coulomb losses. Replacing this in Eq. \ref{eq:Tlife} gives the typical lifetime of relativistic electrons that can be accumulated in the ICM on cluster lifetime and then re-accelerated during mergers:

\begin{equation}
T_{\rm max}\,[{\rm Gyr}]=4\left \{\frac{\beta_z}{3}\sqrt{\frac{3\tilde{n}}{\beta_z}}+\tilde{n}\sqrt{\frac{\beta_z}{3\tilde{n}}}\right\}^{-1} \, .
\label{eq:Tmax}
\end{equation}
Considering that between $z=0.2$ and $z=0.7$ the virial density (and hence the thermal gas density $n$) increases by a factor of $\sim2$, the maximum lifetime ratio for electrons radiating at these two redshifts is $\sim 2.7$, considering $B \sim 1-2 ~\mu$G.


\bigskip

\noindent Several sources can inject electrons in the ICM including AGN, galaxies (galactic winds, GW) and cosmological shocks.
Most of the energy of the relativistic plasma of GW and shocks is in the form of supra-thermal and relativistic protons with only a negligible fraction in the form of electrons (e.g. \cite{brunetti+jones14}).


\bigskip

\noindent AGN are thought to be the dominant sources of relativistic electrons in the ICM (e.g. \cite{brunetti+jones14}); this hypothesis is also supported by recent LOFAR observations, where sources of fossil radio plasma from old AGN-electrons are commonly found in galaxy clusters \cite{degasperin+17,mandal+20}. The injection rate of electrons from a number of AGN in clusters, $N_{\rm AGN}$, is $Q_e \sim N_{\rm AGN} n_{\rm rel} V / \Theta$, where $n_{\rm rel}$ is the number density of radio-emitting electrons in the jets and lobes, $V$ is the volume of the radio plasma and $\Theta$ is the AGN life-time. The number density of radio-emitting  electrons in the AGN lobes can be estimated assuming equipartition between particles and magnetic fields, i.e. $n_{\rm rel}\propto \sqrt{L_{\rm syn} / V}$, with $L_{\rm syn}$ the AGN synchrotron luminosity. This gives $Q_e \propto \sqrt{L_{\rm syn} V} N_{\rm AGN}/\Theta$. The linear size of radio galaxies is observed to decline with redshift, probably due to the higher density of the medium in which radio lobes and jets expand \cite{neeser+95,blundell+99}. Thus we attempt to account for the expected change of $V$ with redshift assuming approximate pressure equilibrium of the lobes with the surrounding medium. This would imply $P_{\rm ICM} \propto n_{\rm rel} \propto \sqrt{L_{\rm syn} / V}$ and consequently $Q_e \propto L_{\rm syn} N_{\rm AGN} / ( P_{\rm ICM} \Theta ) 
\approx \rho_L (1+z)^3/ ( P_{\rm ICM} \Theta ) $, where $\rho_L$ is the luminosity density of radio AGN per comoving-volume. The luminosity density increases by a factor $\sim 2-2.5$ from $z\sim 0.2$ to $z\sim 0.7$ (\cite{smolciv+17} their Fig. 5). The increase of $P_{\rm ICM}$ with redshift can be estimated using virial quantities, i.e. $P_{\rm ICM} \propto n_{\rm ICM} T_{\rm ICM} \propto \frac{M_v}{R_v^3}\frac{M_v}{R_v}$, implying that pressure increases by a factor $\sim 2.5$ from $z\sim 0.2$ to $z\sim 0.7$ for clusters with the same virial mass.
We find that $Q_e$ is only about 2-3 times larger in clusters at $z\sim0.7$ if we assume that the life-time of AGN does not change much with redshift. Combining this result with the maximum lifetime of the radiating electrons, that is $\sim 2.7$ times longer at lower redshift, suggests that the budget of seed electrons at low and high redshift is similar.


\bigskip

\noindent Seed electrons in the ICM can also be generated by the decay chain of inelastic proton-proton collisions (e.g. \cite{brunetti+lazarian11,pinzke+17}). The importance of this channel depends on the energy budget of CR in the ICM, which is poorly known. A relevant contribution of re-accelerated secondary particles to radio halos is still allowed by current gamma-ray limits from FERMI-LAT (e.g. \cite{brunetti+17}). For this reason, we also compare the budget of secondary electrons accumulated in galaxy clusters at low and high redshift considering cosmological shocks as the main sources of relativistic protons. In this case, numerical simulations suggest the ratio of CR to thermal ICM pressures ($X=P_{\rm CR}/P_{\rm ICM}$) in high-$z$ clusters is generally smaller than that at low redshift (e.g. \cite{vazza+14}). The ratio of the energy budget of secondary electrons at high and low redshift can be approximately estimated as:
\begin{equation} 
\frac{N_{e,z=0.7}}{N_{e,z=0.2}} \sim 
\frac{X_{z=0.7}}{X_{z=0.2}} \frac{ n_{{\rm ICM},z=0.7} P_{{\rm ICM,}z=0.7}}{n_{{\rm ICM},z=0.2} P_{{\rm ICM,}z=0.2}} \frac{T_{{\rm max},z=0.7}}{T_{{\rm max},z=0.2}} \, ,
\label{secondari}
\end{equation}
implying $N_{e,z=0.7}/N_{e,z=0.2}\sim 0.4-1.8$ if we assume typical CR pressure ratios measured in simulations, i.e. $\sim 0.2-1$ (\cite{vazza+12,vazza+14}, Figs. 14 and Fig. 12 respectively), and the ratios between physical quantities as derived above: $T_{{\rm max},z=0.7}/T_{{\rm max},z=0.2}\sim 1/2.7$, $P_{{\rm ICM},z=0.7}/P_{{\rm ICM},z=0.2} \sim 2.5$ and  $n_{{\rm ICM},z=0.7}/n_{{\rm ICM},z=0.2} \sim 2$. 
Similarly to the case of AGN, this implies that the budget of secondary particles that is available at higher redshift is very similar to that at low redshift. 

\bigskip
\bigskip



\noindent\textbf{Reynolds number estimation.}
In Figure 3b
we have used a simple model of magnetic field amplification to infer combined constraints on the initial magnetic field and on the Reynolds number.


\bigskip

\noindent We follow the method as outlined in \cite{beresnyak12}. We assume that the amplification initially operates in a kinematic regime, where the magnetic field grows exponentially with time:

\begin{equation}
B^2(t) \sim B^2_0 \exp \left (t \, \Gamma \right)  \, ,
\label{bexp}
\end{equation}
where the time-scale of the magnetic growth is $\Gamma^{-1} \sim 30 L_d/\delta v_d$. This depends on the turbulent eddy turnover time 
at the viscous dissipation scale, $L_d \sim L_{\rm inj}\, {\rm Re}^{-3/4}$,  $\delta v_d \sim v_t (L_d/L_{\rm inj})^{1/3}$. The factor $\sim 30$ is derived from simulations \cite{cho14,beresnyak+miniati16},  and accounts for the less effective stretching of the field lines due to the turbulent diffusion. When the magnetic and the kinetic energy densities become comparable at the viscous dissipation scale ($B(t)^2/8 \pi \sim 1/2 \rho \delta v_d^2$), the turbulent dynamo transits to a phase  where the magnetic field grows linearly with time. This phase is known as the non-linear regime because the magnetic field (Reynolds stress) becomes dynamically important. From Eq. \ref{bexp}, this transition occurs at a time:

\begin{equation}
t_* \sim 60 \frac{L_{\rm inj}}{v_t} {\rm Re}^{-1/2} \ln  \left ( \frac{B_*}{B_0} \right )
\label{timeexp}
\end{equation}
where $B_* \sim \sqrt{4 \pi \rho}\, v_t {\rm Re}^{-1/4}$.
According to simulations, in the non-linear regime, i.e. for time $t \geq t_*$, we assume that a constant fraction of the turbulent kinetic energy flux, $\eta_B$, is channeled into magnetic field energy ($\eta_B=0.05$ \cite{beresnyak12}). Therefore, the evolution of the magnetic field with time is:


\begin{equation}
B^2(t) \sim B^2_* + 4 \pi \frac{\rho v_t^3}{L_{\rm inj}} \eta_B (t - t_*) \, .
\label{evolution}
\end{equation}


\bigskip

\noindent In the calculations shown in Figure 3b,
we assume a continuous injection of turbulence with $v_t = 500$ km s$^{-1}$ at the injection scale $L_{\rm inj}=1$ Mpc, and $n=\rho/m_p=3\times10^{-3}$ cm$^{-3}$ (with $m_p$ the proton mass). We note that a situation of continuous injection of turbulence with these parameters is appropriate for a dynamically active and massive cluster \cite{dolag+05,donnert+18,hitomicoll18,markevitch+vikhlinin07} and consequently calculations would overestimate the magnetic field strength if extrapolated for a time scale that is much larger than that of a typical merger.

\noindent 
An important parameter in the model is the time when the turbulent dynamo starts, as it fixes the time available for magnetic field amplification up to the epoch of observation. Here, we assume that the turbulent dynamo starts when the clusters have formed a quarter of their mass. This implies that the time available for the amplification is of about 3.7~Gyr \cite{fakhouri+10,giocoli+12}. Assuming a turbulent velocity of $v_t = 500$~km~s$^{-1}$ at an earlier stage in the cluster life-time is optimistic, as the cluster would be too small. Since the turbulent energy flux is proportional to $v_t^3$, the effect of much weaker turbulence at earlier epochs is not relevant for the magnetic field evolution. Conversely, if we assume that the dynamo started later in the cluster life-time, e.g. when it has assembled half of its mass, it results in more stringent Reynolds numbers, since the time available for the  amplification is shorter (about 2.7~Gyr). Lower turbulent velocities and lower number densities would imply higher values of the Reynolds number. On the other hand, higher turbulent velocities would imply a smaller value of Re. However, even if we consider $v_t=800$ km s$^{-1}$ \cite{dolag+05,donnert+18}, i.e. the case of a large turbulent pressure of about 30--40\% the ICM pressure, we obtain Reynolds numbers that are only three times smaller than the case with $v_t = 500$ km s$^{-1}$. We also mention that in real clusters turbulence is induced in a medium that is highly stratified due to gravity \cite{roh+19}. In this case, the transport of magnetised turbulent eddies toward the outskirts could make the turbulent dynamo slightly less efficient than in our model, implying that our constraints on the Reynolds number are conservative.


