Anton Schneider

Chu-Chun Chun Fu

Frédéric Soisson

Cyrille Barreteau

Atomic diffusion in α-iron across the

Keywords:

come

Atomic diffusion plays a central role dictating the kinetics of numerous physical processes in solids, such as surface and interfacial segregation, precipitation and phase transitions. Iron based alloys, being the basis of steels, are certainly one of the most studied systems from both theoretical and experimental points of view. Experimental data for solute and solvent diffusion in iron alloys are usually known at rather high temperatures only (above 750K) [1][2][3][4][5][6][7][8][9][10][START_REF] Speich | The electron Microprobe[END_REF][START_REF] Rothman | [END_REF][13][14], whereas a quantitative modeling of the diffusion coefficients as functions of temperature in these systems is not obvious, since the effects of thermal magnetic excitations and magnetic order-disorder transitions need to be properly described.

Even the simplest case of self-diffusion in bodycentered-cubic (bcc) iron via the vacancy mechanism is still the focus of various recent modelling efforts [15][16][17][18][19]. To the best of our knowledge, there were very few theoretical studies on solute diffusion in iron across the Curie point [15,16,18]. To address these properties at a thermal-vacancy regime, vacancy formation and migration free energies should be determined, as the diffusion of Fe and most substitutional solutes in iron is ruled by first nearest neighbor (1nn) atom-vacancy exchanges.

Density functional theory (DFT) calculations provide an accurate estimation of the vacancy properties in the ground-state ferromagnetic (FM) bcc iron [15-17, 20, 21]. But, theoretical determination of these properties at increased temperatures, when magnetic excitations emerge, becomes non-trivial [22][23][24][25]. So far, many previous studies have estimated the temperature evolution of the diffusion activation energy Q(T ) using the Ruch model [26], with that, the diffusion coefficients

D(T ) = D 0 exp( -Q(T ) k B T
) can be obtained. At a given temperature, Q(T ) is actually the magnetic free energy of activation, which includes the contributions from both the vacancy formation and the various atomvacancy exchange barriers. In the case of self-diffusion, it can be written as a sum of the vacancy formation and the vacancy-atom exchange magnetic free energy (Q(T ) = G f mag + G m mag ). These values are called "magnetic" free energies because only the magnetic entropy is included.

The Ruch model [26] proposes:

Q(T ) = Q P M (1 + αS 2 ) (1)
where Q P M is the activation enthalpy of a perfect paramagnetic (PM) state, where no more magnetic correlation is present. S is the magnetic order parameter (the reduced magnetization) and α is a scalar. Note that to apply the Ruch model, the temperature dependence of S, and either Q P M or the α parameter should be known, in addition to Q F M (the activation enthalpy of the ideal FM state). An example of such application of the Ruch model can be found in Ref. [27].

For a perfect PM state, the vacancy formation and migration enthalpies, and therefore Q P M were often estimated via DFT, by adopting for instance the disordered local moment (DLM) approach within a collinear approximation [17,23], or via an expansion on a set of spin spirals [16,22]. In any case, the magnetic short-range order (MSRO) was generally not considered. Such PM states are therefore expected only at extremely high temperatures. The temperature dependence S for a pure or a dilute system can be easily provided by experiments or simple Ising or Heisenberg models. It is however more difficult to obtain for concentrated alloys with any microstructure.

Concerning intrinsic approximations of the Ruch model [26], it is derived from the Ising model, and due to a mean field approximation, the MSRO effect is absent. Therefore, Q(T ) = Q P M immediately above the Curie point.

Another approach used to determine the temperature evolution of the diffusion properties is the spin-lattice dynamics [18,28], employing empirical potentials and Heisenberg-interaction terms. In particular, a recent study of Wen et al. [18] reported a detailed investigation of self-diffusion in bcc iron. A major advantage of such an approach consists in the natural inclusion of the combined phonon-magnon effects. However, in practice an accurate potential is not obvious to parameterize, especially for magnetic alloys with structural defects. Furthermore, such spin-lattice dynamics simulations can hardly reach very long time scales and large simulation systems (typically a few tens of nanoseconds and 16000 atoms as in Ref. 18).

The present study aims at proposing an efficient and quantitative modelling approach that enables a continuous prediction of diffusion properties versus temperature, including explicitly spin and atomic variables. A DFTbased effective interaction model (EIM) coupled with onlattice Monte Carlo (MC) simulations is adopted. We consider the case of self-and Cu diffusion in bcc iron to illustrate the ability of the methodology to predict diffusion properties, which can be transferred to other magnetic metal alloys.

Some previous studies have already proposed EIMs with both magnetic and chemical variables [22,[29][30][31][32][33][34], but without considering defects and, therefore, not able to study diffusion. On the other hand, more conventional EIMs have been developped to study thermodynamic and kinetic properties without explicitly including spin variables [35][36][37][38]. The present EIM (Eq. 2) consists in a magnetic part with a Landau-Heisenberg form as in Ref. [30,32,33,39]. This allows to account for both longitudinal and transversal excitations of spins. In addition, the pairwise terms (V ij ) capture chemical interactions between atoms.

H = N i (A i M 2 i + B i M 4 i ) + N i P n Zn j J (n) ij M i • M j + N i P n Z j V (n) ij (2)
where Z n is the coordination of the n-th neighbor shell and M i is the magnetic moment of the ith atom. A i and B i are the magnetic on-site coefficients.

J (n) ij and V (n) ij
denote respectively the magnetic exchange-coupling and the chemical-interaction parameters for i and j atoms being n-th nearest neighbors.

First, an EIM for pure iron in a bcc lattice is parameterized on DFT [40][41][42][43][44][45][46][47][48][49][50][51] data. Fitting the magnetic parameters using DFT consists in evaluating the energy difference between systems with similar atomic configurations but distinct magnetic configurations (See [52] for details). We checked that the Curie temperature is correctly reproduced (T C = 1050 K, the experimental value being 1044 K [START_REF] Keffer | Handbuch der Physik[END_REF]). Then, to include the presence of a vacancy (EIM V ), the on-site A and B parameters are modified for atoms located at the first and second nearestneighbor (1nn and 2nn) sites of the vacancy, in order to reproduce the change of their magnetic moment magnitude, while, for simplicity, the J ij remain unchanged. Vacancy formation energies for distinct magnetic spin configurations around the vacancy predicted by DFT are successfully captured by this simple model [52]. In order to simulate the atomic migration, another pure-iron derived EIM SP is also constructed to describe the energetics of an Fe atom at a saddle-point position. In this case, the on-site and J ij parameters of the saddle point atom and their 1nn and 2nn atoms are modified based on DFT data. The atom-vacancy exchange barriers are then determined by the energy calculated using the EIM V and the EIM SP . Note that such way of barrier determination was intensively applied in previous studies using non-magnetic interaction models [35][36][37].

For Cu-diffusion in bcc iron, a Cu atom and a vacancy should be included in the iron system. Similarly, an EIM with all the atoms at lattice positions and another one with an atom (Fe or Cu) at a saddle-point site are parameterized on DFT data on Cu-vacancy binding energies and atom-vacancy exchange barriers with various spin configurations. The numerical parameters of the various EIMs are given in Ref. [52].

The lattice vibrational effects (vibrational entropies and attempt frequencies) are not intrinsically accounted in the present EIM-Monte Carlo set-up but calculated separately by DFT [52]. The magnon-phonon effects [START_REF] Koermann | [END_REF] are therefore not considered.

The tracer diffusion coefficients can be expressed with the Einstein's formula [55][56][57] with < r 2 > and t being the mean square displacement of the tracers and the corresponding physical time:

D * = < r 2 > 6t (3)
For the self-diffusion case, it can also be written in terms of the vacancy concentration and the migration barrier at a given T [START_REF] Adda | La diffusion dans les solides[END_REF][START_REF] Mehrer | Diffusion in Solids[END_REF] as:

D F e * F e = a 2 f 0 C v ν 0 exp( -G m mag k B T ) ( 4
)
where a is the lattice constant, f 0 is the self-diffusion correlation factor (0.727 for a bcc lattice [START_REF] Leclaire | [END_REF]), C v is the equilibrium vacancy concentration, ν 0 is the attempt frequency, G m mag is the magnetic free energy barrier for the vacancy-Fe exchange (vacancy migration), and k B and T are respectively the Boltzmann factor and the absolute temperature. Here

C v = exp( -G f k B T
), with G f being the vacancy formation free energy. Both magnetic and vibrational entropies are considered in this study, and the latter is calculated via DFT at the FM state.

Similarly, the solute (Cu) tracer diffusion coefficient in Fe at the dilute limit can be written [START_REF] Adda | La diffusion dans les solides[END_REF][START_REF] Mehrer | Diffusion in Solids[END_REF] as:

D Cu * F e = a 2 C 1nn v f 2 ν 2 exp( -G m,Cu mag k B T ) (5)
where C 1nn v is the equilibrium vacancy concentration at a 1nn site of the solute, ν 2 is the vacancy-Cu exchange attempt frequency, f 2 is the solute diffusion correlation factor and G m,Cu mag is the magnetic free energy barrier for the vacancy-solute exchange.

We propose a Monte Carlo method which allows to determine the vacancy formation magnetic free energy as a function of temperature. Two separate Fe subsystems are considered with the first one frozen at the FM state, while the magnetic configuration of the second one is allowed to evolve according to temperature. The vacancy is allowed to visit each site of the two systems via the Metropolis algorithm. Then, based on the relative number of visits to the two systems and the vacancy formation energy at the FM state, which is known, the magnetic free energy of vacancy formation versus temperature can be obtained (More details in Ref. 52). Note that, as mentioned in Ref. [18,61,62], a quantum treatment of spins is necessary for a correct prediction of the magnetic entropy, at low temperatures. We therefore adopted the Bose-Einstein statistics in our spin-MC simulations [52] up to the Curie point, following the quasi-harmonic approach of Ref. 61, 63. Then, the Fe and Cu diffusion coefficients are obtained by directly simulating the tracer diffusion experiments with MC simulations [64,65]. We compute the mean square displacement of the tracers (< r 2 >) and the physical time at each temperature (eq. 3).

The physical time t is re-scaled in order to consider the equilibrium vacancy concentration instead of the actual vacancy concentration of the simulation (the diffusion coefficient is multiplied by the factor C v /C M C ). [35,36] During these MC simulations, at each T , we start performing 5 • 10 8 spin Metropolis MC steps to reach the equilibrium magnetic state, then 600 spin steps are performed after each atomic MC step, consisting in a 1nn atom-vacancy exchange based on a time residence algorithm. For simplicity, we assume the typical time spent for one atom-vacancy exchange is sufficiently short, so that all the atomic spins remain frozen while going from the initial to the saddle-point state. However, we have determined that considering another assumption has a negligible effect on the results. Indeed, similar simulations were performed assuming the opposite, being that the spin-variation time is much shorter than the lifetime of both the initial and the saddle-point states, and very close migration barrier were obtained [52]. This test suggests that these properties are not sensitive to the detailed way of implementing the characteristic time of spin variations, and contribute to support the validity of our results.

Via the same MC simulations, we also obtained the magnetic free energy of vacancy migration (G m mag ) in the self-diffusion case, and the magnetic free energy barrier for Cu-vacancy (G m,Cu mag ) and the distinct Fe-vacancy exchanges in the Cu diffusion case. Fig. 1 shows the obtained magnetic free energy for vacancy formation and migration in pure iron, comparing results applying the Bose-Einstein distribution and the Boltzmann statistics for spin-MC at low temperatures. The difference between the two curves is significant, the slope approaching to zero near T = 0 only with the former approach. Therefore the quantum statistics is necessary to obtain a correct low temperature behavior of the vacancy formation and migration entropies (Fig. 1), in agreement with previous spin-lattice dynamics data [18].

The vacancy formation and migration energies at perfect FM and PM states are listed in Table I, together with the resulting activation energies. At the low temperature limit, we reproduce closely the DFT energies (see I). The asymptotic PM energies are obtained at 2000 K. These values for pure Fe are in good agreement with previous DFT and available experiments values [16,52,66] In the FM state, the values for iron and for Cu diffusion are clearly different due to the Cu-vacancy attraction (0.24 eV and 0.17 eV for respectively 1nn and 2nn distances). However, such differences become significantly smaller in the PM state, indicating the dominance of the magnetic disorder over the chemical effect in the very dilute Fe-Cu system. As shown in Fig. 2, the present approach predicts the self-and Cu-tracer diffusion coefficients as functions of temperature, in excellent agreement with experimental studies [1-5, 7, 8, 10-14]. Especially, the sudden deviation from Arrhenius law near the Curie temperature is consistently predicted without any additional assumption. The change of slope (activation energy, Q) between the ferromagnetic and the paramagnetic regimes is also successfully predicted.

FeF M FeP M CuF M CuP M H f 2.
The self-diffusion coefficients obtained with purely classical statistics are also shown for a comparison. It reveals that using the Boltzmann distribution at low T in the spin-MC simulations significantly underestimates the local acceleration of diffusion around the Curie point.

Fig. 2 also shows that the kinetic correlation factor f 2 for Cu diffusion increases with temperature up to an asymptotic limit of 0.73, which is the f 0 value in pure bcc iron. To clarify the role of magnetic disorder on the f 2 , we have performed similar MC simulations for Cu diffusion, but imposing a perfect FM order for all atomic-MC temperatures. The results show that the kinetic correlation factor of Cu diffusion increases more slowly when magnetic disorder is absent, This comparison together with the smaller difference between different barriers at the PM than at the FM state (Table I and Fig. 2) suggest the dominance of the magnetic disorder over the chemical interaction effect in this very dilute Fe-Cu system. Indeed, the magnetic free energy of binding between a vacancy and a Cu atom at 1nn decays from 0.24 eV in the FM state (consistently with our DFT results) to 0.09 eV in the PM state. It is worth noting that if applying the Ruch model using the currently obtained Q P M and Q F M and S(T ), we obtain very similar diffusion coefficients. For a closer comparison between results from our method and by using the Ruch model for an interpolation, Fig. 3 shows the respective data for the activation free energy, and the magnetic free energy of vacancy formation and migration. As can be seen, both methods give very close values, especially above the Curie temperature. The largest differences occur between the formation free energies at T < T C . In any case, the discrepancies are smaller than 10 %. This good agreement between our approach and the Ruch model shows that the effects of magnetic short range order on vacancy properties are rather limited. In the case of Cu solute diffusion, the Ruch model also provides very close results. [52] These comparisons suggest that the simple Ruch model allows a very good description of the temperature evolution of activation, formation and migration free energies in a ferromagnetic system. However, it should be noted that the Ruch interpolation requires an accurate knowledge of the asymptotic (FM and PM) energetic values, and the temperature evolution of the magnetization. They are generally not obvious to obtain both experimentally or from ab initio calculations, especially for alloys beyond the dilute limit. Also, if estimating the PM energies for highly itinerant magnetic systems, such as fcc Ni or bcc Cr via DFT, local-magnetism constraints should be applied, which can be extremely time consuming [23,68].

In summary, we propose an approach to efficiently predict atomic diffusion properties in iron, by performing on-lattice Monte Carlo simulations using effective interaction models parameterized on DFT data. These EIMs contain explicitly both chemical and magnetic variables, and the presence of a vacancy.

This approach naturally accounts for the interplay between magnetic and chemical degrees of freedom.

It is shown to successfully predict the temperature evolution of vacancy formation and migration magnetic free energies, and the tracer diffusion coefficients of Fe and Cu in bcc iron, across the Curie temperature. This approach is also ready to address the diffusion as a function of solute concentrations, and in the presence of nonequilibrium vacancies. The same approach is fully transferable to other magnetic metal systems.

At variance with the DFT-Ruch method, the current approach predicts properties for all temperatures regardless of the magnetic state. The crucial issue is to accurately parameterize the EIMs. On the other hand, it allows to reach to a calculation time of several orders of magnitudes longer than the spin-lattice dynamics. Therefore, it is also promising to address more complex kinetic processes than the atomic diffusion, such as the ordering, precipitation or segregation.

This work was partly supported by the French-German ANR-DFG MAGIKID project. Ab-initio calculations were performed using DARI-GENCI resources under the A0070906020 project and the CINECA-MARCONI supercomputer within the SISTEEL project.

FIG. 1 :

 1 FIG. 1: Left panel: Magnetic free energy, enthalpy and entropy of formation. Right panel: Magnetic free energy, enthalpy and entropy of migration.

  FIG. 2: Upper-left panel: Self-diffusion coefficients of Fe versus T, with and without quantum effects and comparison with experimental data obtained from [1-5, 7, 8, 10]. Bottom left panel: Diffusion coefficients of Fe and Cu, compared with experimental results, obtained from [1-5, 7, 8, 10] for Fe and [4, 14, 67] for Cu. Upper-right panel: Kinetic correlation factors of Fe self-diffusion (f0) and Cu diffusion (f2). The f2 kinetic correlation with the frozen FM state (see text) is also displayed. Bottom right panel: Migration magnetic free energies of various jumps.

FIG. 3 :

 3 FIG. 3: Activation, formation and migration free energies compared with the Ruch model.


S A Torchinsky

J.-Ch Hamilton
email: hamilton@apc.in2p3.fr

M Piat
email: piat@apc.univ-paris7.fr

E S Battistelli

C Chapron

G D'alessandro

P De Bernardis

M De Petris

M M Gamboa Lerena

M González

L Grandsire

S Marnieros

S Masi

A Mennella

L Mousset

J D Murphy

C O'sullivan

D Prêle

G Stankowiak

A Tartari

J.-P Thermeau

F Voisin

M Zannoni

P Ade

J G Alberro

A Almela

G Amico

L H Arnaldi

D Auguste

J Aumont

S Azzoni

S Banfi

B Bélier

A Baù

D Bennett

L Bergé

J.-Ph Bernard

M Bersanelli

M.-A Bigot-Sazy

J Bonaparte

J Bonis

E Bunn

D Burke

D Buzi

F Cavaliere

P Chanial

R Charlassier

A C Cobos Cerutti

F Columbro

A Coppolecchia

G De Gasperis

M De Leo

S Dheilly

C Duca

L Dumoulin

A Etchegoyen

A Fasciszewski

L P Ferreyro

D Fracchia

C Franceschet

K M Ganga

B García

M E García Redondo

M Gaspard

D Gayer

M Gervasi

M Giard

V Gilles

Y Giraud-Heraud

M Gómez Berisso

M Gradziel

M R Hampel

D Harari

S Henrot-Versillé

F Incardona

E Jules

J Kaplan

C Kristukat

L Lamagna

S Loucatos

T Louis

B Maffei

W Marty

A Mattei

A May

M Mcculloch

L Mele

D Melo

L Montier

L M Mundo

J A Murphy

F Nati

E Olivieri

C Oriol

A Paiella

F Pajot

A Passerini

H Pastoriza

A Pelosi

C Perbost

M Perciballi

F Pezzotta

F Piacentini

L Piccirillo

G Pisano

M Platino

G Polenta

R Puddu

D Rambaud

E Rasztocky

P Ringegni

G E Romero

J M Salum

C G Scóccola

A Schillaci

S Scully

S Spinelli

M Stolpovskiy

A D Supanitsky

P Timbie

M Tomasi

G Tucker

C Tucker

D Viganò

N Vittorio

F Wicek

M Wright

A Zullo

QUBIC III: Laboratory Characterization

A prototype version of the Q & U Bolometric Interferometer for Cosmology (QUBIC) underwent a campaign of testing in the laboratory at Astroparticle Physics and Cosmology in Paris. We report the results of this Technological Demonstrator which successfully shows the feasibility of the principle of Bolometric Interferometry. Characterization of QUBIC includes the measurement of the synthesized beam, the measurement of interference fringes, and the measurement of polarization performance. A modulated and frequency tunable millimetre-wave source in the telescope far-field is used to simulate a point source. The QUBIC pointing is scanned across the point source to produce beam maps. Polarization modulation is measured using a rotating Half Wave Plate. The measured beam matches well to the theoretical simulations and gives QUBIC the ability to do spectro imaging. The polarization performance is excellent with less than 0.5% cross-polarization rejection. QUBIC is ready for deployment on the high altitude site at Alto Chorillo, Argentina to begin scientific operations.

Introduction

The detection of B-mode polarization in the Cosmic Microwave Background (CMB) is the subject of a worldwide effort due to its importance as a confirmation of the Inflationary model of Cosmology. A clear detection of polarization B-modes in the CMB is evidence of primordial gravitational waves expected during the inflationary phase in the earliest moments of the Universe. For an overview see [START_REF] Hu | A CMB polarization primer[END_REF]. Polarization in the CMB is expected to be a small fraction of the CMB signal itself while the B-mode component of polarization will be several orders of magnitude smaller again. As a result, the measurement of B-mode polarization is a difficult exercise of extracting a signal buried deep within other signals (the CMB E-mode but also foreground contamination and lensing) and noise. This requires very careful instrument design and implementation, with sophisticated data analysis to separate unwanted contributions from the B-mode signal.

The Q and U Bolometric Interferometer (QUBIC) is designed with particular attention to the limitation and control of systematics [START_REF] Aumont | QUBIC Technical Design Report[END_REF]. See also in this series of papers [START_REF] O'sullivan | QUBIC VIII: Optical design and performance[END_REF] for the optics design, [START_REF] Masi | QUBIC V: Cryogenic system design and performance[END_REF] for the cryogenics design, and [START_REF] Piat | QUBIC -IV: Performance of TES Bolometers and Readout Electronics[END_REF] for the readout electronics. QUBIC uses the technique of interferometry which leads to the possibility of doing "self-calibration" in order to have exquisite control of instrument systematics [START_REF] Bigot-Sazy | Self-calibration: an efficient method to control systematic effects in bolometric interferometry[END_REF]. This technique also provides a method for doing spectral imaging. Bolometric Interferometry is the marriage of techniques bringing together the great sensitivity and large bandwith of bolometers and the instrumental control and high fidelity imaging of aperture synthesis. Using this innovative approach, any residual systematics in the data will be largely independent from those in other experiments, thus providing a uniquely powerful dataset in the context of the worldwide experimental effort.

An imaging interferometer measures "visibilities" which are the complex (amplitude and phase) correlations between each antenna pair (baseline). In radio astronomy, the visibilities are recorded directly. A "correlator" digitizes the signals and multiplies pairs of signals to produce a stream of complex numbers, each of which corresponds to the cross correlation product of an antenna-pair. Channelization of the bandpass permits signal processing of individual, very narrow bands, and for each channel the signal is nearly monochromatic. In radio astronomy, large bandwidths are achieved by adding more digital electronics.

A bolometric interferometer takes advantage of the high sensitivity and large bandwidth of bolometers while also benefitting from the calibration technique possible with an imaging interferometer. The spatial sampling of the sky is generated by placing a cluster of back-toback horns that behave effectively as electromagnetic nozzles. This horn cluster creates the u -v sampling of the aperture plane equivalent to what is done by a distribution of antennas in a radio array. For the bolometric interferometer, instead of sampling the signals and computing the cross correlations between antenna pairs, the interference pattern is imaged.

A single image of the interference pattern has all the information convolved together resulting in observing the sky through a synthesized beam. The shape of this synthesized beam is given by the combination of all individual baselines (all pairs of horns). The bolometric interferometer end up being a synthesized imager observing the sky through its synthesized beam just the same way as a classical imager observes the sky through the beam formed by the telescope. For calibration and instrumental systematics studies it is however crucial to extract the individual visibilities. By blocking all horns except two, we measure the interference pattern of that baseline [see 6]. For example, using a 20 × 20 cluster of horns as in the QUBIC Final Instrument (QUBIC-FI), there are 400 horns making n(n -1)/2 = 79800 baselines which are needed to be observed individually for self-calibration. The QUBIC Technological Demonstrator (QUBIC-TD) under test in the lab at APC has a smaller horn array with only 64 horns in a square array giving 2016 baselines.

Performance of the bolometric interferometer improves as the number of baselines increases. A larger cluster of horns provides more baselines, but this in turn must be sampled by a larger array of detectors in the focal plane. As a result, the overall sensitivity of the bolometric interferometer is indirectly a function of the number of detectors in so much as the focal plane must be large enough to match the horn cluster.

An additional feature of bolometric interferometry is the ability to do spectroimaging. The synthesized beam varies with frequency [for details see 3]. In particular, the beam secondary lobes move inward as frequency increases while the central lobe remains in the same place. As a result, the Time Ordered Data effectively samples different frequencies as the beam passes over the same point in the sky. This frequency selectivity can be deconvolved in the data post processing. The spectral resolution improves with the number of baselines as the synthesized beam has finer secondary lobes with a larger cluster of horns. Spectroimaging is an innovative feature of bolometric interferometry which gives QUBIC an important advantage over other CMB imagers [see 7, 8, for details] and the possibility to do this with QUBIC has been confirmed (see Section 7).

This paper is organized as follows. The laboratory setup is described in Section 2 with particular attention given to the placement and alignment of the calibration source. Section 3 presents the measured bandpass of QUBIC. Section 4 shows the measured response to modulated polarization using the rotating Half Wave Plate. The QUBIC design provides extremely pure polarization performance with 99.9% cross polarization rejection. Section 5 shows the measurement of fringes on the focal plane using the mechanical switches in the horn array to select baselines. Section 6 shows an early application of self testing with the independent determination of the physical orientation of the horn array using the fringe measurements. In Section 7 we show measurements of the QUBIC synthesized beam, and in Section 8 we demonstrate the results of the QUBIC mapmaking pipeline by generating a "sky" map of the calibration source using the calibration information determined from previous measurements. Finally, some concluding remarks are given in Section 9.

Calibration Source Setup

Characterization of the QUBIC instrument is done primarily using a frequency tunable monochromatic point source in the far-field. This permits measurement of the bandpass (Section 3), the polarization performance (Section 4), the measurement of interference fringes (Section 5 and Section 6), and the beam Point Spread Function (PSF -Section 7 and Section 8). This section describes the optical setup used to characterize QUBIC.

The setup is shown in the sketch of Figure 1. The calibration source points at a flat mirror which redirects the beam into the QUBIC cryostat window. There is an 11 m optical path, putting the calibration source comfortably in the far field. The flat mirror is an aluminium sheet mounted on a scaffolding at a height of 3.5 m (photo in Figure 2). The tilt angle of the flat mirror can be adjusted by a long screw permitting very precise selection for the correct tilt angle. Alignment of the system was accomplished using a laser temporarily mounted at the window of the QUBIC cryostat, pointing normal to the window (photo Figure 3). The laser light is reflected from the flat mirror to the calibration source where a small flat mirror was fitted to the front of the calibration source feedhorn. The laser reflects from the mirror at the feedhorn mouth and returns to the large flat mirror and finally to the QUBIC cryostat window. In the photo (Figure 3), one can clearly see the spot of the laser on the corner of the laser mounting structure on the window. The alignment is therefore precise to within a fraction of a degree which is well within the tolerance necessary to have the calibration source visible to the QUBIC-TD.

The calibration source system is composed of a number of electronics aside from the source itself. The calibration source was purchased from VDI electronics and is composed of a Gunn oscillator running at frequencies around 10 GHz. Two multipliers in a chain follow the synthesiser multiply the base frequency by 12 resulting in frequencies around 150 GHz. The range of the system is between 130 GHz and 170 GHz. The source can be commanded via USB connection to select the frequency. This is done by a nearby Raspberry Pi mini computer (see photo Figure 4).

A signal generator provides a square wave at around 1 Hz which is used to modulate the calibration source. The Raspberry Pi configures the signal generator, selecting the amplitude, offset, frequency, shape, and duty cycle.

The power output from the calibration source is monitored at a port near the calibration source feedhorn which provides a voltage which is a function of the output power. This voltage is sent to an Analog-to-Digital Converter (ADC) board integrated in the Raspberry Pi mini computer which stores the data along with a timestamp for each sample. There are approximately 300 samples per second.

A Command Line Interface written in python is used to configure the calibration source setup. This system accepts commands via socket and can be easily interfaced by the Graphical User Interface called QubicStudio.

For the result shown in Figure 5, the source was set at 150 GHz and was modulated with a square wave with 3 second period (0.333 Hz) and a 33% duty cycle. Figure 5 shows the very clear signal in many of the TES detectors in the array. Figure 6 shows the signal seen on one TES together with the modulation signal measured by the power monitor. There is a very clear correlation between the two. Photo of the electronics associated with calibration source setup. This includes a signal generator for modulating the calibration source with a square wave. The output power of the calibration source is sent through an amplifier and then digitized by an ADC board integrated in the Raspberry Pi (seen suspended on the cables below the calibration source).

Spectral Response

The spectral response was measured using the calibration source (see Section 2). The calibration source was setup with a sine wave modulation at 1 Hz and stepped through frequencies across the band from 110 GHz to 180 GHz. After synchronized demodulation at the bolometer, this provides an extremely high signal-to-noise measurement. Figure 7 shows the measured profile which corresponds to the expected profile of the filter provided by Cardiff. The plot shows the average of all operational TES [181 detectors, see 5], as well as the median, and the points filtered for minimum 3σ detection. All curves agree within the bandpass.

Half Wave Plate Polarization Rotation Test

A functionality test was carried out of the Half Wave Plate rotator mechanism, and at the same time, the calibration source was operating. Measurements were taken at each of the 7 evenly spaced positions of the Half Wave Plate from 0 • to 90 • (spacing of 15 • ).

Figure 8 shows the signal, after Fourier filtering, for TES #95 at the different HWP positions. The peak-to-peak amplitude clearly varies with HWP position. Each position is shown in a different colour. Figure 9 shows a zoom of the calibration source signal measured by TES#95 with the HWP in position #1. The amplitude at each position is plotted and fitted to a sine curve, as shown in Figure 10.

The cross-polarization contamination at 150 GHz is compatible with zero to within 0.18%. The signal in each case was measured by the RMS of the TES data in a given position while the source was modulated. The measurement was the source modulation amplitude together with the RMS of the noise, quadratically added. As a result the minimum value is dominated by noise. The maximum signal occurs when the HWP is between position 6 and 7.

As polarization is selected before the horn array, the cross-polarization leakage is independent of the location in the focal plane. Therefore, the cross-polarization leakage is expected to be the same for all TES. The example shown in Figures 10 and11 are for TES#95 which has the best signal-to-noise for this measurement due to the pointing.

The measured cross-polarization contamination was measured with the calibration source tuned to 150 GHz. This is therefore a measurement at a single, essentially monochromatic frequency. Future measurements will measure the cross-polarization contamination at frequencies across the band in order to determine the integrated response of the system for the full bandpass.

Self-Calibration

Self Calibration is a technique developed for aperture synthesis in radio interferometry. This technique evolved from the original idea in the 1970's of "phase-closure" [START_REF] Högbom | Aperture Synthesis with a Non-Regular Distribution of Interferometer Baselines[END_REF][START_REF] Fort | A Method of Obtaining Brightness Distributions from Long Baseline Interferometry[END_REF][START_REF] Readhead | The mapping of compact radio sources from VLBI data[END_REF] to become in the early 1980's "self-calibration" [START_REF] Cornwell | A new method for making maps with unstable radio interferometers[END_REF][START_REF] Pearson | Image Formation by Self-Calibration in Radio Astronomy[END_REF]. See [START_REF] Cornwell | Self-Calibration[END_REF] for a detailed overview. Precise knowledge of the calibration source is not required, as long as it is a strong and stable point source. The large number of baseline visibilities allows us to solve for many unknowns, including the gain and phase corrections required, without having knowledge of the source amplitude.

In order to advance towards the full analysis of self-calibration, a key component is the ability to measure fringes with QUBIC. If fringes can be measured with a horn pair, then the full analysis can be done once the fringe measurement is done for all horn pairs. By measuring the fringe pattern of a single pair of horns we demonstrate the feasibility of doing self-calibration with the bolometric interferometer.

The QUBIC Technological Demonstrator successfully measured fringes between a pair of horns (Figure 13 left). This is the derived image after analysing measurements of images with all horns open, with two horns closed, with one horn closed, and with the other horn closed. This is the equivalent of having all horns closed except the two, as shown in [START_REF] Bigot-Sazy | Self-calibration: an efficient method to control systematic effects in bolometric interferometry[END_REF]. The fringes are expected to be fainter in proportion to the distance to the centre of the focal plane (Figure 13 right). This is not the case here because of saturation of the TES detectors. As a result, the fringe amplitude appears to be relatively constant, or near zero where saturated detectors were subtracted from one another. The problem with saturation was due to the cryogenic system which did not reach the optimum temperature for the TES. Future tests will work below the saturation point of the TES. An initial analysis using the fringe measurement is carried out in Section 6. The measurement is used to verify the orientation of the horn array.

Verification of the Horn Cluster Orientation with Fringe Measurements

The orientation and spacing of the fringe pattern seen in the focal plane is a function of the selected baseline which is determined by which horns are open or closed. In this section, we describe the use of the fringe measurement at different baselines to confirm the mechanical orientation of the horn cluster.

The fringe pattern described in the previous section was measured for the baseline formed by horns 25-57 (see Figure 13 above). The position of the horn cluster at that time was verified using pictures taken when it was mounted (see Figure 14). Figure 15 shows the reference frames used on the instrument.

Simulations of interferometry patterns obtained on the focal plane for different horn configurations are done using geometric optics and calculating the interference at infinity. This is computed as if the horn plane and the focal plane were parallel. Simulations are also performed in Maynooth with GRASP (a commercial optical modelling software) and MODAL (an 'in-house' software developed and maintained by Marcin Gradziel at Maynooth) [see 3, for details about optics modelling]. While both software packages model systems using the technique of Physical Optics, there are differences in how systems are implemented which make them complementary. MODAL can calculate the coupling of beams from free-space to feedhorns, while GRASP considers feedhorns to be point sources. However, GRASP can describe surfaces using a grid of points which permits the modeling of real measured mirror surfaces, and then using the Physical Theory of Diffraction, GRASP can model mirror edge effects. The two packages can exchange data through common file formats which permits the analysis of a complete optical system, using the best package for each subsystem and then combining the results. They can also be used independently and compared as a verification of modelling results.

A model was created in GRASP to match the MODAL model. This model plots the interference pattern of any horn baseline. The comparison of the output of the models for a given baseline is shown in Figure 16 as plotted by our Python-based QUBIC simulation software using the output files from the MODAL simulation model. The result shows that MODAL and the GRASP simulations are in agreement. Additional baselines were tested and also had matching patterns. Having these two independent methods in agreement is a good check that the horn setup, instrument data, software, and simulations are working as expected. Finally, the optics simulations can be used to confirm the physical horn orientation in the optical path. The orientation between all the simulations and the measurement agree. A π/2 rotation difference between measured results and simulations would be easily detected. However, it is not straightforward to detect a π rotation because the differences would be very small (see Figure 17). Such a rotation might be detectable by fitting the measurement with the simulations from Maynooth, and is a work in progress.

The main concern which motivated this study was to confirm that the horn switch electronics [START_REF] Cavaliere | QUBIC VII: The feedhorn-switches system of the technological demonstrator[END_REF] was operating the switches as expected. That is, that the horn commanded to be open or closed was in the correct physical location in the horn cluster. This was confirmed.

Synthetic full Beam reconstruction

The QUBIC-TD synthetic beam maps were measured at five frequencies in the range 130 GHz to 170 GHz by tuning the VDI calibration source (see Section 2) in steps of 10 GHz from 130 GHz. For each frequency measurement the calibration source was modulated at a period of 1 second with a sinusoidal profile, and the TES signal was demodulated in post processing. QUBIC-TD was configured to point at an elevation angle of 35 • and to scan across azimuth at a constant rate from -25 • to +25 • . The elevation angle was then increased by 0.2 • and a new azimuth scan at constant rate was done in the reverse direction from +25 • to -25 • . Each azimuth scan takes nearly 8 minutes, and an entire beam map at a given frequency is completed in 22 hours and 30 minutes. This matches well with the cryogenic hold time of the cryostat. The entire beam map measurement was preprogrammed in a script executed by the QUBIC instrument control software QubicStudio.

The result of the beam mapping measurement campaign is a series of maps for each TES pixel in the detector array. This is 244 maps for each frequency for a total of 1220 maps. Figure 18 shows example maps for TES#93 (ASIC 1) at three frequencies. The main lobe and secondary lobes are clearly visible and match well with theory. The secondary lobe locations depend on frequency while the main lobe is always at the same place. Figure 19 shows a cut across the beam maps of Figure 18. The secondary lobes are closer to centre as frequency increases, as expected. This is a key feature of the bolometric interferometer which makes spectroscopic imaging possible.

Map Making with Measured Synthesized Beams

Using the relative location and amplitude of all the peaks in the synthesized beam for each of the TES, we can now project the data onto the sky using optimal map-making to deconvolve from the effect of the multiple peaks. When observing the sky, this will result in an unbiased CMB map as has been shown using simulations [see 8, in this series of papers]. We have performed this with the calibration data in order to obtain an image of the point-like calibration source we have been using. The resulting image is shown in Figure 20 exhibiting the expected point-source shape with a FWHM of 0.68 degree in excellent agreement with expectations. The successful mapmaking with the measured synthesized beam is effectively an endto-end checkout of the entire QUBIC system. In order for this exercise to be possible, all subsystems, interfaces between subsystems, and all associated software must be functioning correctly. This includes scientific and housekeeping data acquisition, telescope pointing control, and control and synchronisation of all subsystems. The software needed to make this measurement includes the system control software, the data acquisition software, data archiving and reading, and finally data analysis software together with comparison to system simulation software. The resulting map of Figure 20 shows that all subsystems are functioning correctly, and all subsystems are correctly managed and synchronized together into the overall system.

Conclusion

QUBIC underwent a long campaign of testing in the laboratory at APC in Paris. These measurements confirm that Bolometric Interferometry is a viable method for the measure of CMB B-mode polarization. In particular, using an artificial millimetre-wave source in the telescope far-field, a number of fundamental requirements were successfully demonstrated. The measured synthesized beam matches well with theory, and demonstrates that spectral imaging is possible. The spectral response matches the expected bandpass given the optical components in the optical chain. The polarization performance has less than 0.2% cross polarization contamination at 150 GHz which will make QUBIC the best instrument for polarization purity among the currently running CMB B-mode experiments.

The emphasis on QUBIC design has been on the suppression and control of systematic effects. QUBIC employs bolometric interferometry and a conservative polarization optical design which together has the advantage of exquisite control of systematics while still giving the high sensitivity inherent in the use of wide band bolometers.

Deployment of QUBIC on the scientific site of Alto Chorillo at 5000 m altitude in Argentina is expected during 2020. QUBIC will provide an extremely clean polarization map of the sky and together with spectro imaging, will have excellent separation of foreground sources.

Figure 1 .

 1 Figure 1. Layout for the calibration source relative to the QUBIC instrument.

Figure 2 .

 2 Figure 2. left: Photo of QUBIC looking along the line-of-sight from the calibration source. The reflection of the window is clearly visible in the flat mirror. right: The flat mirror is mounted on a scaffold at a height of 3.5 m and has a finely adjustable tilt angle using a system composed of a long screw.

Figure 3 .

 3 Figure 3. Photo of the alignment procedure. The laser is mounted orthogonal to the cryostat window and shines towards the flat mirror, sending the light to the calibration source feedhorn at the other side of the room. A small mirror fitted to the mouth of the feedhorn sends the laser light back on the same path to the cryostat window. The laser spot is clearly visible on the laser mount structure.

Figure 4 .

 4 Figure 4.

Figure 5 .

 5 Figure 5. The first measurement of the calibration source by QUBIC. The signal is clearly seen in most pixels, corresponding to the good pixels in the array [see 5, for a discussion on the bolometer array performance]. The black "pixels" in the top-right are empty positions. The QUBIC Final Instrument will have four arrays equivalent to this one in order to make a roughly circular focal plane for each frequency channel.

Figure 6 .

 6 Figure 6. Overlay of the signal detected by TES#54 of the QUBIC array (blue curve) together with the modulated calibration source signal as measured by the calibration source power monitor (red curve).

Figure 7 .

 7 Figure 7. Bandpass of the QUBIC Technological Demonstrator. This was measured by stepping through the frequencies of the calibration source and measuring the relative power on the TES at each frequency.

Figure 8 .

 8 Figure 8. Peak-to-peak amplitude of the detected calibration source for different Half Wave Plate positions as measured with TES#95

Figure 9 .

 9 Figure 9. Peak-to-peak amplitude of the detected calibration source for the Half Wave Plate in position #1 as measured with TES#95. The blue curve is the response measured by the bolometer and the orange curve is the power monitor of the calibration source (see Section 2).

Figure 10 .

 10 Figure 10. Amplitude of the detected calibration source at 150 GHz as measured with TES#95 for the Half Wave Plate in the different positions and fitted to a sine curve.

Figure 11 .

 11 Figure 11. Same as Figure10but with the calibration source set at various power output levels. The effect of saturation is evident for the source at high power. The Half Wave Plate angle of maximum signal damping is consistent between all the measurements.

Figure 12 .

 12 Figure 12. Distribution of the trough in the amplitude for various source output power (see Figure 11).

Figure 13 .

 13 Figure 13. left: Fringes measured on the QUBIC TES array. Fringes are clearly visible as bright, diagonal lines across the detector array. right: Simulation of the fringe pattern expected to be measured with a baseline between horns 25 and 57. The fringe lines are in the same place as the measured image (left). The amplitude of the fringe lines are fainter with distance from the centre of the focal plane. This was not measured because of saturation of the detectors.

Figure 14 .

 14 Figure 14. Pictures showing the the horn cluster orientation. left: The On Axis Focal Plane (ONAFP) reference frame projected onto the face of the horn cluster. right: The Global Reference Frame (GRF) is co-planar to the face of the horn cluster.

Figure 15 .

 15 Figure 15. QUBIC reference frames. left: The face of the horn cluster is overlaid on the optics schematic. right: The fringe measurement in the negative-x and negative-y quadrant (quadrant 3).

Figure 16 .

 16 Figure 16. Fringes obtained on the focal plane for given baselines using our Python-based QUBIC simulation software (optical geometry) using results from the MODAL simulatin. The corresponding GRASP simulation is shown on the right.

Figure 17 .

 17 Figure 17. Fringes obtained on the focal plane for 2 baselines, symmetric under a π rotation.

Figure 18 .

 18 Figure 18. Maps of the synthetic beam for TES#93 (ASIC 1). The multilobe synthetic beam is shown here for 130 GHz, 150 GHz, and 170 GHz. An animated version of this plot with comparison to simulation can be viewed online at https://box.in2p3.fr/index.php/s/bzPYfmtjQW4wCGj

Figure 19 .

 19 Figure 19. A cut across the synthetic beam maps shown in Figure 18. The plot on the left shows the coincident central lobe for each frequency and the secondary lobes which are closer to centre as frequency increases. The image on the right indicates the orientation of the slice.

Figure 20 .

 20 Figure20. Calibration data with the source at 150 GHz projected on the sky using our map-making software to deconvolve from the multiple peaked synthesized beam. Small residuals of the deconvolved multiple peaks can be seen and are due to unmodeled non-linearity of the detectors due to saturation with the strong calibration source.

Acknowledgments QUBIC is funded by the following agencies. France: ANR (Agence Nationale de la Recherche) 2012 and 2014, DIM-ACAV (Domaine d'Intérêt Majeur-Astronomie et Conditions d'Apparition de la Vie), CNRS/IN2P3 (Centre national de la recherche scientifique/Institut national de physique nucléaire et de physique des particules), CNRS/INSU (Centre national de la recherche scientifique/Institut national et al de sciences de l'univers). Italy: CNR/PNRA (Consiglio Nazionale delle Ricerche/Programma Nazionale Ricerche in Antartide) until 2016, INFN (Istituto Nazionale di Fisica Nucleare) since 2017. Argentina: MINCyT (Ministerio de Ciencia, Tecnología e Innovación), CNEA (Comisión Nacional de Energía Atómica), CON-ICET (Consejo Nacional de Investigaciones Científicas y Técnicas). D. Burke and J.D. Murphy acknowledge funding from the Irish Research Council under the Government of Ireland Postgraduate Scholarship Scheme. D. Gayer and S. Scully acknowledge funding from the National University of Ireland, Maynooth. D. Bennett acknowledges funding from Science Foundation Ireland.


Sanne Govaert
email: sanne.govaert@ugent.be

Camille Meeussen

Thomas Vanneste

Kurt Bollmann

Jörg Brunet

Sara A O Cousins

Martin Diekmann

Bente J Graae

Per- Ola Hedwall

Thilo Heinken

Giovanni Iacopetti

Jonathan Lenoir

Sigrid Lindmo

Anna Orczewska

Michael P Perring

Quentin Ponette

Jan Plue

Federico Selvi

Fabien Spicher

Matteo Tolosano

Pieter Vermeir

Florian Zellweger

Kris Verheyen

Pieter Vangansbeke

Pieter De Frenne

Edge influence on understorey plant communities depends on forest management Edge influence depends on management

Keywords: Edge effects, edge influence, forest specialists, generalists, herbaceous layer, patch contrast, plant biodiversity, species richness, thinning, understorey

Does the influence of forest edges on plant species richness and composition depend on

Introduction

Forests harbour more than two-thirds of terrestrial biodiversity (Millennium Ecosystem Assessment, 2005). In Europe, dramatic deforestation in the past has led to habitat loss and fragmentation, both being major drivers of biodiversity loss (Secretariat of the Convention on Biological Diversity (CBD), 2010). Habitat fragmentation creates new forest edges and subsequently leads to a higher edge to interior ratio [START_REF] Saunders | Biological consequences of ecosystem fragmentationa review[END_REF]. At least 20 % of the world's forest area lies within 100 m of a forest edge [START_REF] Haddad | Habitat fragmentation and its lasting impact on Earth's ecosystems[END_REF], indicating the need to further our understanding of how edges affect forest biodiversity and functioning. Compared to forest interiors, forest edges experience more solar radiation, higher wind speeds and faster air mixing, resulting in higher light availability, decreased soil moisture, higher maximum and lower minimum temperatures, and increased diurnal and seasonal variability in temperatures [START_REF] Chen | Microclimate in forest ecosystem and landscape ecology -Variations in local climate can be used to monitor and compare the effects of different management regimes[END_REF][START_REF] Gehlhausen | Vegetation and microclimatic edge effects in two mixed-mesophytic forest fragments[END_REF][START_REF] Matlack | Microenvironment variation within and among forest edges in the eastern United States[END_REF][START_REF] Schmidt | Microclimate and matter dynamics in transition zones of forest to arable land[END_REF][START_REF] Tuff | A framework for integrating thermal biology into fragmentation research[END_REF]. In addition to altered microclimatic conditions, forest edges, as opposed to forest interiors, are also characterized by a higher seed influx of non-forest species [START_REF] Devlaeminck | Inflow of seeds through the forest edge: evidence from seed bank and vegetation patterns[END_REF], differences in disturbance regimes, and higher nitrogen (N) and carbon (C) stocks [START_REF] Remy | Strong gradients in nitrogen and carbon stocks at temperate forest edges[END_REF]. These biotic and abiotic factors all influence understorey plant communities, which contain more than 80 % of total plant species richness in temperate forests and are essential for several ecosystem functions such as nutrient cycling, carbon dynamics and tree regeneration [START_REF] Gilliam | The ecological significance of the herbaceous layer in temperate forest ecosystems[END_REF]Landuyt et al., in press;[START_REF] Whigham | Ecology of woodland herbs in temperate deciduous forests[END_REF].

In general, higher understorey species richness occurs at forest edges because those environmental conditions are intermediate between the forest interior and the matrix. But not all species respond similarly, and species can be grouped according to their ecological preferences [START_REF] Gehlhausen | Vegetation and microclimatic edge effects in two mixed-mesophytic forest fragments[END_REF][START_REF] Guirado | Understorey plant species richness and composition in metropolitan forest archipelagos: effects of forest size, adjacent land use and distance to the edge[END_REF][START_REF] Murcia | Edge effects in fragmented forests: implications for conservation[END_REF]. For example, a study in temperate forests in northern France found that edge species were more often thermophilous (warmth-loving), basophilous (alkaline-loving), nitrogen-and light-demanding, and species reproducing by seed only (and not clonally) [START_REF] Pellissier | Understorey plant species show long-range spatial patterns in forest patches according to distanceto-edge[END_REF]. In contrast, interior species were more often slow-colonizing species, acidophilous and species reproducing both by seeds and vegetatively [START_REF] Pellissier | Understorey plant species show long-range spatial patterns in forest patches according to distanceto-edge[END_REF]. The effect of distance to edge, hereafter also referred to as 'edge influence', on species richness has been related to mechanisms such as past land-use change [START_REF] Berges | Past landscape explains forest periphery-to-core gradient of understorey plant communities in a reforestation context[END_REF], changes in light regime (Honnay, [START_REF] Hijmans | raster: Geographic Data Analysis and Modeling. R package version 2.8-19[END_REF] and altered soil moisture content [START_REF] Gehlhausen | Vegetation and microclimatic edge effects in two mixed-mesophytic forest fragments[END_REF].

While edge influence has been well studied in dense forests in single regions or landscapes, less is known on how edge influence interacts with forest management across forest types in Europe. Similar to habitat fragmentation, forest management also affects forest microclimate, and thus light availability, temperature, soil moisture and nutrient availability [START_REF] Ash | Changes and variability in the field layer of a coppiced woodland in Norfolk, England[END_REF][START_REF] Grayson | Understory light regimes following silvicultural treatments in central hardwood forests in Kentucky, USA[END_REF]. Understorey plant communities respond to changes in light and soil disturbances which can be induced by forest management practices [START_REF] Ash | Changes and variability in the field layer of a coppiced woodland in Norfolk, England[END_REF][START_REF] Aude | Vegetation in Danish beech forests: the importance of soil, microclimate and management evaluated by variation partitioning[END_REF][START_REF] Decocq | Plant diversity in a managed temperate deciduous forest: understorey response to two silvicultural systems[END_REF][START_REF] Scolastri | Old coppice versus high forest: the impact of beech forest management on plant species diversity in central Apennines (Italy)[END_REF][START_REF] Strubelt | Inter-annual variation in species composition and richness after coppicing in a restored coppice-with-standards forest[END_REF][START_REF] Widenfalk | Plant species richness in managed boreal forests-Effects of stand succession and thinning[END_REF]. It can be expected that forest management, which influences the overall forest structure and edge physiognomy, interacts with edge influence on understorey species richness and composition by reducing the edge contrast with the surrounding matrix [START_REF] Cadenasso | Effect of edge structure on the flux of species into forest interiors[END_REF][START_REF] Harper | Edge influence on forest structure and composition in fragmented landscapes[END_REF][START_REF] Ries | Ecological responses to habitat edges: mechanisms, models, and variability explained[END_REF]. For example, open, thinned forests might be more susceptible to the penetration of wind and seeds into the forest interior due to lower edge contrast. Moreover, energy exchange in open forests is also determined by the lower albedo of dark surfaces in the understorey, such as the soil surface, litter and tree trunks, instead of green leaves of the canopy cover [START_REF] Wright | Edge microclimate of temperate woodlands as affected by adjoining land use[END_REF].

Here we studied understorey vegetation biodiversity and composition responses to distance to the forest edge in ancient forest stands (i.e. continuously forested since the first available land use maps and thus never converted to another land use type ever since) with different management types (unthinned, thinned and recently thinned) in 45 edge-to-interior forest transects across Europe. Our aim was to disentangle the effects of the management type and the distance to the edge on alpha diversity of understorey plants in multiple European regions, while accounting for environmental characteristics. We expected to find (1) higher species richness at edges, with forest specialists and generalists responding differently to management intensity and ( 2) interactions between distance-to-edge influence and forest management on richness. For example, the edge-to-interior gradient is expected to become weaker in increasingly managed (thinned) forests and thus have less impact on understory plant diversity than in recently thinned forests.

Materials and methods

Study region and set up

Our study area includes the temperate, Mediterranean and boreonemoral forest biomes of Europe. We focus on deciduous broadleaved forests, which are hotspots of biodiversity and a widespread and ecologically important forest type in most of Europe [START_REF] Brus | Statistical mapping of tree species over Europe[END_REF]. Study sites are situated in typical European fragmented landscapes, with forest patches surrounded by a matrix of arable fields and grasslands. Forest stands were selected across nine regions along a latitudinal gradient from Italy to Norway to capture much of the macroclimatic variation across Europe (Fig. 1a). In the south, middle and north of this latitudinal gradient, more specifically in Italy, Belgium and Norway, an elevational gradient with three levels (low, medium and high) was included to capture macroclimatic variation due to elevation. Study sites of the elevational gradient were selected to have a similar understorey composition as in the lowlands of the same region (with understorey plant species typical for Natura2000 habitat types 9120 and 9130). A total of 15 sites were thus selected: six regions without an elevational gradient, and three regions with three levels of elevation (Fig. 1). A study site description can be found in Appendix S1.

Forest stands were selected in a standardized way to increase comparability of sites. The stands, with a minimum forest area of 1 ha, were deciduous and mainly dominated by oaks (Quercus robur, Q. petraea, Q. cerris), Fagus sylvatica, Betula pubescens, Populus tremula, Ulmus glabra, Alnus incana and Carpinus betulus. All sites had the same type of land-use history (ancient forest) and intermediate soil moisture (mesic).

At each site, three forest stands with different forest management were selected: unthinned (1), thinned [START_REF] Tuff | A framework for integrating thermal biology into fragmentation research[END_REF] and recently thinned (3) (Fig. 1b) (see Appendix S2: site selection protocol).

1. The first type (unthinned) is typically a dense forest, with a well-developed shrub layer, high basal area (mean ± SE was here 28.8 ± 1.5 m²/ha) and high canopy cover (openness 5.8 ± 0.6 %, mean of three densitometer measurements, Appendix S3), and not recently thinned for at least 10 years and mostly > 30 years ago , indicated by e.g. the absence of cut tree stumps. 2. The second type is intermediate to types one and three (openness 6.5 ± 0.6 %, basal area 31.4 ± 1.9 m²/ha; Appendix S3); even-aged and regularly thinned, but not recently (probably 5-10 years ago). 3. The third type is typically an open, even-aged and recently thinned forest (probably within 4 years of sampling), indicated by the presence of cut tree stumps, and characterised by the absence of a shrub and subdominant tree layer, with low basal area (21.6 ± 1.3 m²/ha) and low canopy cover (mean openness 14.8 ± 2.1 %, Appendix S3).

We established 100 m transects from the southern forest edge to the forest interior of each forest stand at each site. We thus established 45 transects in total (15 sites x 3 management types). The transect started at the hypothetical line of tree stems at the edge of the forest stand (0 m). These edges are outer edges, bordering a matrix of agricultural land and were created by ancient deforestations (appendix S1). Because edge orientation affects understorey species richness and microclimatic conditions [START_REF] Didham | Edge structure determines the magnitude of changes in microclimate and vegetation structure in tropical forest fragments[END_REF]Honnay et al., 2002;[START_REF] Matlack | Microenvironment variation within and among forest edges in the eastern United States[END_REF][START_REF] Orczewska | Floristic analysis of the two woodland-meadow ecotones differing in orientation of the forest edge[END_REF], we standardized this by only locating the transects at the south-facing edge (or south-western or south-eastern), bordering with either grassland or arable land (appendix S1).

The transect was installed perpendicular to the edge and at least 100m from any another edge.

Elements that affect the microclimate at the edge or along the transect (such as water bodies, buildings and wide, tarmacked roads) were avoided.

Along each transect, five 3 by 3 m² quadrats were installed at an exponentially increasing distance from the forest edge because of the exponential change in microclimatic condition close to the edge [START_REF] Chen | Microclimate in forest ecosystem and landscape ecology -Variations in local climate can be used to monitor and compare the effects of different management regimes[END_REF][START_REF] Didham | Edge structure determines the magnitude of changes in microclimate and vegetation structure in tropical forest fragments[END_REF], resulting in a total of 225 plots. The centres of the plots were thus at distances of 1.5 m, 4.5 m, 12.5 m, 35.5 m and 99.5 m from the edge (Fig. 1c). However, in six plots the centres were located at different distances (for example at 18m instead of 12.5m), for instance, to avoid the influence of forest paths on the vegetation.

Plant biodiversity variables: understorey

Vegetation surveys were conducted at the peak of vegetation biomass from May until early July 2018, depending on the regional phenology (Appendix S1). All vascular plant species were identified and their percentage ground cover visually estimated in teams of two persons. The herb layer included all vascular species, both woody plants smaller than 1 m and non-woody plants, as well as lianas. The shrub layer was defined as all woody species with a height between 1 and 7 m and the tree layer as all trees reaching heights more than 7 m. 2019) species can shift forest guilds over the geographical gradient, and thus the forest guild a certain species was assigned dependent on the region it was observed in. Individuals only determined to the genus level were excluded from this categorisation, as well as seven species (out of 383 taxa) for which no data was available to classify them as forest specialists or generalists (Appendix S4).

Environmental predictor variables

Light transmission differences due to overstorey species identity were accounted for by means of the shade casting ability (SCA) index. The SCA index is a species-specific, expert-based index that varies from 1 to 5, indicating low to high shade casting ability of the canopy species (Verheyen et al., 2012). The SCA of the canopy (the shrub and tree layer combined) was calculated as a cover weighted mean of the scores listed in Appendix S5.

In all 225 plots, topsoil samples (0-10 cm depth, diameter 30mm) were collected for chemical analyses (pH and soil nutrient concentration) and between 10-20 cm depth for texture analysis (% silt, clay and sand). Within each plot, five random subsamples were taken after removal of the litter layer, and pooled together. The 0-10 cm soil samples were dried to constant weight at 40 °C for 48 h, ground and sieved over a 2 mm mesh. These samples were analysed for pH-H2O by shaking a 1:5 ratio soil/H2O mixture for 5 min at 300 rpm and measuring with a pH meter Orion 920A with pH electrode model Ross sure-flow 8172 BNWP, Thermo Scientific Orion, USA. To measure the total concentration of C and N in the soil, subsamples were combusted at 1200 °C and the gases were measured by a thermal conductivity detector in a CNS elemental analyzer (vario Macro Cube, Elementar, Germany). Bioavailable phosphorus (P) which is available for plants within one growing season [START_REF] Gilbert | Available soil phosphorus in semi-natural grasslands: Assessment methods and community tolerances[END_REF] was measured by extraction in NaHCO3 (POlsen; according to ISO 11263 (1994)) and colorimetric measurement according to the malachite green procedure [START_REF] Lajtha | Soil phosphorus: characterization and total element analysis[END_REF]. Total calcium (Ca), potassium (K) and magnesium (Mg) were measured by atomic absorption spectrophotometry (AA240FS, Fast Sequential AAS) after complete destruction of the soil samples with HClO4 (65 %), HNO3 (70 %) and H2SO4 (98 %) in teflon bombs for 4 h at 150 °C. Texture analysis was performed by sieving and sedimentation with a Robinson-Köhn pipette according to ISO 11277 (2009).

Landscape and topographic characteristics were derived from satellite-based global tree cover data with a spatial resolution of 30 m [START_REF] Hansen | High-Resolution Global Maps of 21st-Century Forest Cover Change[END_REF] and a pan-European digital elevation model (DEM) with a spatial resolution of 25 m, using Copernicus data and information from the European Union (EU-DEM, 2018). The proportion of forest cover, hereafter proportion forest, was assessed within a circular buffer area with a radius of 500 m and measured as the percentage of area covered by a minimum tree cover of 20 % [START_REF] Hansen | High-Resolution Global Maps of 21st-Century Forest Cover Change[END_REF]. Forest edge density was computed by transforming the forest cover mask into contour lines and summing up their length within the 500 m buffer area, using the rasterToContour function in the R package raster [START_REF] Hijmans | raster: Geographic Data Analysis and Modeling. R package version 2.8-19[END_REF].

The mean annual temperature (°C) and De Martonne's aridity index were used to take into account macroclimatic variation from latitudinal and elevational differences. De Martonne's aridity index was calculated by dividing the mean annual precipitation (mm) by mean annual temperature plus 10°C [START_REF] De Martonne | Une nouvelle function climatologique: L'indice d'aridité[END_REF]. Temperature and precipitation data were extracted for each plot from the CHELSA database for 1979-2013 [START_REF] Karger | Data Descriptor: Climatologies at high resolution for the earth's land surface areas[END_REF].

Data analysis

In total, six response variables were tested comprising measurements for alpha diversity of the understorey layer. The number of forest specialist species (1), as well as generalist richness ( 2) and total species richness (3) per plot, were calculated. In addition, Pielou's evenness index (4) [START_REF] Pielou | The measurement of diversity in different types of biological collections[END_REF] and the proportion of forest specialists compared to the total number of species (5) were also determined. As two forest guilds (1.1 and 1.2) were grouped as forest specialists, the number of closed forest species (6) was also calculated to be able to compare the closed forest species and forest specialists.

Two out of 225 plots did not contain any understorey species. These two plots were omitted for the analysis on evenness and proportion of specialists.

To take into account the hierarchical structure of the data, (generalized) linear mixed-effects models were used with transect nested within region as random effect (random intercepts).

Poisson error distributions were applied for the count data (such as species richness), while a binomial distribution was used for the proportion of forest specialists. Evenness was initially modelled with a binomial distribution as well, but this models had convergence issues, and therefore a Gaussian error distribution was applied for the model presented here. Models were built in a two-step process [START_REF] Perring | Global environmental change effects on plant community composition trajectories depend upon management legacies[END_REF]. First, a model including latitude, elevation (as continuous variables), soil variables (pH, C:N ratio, % silt, % clay, Olsen-P and K), SCA, the proportion of forest and edge density was assessed for each of the eight response variables. Mg, Ca and sand (%) were highly correlated with other explanatory variables and thus left out of the analysis (Appendix S6). Olsen-P and K were log-transformed to symmetrize the skewed distribution and eradicate influential outliers. For the sake of simplicity, quadratic terms were not included to keep the number of explanatory variables as low as possible. The most parsimonious model (single best model) was selected based on the Akaike Information Criterion (AIC) with the dredge-function of the package MuMIn [START_REF] Barton | MuMIn: multi-model inference[END_REF]. Second, the focal explanatory variables, that is, the distance to the edge and the forest management type, were added to the simplified model as a two-way interaction, resulting in the final model. Parallel to these analyses, the same procedure was followed for analyses taking macroclimatic variables into account by replacing latitude and elevation with the mean annual temperature and De Martonne's aridity index.

The distance to the edge was log-transformed to linearize the relationship. The figures presented show the back-transformed results. All continuous explanatory variables were scaled (ztransformation) to facilitate comparisons. A protocol for data exploration and model evaluation was followed precisely [START_REF] Zuur | A protocol for data exploration to avoid common statistical problems[END_REF]. All analyses were performed in R [START_REF] Core | R: A language and environment for statistical computing[END_REF], with packages lme4 [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF] and MuMIn [START_REF] Barton | MuMIn: multi-model inference[END_REF].

In addition, the magnitude of edge influence (MEI) and distance of edge influence (DEI) were calculated per management type for specialists, generalist and total species richness. For each transect, MEI was estimated for the four edge plots as (e -i)/(e + i) where e is the response value at the edge (either at 1.5 m, 4.5 m, 12.5 m or 35.5 m), and i the response value at the forest interior (99.5 m), resulting in four MEI values per transect. We report the mean and standard error of MEI calculated per management type. Next, we calculated the distance at which MEI was significantly different from zero per management type (this is the DEI), by means of a randomization test of edge influence (RTEI) with blocking [START_REF] Harper | Quantifying distance of edge influence: a comparison of methods and a new randomization method[END_REF]. A routine provided by [START_REF] Dodonov | The role of edge contrast and forest structure in edge influence: vegetation and microclimate at edges in the Brazilian cerrado[END_REF] was used: [START_REF] Harper | Edge influence on forest structure and composition in fragmented landscapes[END_REF] For each distance to the edge, the observed MEI was calculated relative to the interior plot (99.5m). ( 2) A data-set was created with the edge and interior values. (3) One value was randomly assigned to be the new edge value, and the leftover as the interior value. (4) We calculated randomized MEI based on the randomized edge and interior values. (5) Steps 3-4 were repeated 5000 times for a significance level of 1%, resulting in a distribution of randomized MEI values. ( 6) Two times the percentile of the observed MEI within the distribution of the randomized MEIs was used as the p-value (alpha = 0.01) for this distance. This routine was repeated per distance to the edge per management type. DEI was estimated as the farthest distance from the edge that was preceded by no more than one nonsignificant value [START_REF] Dodonov | The role of edge contrast and forest structure in edge influence: vegetation and microclimate at edges in the Brazilian cerrado[END_REF].

Results

In total, 351 understorey plant species were identified to the species level, and 32 only to the genus level. Of all species, 97 were forest specialists across all regions and 224 generalists, while 23 species switched forest guilds across regions (Appendix S4). Recently thinned forest stands were in general richer in generalists compared to unthinned forest stands (on average 13.3 and 11.4 species, respectively) (Fig. 2). The total number of species as well as the number of generalists decreased from the edge to the interior (Fig. 2), but for generalists this pattern was dependent on the management type: in unthinned forests, no effect of the distance to the edge was detected (Fig. 2c). However, in thinned and recently thinned forests, the number of generalist species was highest at the edge, and declined exponentially towards to forest interior. Generalist richness decreased more strongly from edge to interior in thinned than in recently thinned forests.

In thinned (and recently thinned) forests, on average 8.3 (8.3) generalist species were found at the edge, compared to only 4.5 (5.7) generalist species at the forest interior. Because of the decreasing generalist and stagnant specialist numbers from the edge to interior, the proportion of specialist species increased (Table 1). Neither the distance to the edge nor the management type influenced the number of forest specialists (Fig. 2b).

Richness of forest specialists, generalists and total species richness increased with the increasing proportion of forest in a radius of 500 m (Fig. 3a). Furthermore, the community-weighted SCA of the tree and shrub layer positively affected the proportion of forest specialists, but negatively affected total species richness and generalist richness (Fig. 3c). In contrast, soil pH was positively correlated with the total number of plant species as well as with the species richness of forest specialist (Fig. 3b). Finally, the proportion of forest specialists was the highest when the amount of edge in a radius of 500 m was minimal (Table 1). Evenness increased with percentage silt in the topsoil (Table 1). Analysing the closed forest species richness alone resulted in the same findings as the results on specialist richness (Appendix S7).The macroclimatic analysis, with mean annual temperature and De Martonne's aridity index replacing latitude and elevation, did not differ much with the initial analysis (reported in Appendix S8).

The MEI and DEI results were in line with the outcomes of the previous results. MEI and DEI varied across management types and species richness variables (Fig. 4). MEI was greatest in (recently) thinned forests for generalists and total richness, while it was close to zero for forest specialist richness (Fig. 4a). The edge influence was not significant for forest specialists in all three management types (Fig. 4b). Total species and generalist richness were significantly higher at the first 1.5-4.5 m and 18-35.5 m respectively in (recently) thinned forests. In unthinned forests, the MEI was not significantly different from zero for any of the three response variables.

Discussion

Forest management interacted with edge-to-interior gradients on generalist species richness in ancient forests along a latitudinal gradient across Europe. Contrary to our hypothesis, the edge-tointerior gradient in (generalist) species richness was stronger in thinned forests than in unthinned forests, where richness gradients were weaker. Forest edges in unthinned forests were likely more densely vegetated along the entire vertical gradient compared to those in thinned and recently thinned forest. This 'vegetation wall' may thus have served as an effective barrier that buffered the forest interior from seed inflow from the matrix [START_REF] Devlaeminck | Inflow of seeds through the forest edge: evidence from seed bank and vegetation patterns[END_REF], from wind and solar radiation (Honnay et al., 2002;[START_REF] Matlack | Microenvironment variation within and among forest edges in the eastern United States[END_REF], and the influx from warm air thereby preserving the microclimatic conditions typical of the forest interior. In this study, edge sealing might have occurred at the non-forest side of the 0 m edge plot, as the detected MEI was close to zero. In recently thinned forests, forest edge and interior microclimates did probably not contrast much due to more incoming solar radiation in the forest interior, also resulting in a more gradual gradient which levels off more quickly in species richness. For example, sparse Mediterranean temperate forests in southern Australia are found to have less pronounced edge influence on the microclimate than in closed forests, with even reversed patterns being recorded (e.g. warmer and drier forest interiors compared to edges) [START_REF] Wright | Edge microclimate of temperate woodlands as affected by adjoining land use[END_REF]. In thinned forests, however, interior and edge conditions were presumably most dissimilar, and edges relatively open, resulting in steep edge-to-core gradients in species richness penetrating 35.5 m into the forest. Hence, forest management modifies the contrast in structure and composition, which in turn determines the magnitude and distance of edge influence [START_REF] Harper | Edge influence on forest structure and composition in fragmented landscapes[END_REF].

The distance to the edge had the most significant effects on species richness. Indeed, the decline of total species richness from edge to interior corroborates the findings of [START_REF] Fraver | Vegetation responses along edge-to-interior gradients in the mixed hardwood forests of the Roanoke rives basin, North-Carolina[END_REF], [START_REF] Gehlhausen | Vegetation and microclimatic edge effects in two mixed-mesophytic forest fragments[END_REF], [START_REF] Guirado | Understorey plant species richness and composition in metropolitan forest archipelagos: effects of forest size, adjacent land use and distance to the edge[END_REF]), Honnay et al. (2002), Pellissier et al. (2013), Vallet, Beaujouan, Pithon, Roze, and[START_REF] Vallet | The effects of urban or rural landscape context and distance from the edge on native woodland plant communities[END_REF] and [START_REF] Willi | The modification of ancient woodland ground flora at arable edges[END_REF]. Higher species richness at forest edges compared to forest interiors was mostly driven by increased generalist richness at the edge, while forest specialists richness was constant along the edge-to-interior transects. Forest specialists thrive best in spatially/temporal varying intermediate light conditions and are characterized by a slow demography [START_REF] Whigham | Ecology of woodland herbs in temperate deciduous forests[END_REF], and can thus persist at the edge. At the same time, many generalists can opportunistically colonize edge habitats characterized by increased disturbance regimes [START_REF] Godefroid | Distribution pattern of the flora in a peri-urban forest: an effect of the city-forest ecotone[END_REF], warmer temperatures and a relatively higher amount of available resources, such as light and nutrient inputs [START_REF] Brunet | Understory succession in post-agricultural oak forests: Habitat fragmentation affects forest specialists and generalists differently[END_REF][START_REF] Thimonier | Floristic changes in the herb-layer vegetation of a deciduous forest in the Lorraine plain under the influence of atmospheric deposition[END_REF]. As a result, forest edges were more species-rich compared to interiors in our study (respectively 14 and 11 species per plot, on average), and the proportion of forest specialists was highest in the forest interior, and lowest at the edge. Conversely, some studies also detected decreasing richness from edge to interior for forest specialists, for instance, in oak (Q. robur) and chestnut (Castanea sativa) dominated forests in northern France [START_REF] Vallet | The effects of urban or rural landscape context and distance from the edge on native woodland plant communities[END_REF].

However, their study was performed at a smaller, regional scale and forests were not selected to be continuously forested since the last available land use maps. In addition, [START_REF] Guirado | Understorey plant species richness and composition in metropolitan forest archipelagos: effects of forest size, adjacent land use and distance to the edge[END_REF] observed higher forest species richness at the edge compared to the interior over a distance of 500 m in oak and pine dominated Mediterranean forests in Spain, indicating that distance-to-edge influence is dependent on the specific condition of the studied system.

Furthermore, forest management also influenced species diversity patterns. Recently thinned forest stands harboured overall higher generalist diversity compared to thinned and unthinned stands. In contrast, forest specialist as well as total species richness did not consistently differ between the three tested forest management. Forest management practices, such as thinning, indeed can strongly modify the understorey environmental conditions by enhancing the amount of solar radiation penetrating the canopy, increasing disturbances and modifying resource availability and microclimatic conditions (e.g. light, temperature, humidity and soil moisture) [START_REF] Ash | Changes and variability in the field layer of a coppiced woodland in Norfolk, England[END_REF][START_REF] Grayson | Understory light regimes following silvicultural treatments in central hardwood forests in Kentucky, USA[END_REF]. Changes in the amount of light along with microclimatic stability and soil disturbance are the most important drivers of understory plant diversity following forest management practices [START_REF] Ash | Changes and variability in the field layer of a coppiced woodland in Norfolk, England[END_REF][START_REF] Aude | Vegetation in Danish beech forests: the importance of soil, microclimate and management evaluated by variation partitioning[END_REF][START_REF] Brunet | Herb layer vegetation of south Swedish beech and oak forests -Effects of management and soil acidity during one decade[END_REF][START_REF] Decocq | Plant diversity in a managed temperate deciduous forest: understorey response to two silvicultural systems[END_REF][START_REF] Scolastri | Old coppice versus high forest: the impact of beech forest management on plant species diversity in central Apennines (Italy)[END_REF][START_REF] Strubelt | Inter-annual variation in species composition and richness after coppicing in a restored coppice-with-standards forest[END_REF][START_REF] Widenfalk | Plant species richness in managed boreal forests-Effects of stand succession and thinning[END_REF], which explains why light-demanding generalists can thrive in thinned, open forests. Forest specialists, on the contrary, can tolerate the changed environmental conditions. Forest specialists often lag behind environmental drivers (Hermy, 2015), requiring them to cope with altered environmental conditions [START_REF] Decocq | Plant diversity in a managed temperate deciduous forest: understorey response to two silvicultural systems[END_REF], and show extinction debts [START_REF] Vellend | Extinction debt of forest plants persists for more than a century following habitat fragmentation[END_REF]. On the other hand, most forest specialists depend on intermediate disturbance for long-term survival.

While we focussed on the effects of forest edges, forest management and their interaction, we also took into account environmental variables that potentially influence understorey species richness and composition. The amount of forest in the surrounding landscape was important to predict the number of total species, specialists and generalists. Specialist richness was associated with the forest cover. High forest cover in the landscape may positively influence dispersal and colonisation of forest specialists, as obligate forest specialists are poor colonizers [START_REF] Ehrlen | Dispersal limitation and patch occupancy in forest herbs[END_REF]. What is intriguing is that generalist richness increased similarly with the amount of forest cover in the landscape. This pattern might be associated with land use intensity, as more forest at the landscape level gives an indication about the (semi-)naturalness of the environment, and is likely negatively related to the amount of intensive agriculture and nitrogen input [START_REF] Jamoneau | Drivers of plant species assemblages in forest patches among contrasted dynamic agricultural landscapes[END_REF][START_REF] Paal | Both spatiotemporal connectivity and habitat quality limit the immigration of forest plants into wooded corridors[END_REF][START_REF] Takkis | Drivers of species richness and community integrity of small forest patches in an agricultural landscape[END_REF]. Alternatively, high forest cover in the surrounding landscape is also linked to a higher availability of forest edges (at least in fragmented landscapes), and thus perhaps to more available habitat for generalists. However, the density of forest edge in the landscape, a measure of fragmentation, did not explain additional variability in generalist richness. The density of forest edges did affect the proportion of forest specialists negatively, although this model explained little variation in response, and the numbers of generalists and specialists were not impacted significantly (Table 1). This indicates that the amount of forest in the landscape has positive effects on understorey species richness, while the density of forest edge in the surrounding landscape is of minor importance.

High community-weighted SCA (shade casting ability), a measure for light transmission differences due to overstorey composition, resulted in lower total and generalist species richness. Thus, in our study, highest species richness was detected in forest stands dominated by Betula pubescens, Quercus robur and Populus tremula, species with relatively low SCA. In contrast, understorey species richness was lowest under canopies dominated by species with relatively high SCA such as Fagus sylvatica, Carpinus betulus and Ulmus glabra. Shade-intolerant generalists were suppressed by a lack of light, as the proportion of specialist species increased slightly with increasing SCA levels. The composition of the overstorey has been shown to be important to explain compositional shifts in understorey communities, by changing the availability of light [START_REF] Baeten | Herb layer changes (1954-2000) related to the conversion of coppice-withstandards forest and soil acidification[END_REF].

Of all measured soil variables, only soil pH had a significant influence. Specialist richness increased with increasing soil pH resulting in increasing total species richness, consistent with literature [START_REF] Borchsenius | Vegetation structure and diversity of an ancient temperate deciduous forest in SW Denmark[END_REF][START_REF] De Keersmaeker | Impact of soil fertility and insolation on diversity of herbaceous woodland species colonizing afforestations in Muizen forest (Belgium)[END_REF]. However, the true relationship between pH and richness of total species and generalists is may be more hump-shaped, with an optimum around pH 5.5 (Fig. 3b) [START_REF] Schuster | Changes in species density along the soil pH gradient -Evidence from German plant communities[END_REF]. Our results are in agreement with [START_REF] Brunet | Understory succession in post-agricultural oak forests: Habitat fragmentation affects forest specialists and generalists differently[END_REF], who found that forest specialists are not affected by light (measured as the tree and shrub cover), but by soil pH, while the opposite was found for generalist species. On acidic soils, the share of forest specialists is generally lower (and thus the share of forest generalists higher) than on base-rich soils, because there is a very limited species pool of forest specialists on acidic soils [START_REF] Schmidt | Waldartenlisten der Farn-und Blütenpflanzen[END_REF].

We have shown over a large geographic gradient that thinning as forest management does not affect forest specialist species richness negatively in ancient forests, while generalists showed a contrasting pattern. Our results underpin how forest management affects edge influences on forest plant biodiversity. We recommend to take this interaction between forest management and edgeto-interior gradients better into account in future research as well as in conservation decisions because both modify the microclimate, soil processes and deposition of polluting aerosols. Therefore, this interaction is key to predict the effects of global change on forest plants in landscapes characterized by a mosaic of forest patches and agricultural land, now typical for many parts of Europe.

Figure 1 .Figure 2 .

 12 Figure 1. Study set up. Nine regions were selected across a latitudinal gradient in Europe (a). To study the effect of forest management, we compared unthinned, thinned and recently thinned forests (b). Study plots were distributed with an exponentially increasing distance from the southfacing forest edge to the interior (c)

Figure 3 .

 3 Figure 3. Understorey species richness increases with the proportion of forest cover in the landscape (a), soil pH (b) and decreasing SCA of the overstorey (c). Lines represent model estimates with the other continuous variables set at their observed mean and with forest management set to thinned. Shaded areas indicate 95 % confidence intervals. A small amount of random variation is added to the location of each point along the y-axis of panel c to avoid overplotting (jittering).

Figure 4 .

 4 Figure 4. Magnitude of edge influence (MEI) (a) and distance of edge influence (DEI) (b) per management type per response variable. In (a) the mean is depicted with error bars indicating 95% confidence intervals.

Table 1 .

 1 Effects of forest management type, distance to the forest edge, the interaction between the distance to the edge and forest management, the proportion of forest in a radius of 500 m, soil potassium concentration (K), soil pH, SCA (shade casting ability), soil texture (clay and silt percentage), elevation and forest edge density on species richness, evenness and species composition of forests across Europe. Environmental variables not shown in the table were not included in any final model (latitude, Mg, P and C:N ratio). Values here are χ² of the (generalized) linear mixed models. Symbols denote significance levels (p < 0.05 *, p < 0.01 **, p < 0.001 ***), arrows denote the direction of the continuous effects (positive ↑, negative ↓). The marginal r² is the proportion of variance explained by the fixed factors alone, while the conditional r² is the proportion of variance explained by both the fixed and random factors.

		Management	type	Distance to	edge (m)	Interaction	(management	type x distance	to edge) Proportion	forest (%)	Potassium (K)	(%)	pH		SCA		Clay (%)	Silt (%)	Elevation	(m.a.s.)	Forest edge	density (m/ha)	Variation explained Marginal Conditional r² (%) r² (%)
	Total species richness	7.96	27.05*** ↓ 5.35		5.91* ↑	0.69		4.84*	↑	11.99**	↓					22	72
	Specialist richness	1.81	0.91	0.71		5.14* ↑	2.32		6.55**	↑							16	49
	Generalist richness	10.41*	10.41*	↓ 7.35*		6.06* ↑	0.51				13.68***	↓					17	73
	Proportion specialists	7.61	20.09*** ↑ 4.45								7.27**	↑	1.12		4.9* ↑ 5.34*	↓	3	7
	Evenness	1.14	2.59	0.91				1.01							5.75* ↑			7

Acknowledgements

We are grateful to Evy Ampoorter, Shruti Bakshi, Haben Blondeel, Kim Calders, Filip Ceunen, Kris Ceunen, Robbe De Beelde, Emiel De Lombaerde, Lionel Hertzog, Dries Landuyt, Pierre Lhoir, Audrey Peiffer, Sanne Van Den Berge, Lotte Van Nevel and Mia Vedel-Sørensen for their support during the intense fieldwork campaign, to Luc Willems & Greet De bruyn for the chemical analyses and Abdulwahhab Ghrairi for the texture analysis. "Europe" icon by Denis Sazhin and "Tree" icon (adapted) by Veremeya, from thenounproject.com (Fig. 1).

Funding information

This work was supported by the Research Foundation Flanders (FWO) (project G0H1517N), and the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (ERC Starting Grant FORMICA 757833) and ERC Consolidator Grant (PASTFORWARD 614839). The FWO Scientific research network FLEUR (www.fleur.ugent.be) also supported this work.

Author contributions

PDF, SG, PVG, CM and KV designed the research; all authors collected data; SG performed statistical analyses; SG, with contributions from PDF, PVG and FZ, wrote the paper; all authors discussed the results and commented on the manuscript drafts.

Data accessibility

Raw data and R code is available at Figshare. doi: 10.6084/m9.figshare.c.4665125.v1

List of appendices