Manarshhjot Singh

Anne-Lise Gehin

Belkacem Ould Bouamama

Prognosis and Health Management using Energy Activity

Keywords: Bond Graph, Prognosis and Health Management, Energy Activity, Maintenance

Accurate detection of faults in a dynamic system is very beneficial as this information can be used in wide variety of ways by the machine operators or designers. This advantage becomes many folds when regarding the future condition i.e. time to failure, named remaining useful life, is available in addition to that of the present condition. Thus, prognosis is one of the most useful tools to improve the working of a machine as many critical decisions can be made. Prognosis can be critical for applications which risk loss of life and property. In this paper, a hybrid method, utilizing bond graph and artificial intelligence, is proposed for achieving fault diagnosis and prognosis. The Bond Graph model is used to calculate Energy Activity, which is used as a common metric for both diagnosis and prognosis. While diagnosis is achieved by using Energy Activity with neural network and short time fourier transformation, prognosis is achieved using mathematical form of Energy Activity. The proposed method is checked by simulation on a spring mass damper system undergoing a fault.

INTRODUCTION

Maintenance is one of the most important and one of the most under appreciated aspect for the proper operation of any system. In most of the industries, the maintenance strategy falls under one of the two approaches i.e. Preventive maintenance and Corrective maintenance. Under the preventive maintenance approach, maintenance action is implemented after some fixed intervals of time.

The objective in such approach is to perform maintenance action before the occurrence of the failure itself. On the other hand, in corrective maintenance approach, the maintenance action is performed once the system has failed. Both the above mentioned approaches have certain disadvantages. While there is a high monetary cost associated with the preventive maintenance due to frequent maintenance actions, high cost is also associated to corrective maintenance approach on account of the loss in working time of machine under failure. Therefore, a new Condition Based Maintenance approach is gaining popularity. Under this approach, the operating conditions of the system are continuously monitored. Whenever fault (not failure) is detected, i.e. system is in a degraded state but still in a functioning condition, further action is suggested by the monitoring system. Condition based maintenance is achieved using Prognostics and Health Management (PHM). The sequence of steps involved in PHM are shown in figure 1 and are in detail as follows [START_REF] Atamuradov | Prognostics and health management for maintenance practitioners-review, implementation and tools evaluation[END_REF]:

• Data processing: System is characterised by a set of physical values that have to follow predefined trajectories expressed by system input-output relations. Data processing extracts from the sensor signals the relevant information used for further analysis and system control. • System health estimation:System health estimation groups the FDI procedures allowing to estimate the health state of the system (faulty or non faulty). These procedures consist in three steps.

• Fault Detection: Fault detection investigates the consistency between the actual values of the system outputs provided by the sensors and the predicted values of these outputs obtained from the reference system model. A fault is detected as soon as this consistency, expressed on the form of mathematical expressions called residuals, is not respected. Equation 1gives the general form of a residual where Y measured is the actual value of the system output Y and Y estimated is its estimated value predicted by the model or other reference tables.

Residual = Y measured -Y estimated (1)
• Fault Isolation: Fault isolation consists in finding the faulty component using sensor information and, for examples, logic procedures, signal processing or reference tables. • Diagnosis: Diagnosis gives an interpretation of the nature and the cause of the fault. • Prognosis: Prognosis is a dynamic estimate of the degradation of the system. This deals with calculation of the End of Life of a system, a point in time at which the fault increases to its maximum limit resulting in system failure. Remaining Useful Life (RUL) of the system is expressed by equation 2 where t f ailure is the predicted time where the system cannot continue to operate due to complete failure and t current is the time at which the RUL is calculated.

RU L = t f ailure -t current
(2) This RUL is represented on figure 2 • Decision making: Detecting the occurring fault and estimating the RUL of the system can help in both protecting the system components, the system environment and/or ensuring the continuity of service when possible. Decision making can range from immediate human intervention to implementation of fault tolerant control by putting in priority users safety measures, system protection and continuity of service.

Prognosis is the most important step of PHM and has attracted the attention from researchers all over the world. The quick and accurate prediction of RUL is important because the subsequent action after the appearance of fault depends on the RUL. The various approaches introduced for prognosis fall under one of the following:

(1) Data-based approach (2) Model-based approach (3) Hybrid approach For data driven approach [START_REF] Sutharssan | Prognostic and health management for engineering systems: a review of the data-driven approach and algorithms[END_REF], the failure history from previous failed operation condition is utilised for calculation of the RUL. Various mathematical and statistical techniquesSi et al. ( 2011) can be used achieve prognosis. Neural Networks [START_REF] Byington | Data-driven neural network methodology to remaining life predictions for aircraft actuator components[END_REF], Fuzzy Logics, Principal Component Analysis [START_REF] Zhang | Machinery condition prognosis using multivariate analysis[END_REF], Hidden Markov model [START_REF] Peng | A prognosis method using age-dependent hidden semi-markov model for equipment health prediction[END_REF] etc or their combination [START_REF] Jahromi | Sequential fuzzy clustering based dynamic fuzzy neural network for fault diagnosis and prognosis[END_REF] can also used to improve results. However, for certain applications, especially for critical systems like space shuttles, nuclear plants etc., it is not easy to have a large database of failure information as every system failure can be accompanied by huge loss in life and property. In such systems, it is always preferred to use model-based approach. In modelbased approach the system behavior is known with accurate mathematical equations. The mathematical model can be used with various mathematical tools of ways like dissociation [START_REF] Chelidze | A dynamical systems approach to failure prognosis[END_REF]Cusumano (2004), particle fil-terJha et al. (2016), similarity [START_REF] Roychoudhury | An integrated model-based diagnostic and prognostic framework[END_REF], neural network [START_REF] Singh | Bond graph model for prognosis and health management of mechatronic systems based on energy activity[END_REF]. The availability of mathematical equations can allow faster diagnosis and prognosis without the need of a failure database. The accurate model parameters are seldom known for physical systems, hence hybrid approaches are commonly employed for prognosis. Hybrid [START_REF] Neerukatti | Fatigue life prediction using hybrid prognosis for structural health monitoring[END_REF] approaches are a combination of data-based and modelbased approaches. The voids observed in the model-based techniques are filled using data, and hence a combination model is used for prognosis.

Modern physical systems contain subsystems of different domains like electrical, rotational, hydraulic etc. The management of interaction between such subsystems is a tedious task. Energy being the common currency of exchange between sub-systems, makes it an ideal parameter for study and analyses of complex systems. Energy has been studied in various forms like Analytical Redundancy Relations [START_REF] Bouamama | Derivation of constraint relations from bond graph models for fault detection and isolation[END_REF], Energy Activity In-dexSingh et al. (2018), Junction Activity [START_REF] Jha | New concept of junction activity in a bond graph model: Application for fault identification[END_REF]), Energy Efficiency Hoang et al. (2017) for various stages of Prognosis and Health Management. However, an energy based framework for both diagnosis and prognosis can be of great interest. In this paper, such a framework is developed using Energy Activity. Energy Activity is chosen because the increasing nature of energy activity has potential to capture minute faults over time. The occurrence of the mathematical singularities due to inherent nature of Energy Activity are avoided by using Neural Networks instead of mathematical evaluation. Energy activity is well developed for application using bond graph. Therefore, the current framework is developed using bond graph model. Also, bond graph has been used previously for integrated diagnosis and prognosis frameworks Jha (2015) [START_REF] Yu | Incipient fault diagnosis and prognosis for hybrid systems with unknown mode changes[END_REF]. Hence, bond graph is a proven tool for such framework.

This paper is organised as follows. Section 1 gives a basic introduction to Prognosis and Health Management and it's various processes. In Section 2 the concept of Pseudo Energy Activity is introduced and developed in the context of bond graph modelling. Section 3 gives the detailed methodology of the using pseudo energy activity for for prognosis and health management. In Section 4, the proposed method is simulated using a spring mass damper system and the simulation results are presented. The paper is concluded in Section 5.

PSEUDO ENERGY ACTIVITY FOR PROGNOSIS

Bond Graph

Bond Graph [START_REF] Mukherjee | Bond graph in modeling, simulation and fault identification[END_REF] is a technique used to model the behavior of dynamic systems. Bond Graph is a component based technique therefore the individual components and subsystems can be easily visualised while modelling the dynamics. This also gives a better understanding of the results. The bond graph technique is based on the law of conservation of energy and therefore used power, i.e. time derivative of energy, to capture the dynamics of the system. As energy is the common currency of exchange among the different domains, therefore bond graph technique is suited for modelling of multi-physics system with a combination of different domains of energy like mechanical, electrical, thermal etc.

In bond graph technique power in every component is measured as a product of an agent bringing about the change, called generalised effort e, and the rate of observed change, called generalised flow f . The components in the system can be either energy storing elements or energy dissipators. Components storing generalised kinetic energy are represented as I elements, components storing generalised potential energy are represented as C elements and the energy dissipators are represented as R elements.

The system can also have energy conversion elements i.e. Transformers T F and Gyrators GY . The different components can interact with each other and the system in general, through either a 1 or 0 junction. Components sharing a common 1 junction have same generalised flow whereas those sharing common 0 junction have a common generalised effort. External input to the system can be either an effort source SE or a flow source SF . As discussed above, it is very easy to manage individual components when using bond graphs, therefore use of bond graphs is beneficial for the fault diagnosis. Also, the 0 and 1 junctions impose a structural condition on the behavior of on the various components of the system. This condition always represents a physical phenomenon like kirchhoff's law or physical condition like physical junction etc. Hence, the structural nature of bond graph modelling has been successfully used for fault diagnosis [START_REF] Bouamama | Derivation of constraint relations from bond graph models for fault detection and isolation[END_REF].

Energy Activity

All dynamic components, in any dynamic system, irrespective of the domain interact with energy. The components can either absorb energy during certain phase and release the absorbed energy during the other. Some components just dissipate energy. The amount of energy that they interact with depends on the component values of the system. When a system is under fault, the component value, and hence the energy interaction for that of the system component changes. As energy can not be created or destroyed, the whole system experiences a redistribution of energy as a consequence of the fault. Also, for a system to work properly, every component should play it's role in energy interaction within certain limits. Monitoring this energy interaction in the various components of the system can help in PHM. Bond graph being a component based approach is suited for study of such interaction.

This section deals with the introduction to the known concept of Energy Activity and Energy Activity Index. Furthermore the modified pseudo-Energy Activity and pseudo-Energy Activity Index are introduced for utilization in prognosis and health management.

Energy Activity and Energy Activity Index Energy Activity was introduced in Louca et al. ( 2010) as a tool for physical model reduction. EA is the total energy interaction that a component has with the system over a time period. The major difference between energy and energy activity is that while the energy associated with a component can increase or decrease, the energy activity of the component always increases with time. This makes energy activity a more suited component for prognosis than energy. The equations for energy and energy activity is given by equation 3 and eqution 4 respectively.

Energy = a+∆t a (P ower)dt (3) EA = a+∆t a |P ower| dt (4)
Energy Activity Index of a component during a time interval is the ratio of the energy activity of that component to the total energy activity of the system during this time. Therefore, energy activity index analysis provides the relative picture of the activity of different components on the system. The expression for calculation of energy activity index of component i in a system with n components is given by equation 5.

EAI i = EA i n i=1 EA i (5)
Pseudo Energy Activity and Pseudo Energy Activity Index In the context of bond graph, EA is given by equation 6 and 7. Here e and f are the generalized effort and flow respectively.

EA =

e = φ R f or f = φ -1 R e (8)
Consequently, equation 6 and equation 7 can be written as

EA = a+∆t a |f.φ R f | dt = φ R a+∆t a S 2 dt (9)
For bond graph element C and I, g is an integral or differential relation since:

e(t) = φ C t 0 f (t)dt f (t) = d dt (φ -1 C e(t)) f (t) = φ I t 0 e(t)dt e(t) = d dt (φ -1 I f (t)) (10)
When a component in a system undergoes a fault, physical law imposed by the component is modified since the parameter value is modified (as example, more energy than expected is dissipated by R element or less energy than expected is stored by C element). Parameter modifications due to component faults are not known in advance. This is the aims of the PHM i.e. to detect, locate and quantify them. In this case, the actual Energy Activity of the component is different of its predicted one. To compare the predicted EA to the actual one, the concept of pseudo-Energy Activity (pEA) is introduced. For calculation of pseudo-Energy Activity, the ideal parameter value θ is used instead of the actual parameter value φ. The equation defining pseudo-Energy activity is given by equation 11. Similarly the Energy Activity Index calculated using pseudo-Energy Activity is called pseudo-Energy Activity Index. The expression for calculation of pseudo energy activity index of i th component in a system with n components is given by equation 12.

pEAI i = pEA i n i=1 pEA i (12)
Hereafter, in this paper, unless otherwise stated Energy Activity and Energy Activity Index refer to pseudo-Energy Activity and pseudo-Energy Activity Index.

Residuals generated from Energy Activity Indexes have inherent characteristics which make them ideal for fault detection or isolation in time domain [START_REF] Singh | Bond graph model for prognosis and health management of mechatronic systems based on energy activity[END_REF]. However, the mathematical evaluation of the discussed characteristics is not possible due to singularities in the solution. Hence, a neural network approach is suggested.

PHM PROCESS USING ENERGY ACTIVITY

A general overview of PHM process using Energy Activity is shown in figure 3. A virtual system in fault-free conditions is simulated in parallel with the actual system. Similar outputs from both the systems are used for the remainder of the process. The outputs from the virtual system are used to calculate the reference Energy Activities. However, for calculation of the Energy Activities associated with real components, the output signals must be de-noised. The energy activity indexes of various components corresponding to the actual and virtual system are calculated and subsequently residuals are generated. The process of filtering can modify signal phase & magnitude. This introduces error in the residuals. These errors can be easily removed by analysing the residuals in frequency domain. As frequency of residual signal changes when system moves from a fault-free state to faulty state, residuals are analysed using Short Time Fourier Transformation technique. The results from Short Time Fourier transformation are compared with a pretrained Neural Network to achieve fault isolation. Once fault isolation is performed, fault prognosis can be performed using mathematical equations. The equations are discussed in detail in the following sections.

Neural Network

Neural network is trained to use the residuals from Energy Activity Indexes to classify the fault location. The neural network must be trained before the actual application, using the residuals generated from simulation as training data. The process of generation of Neural Network training data is shown in figure 4.

A fault is simulated by changing the value of a component parameter. Therefore to have an exhaustive database, a range of fault for every component is decided. A faultfree system is simulated in parallel with the faulty system. The generated residual is the difference of energy activity indexes from faulty and fault-free system. Fourier transformation is performed on the residuals to have a frequency picture of the residual. The location of the peaks is recorded in the frequency domain. The noisy peaks i.e. peaks at high frequency and low amplitude are removed.

The remaining peaks are added to the neural network training data.

The neural network is trained under the following conditions

• At a time only one component can have fault.

• Fault is introduced at the start of the simulation.

• For every iteration the fault magnitude remains the same.

Fault Isolation

As discussed earlier, for fault isolation, Fourier transformation is performed on the residuals. However, in a real system, fault can occurs after some period of fault-free operation. Hence, the residuals change from zero to non-zero after some time. In such a case, a direct fourier transformation does not properly capture the change in behavior of residuals. Therefore, short time fourier transformation is used. In short time fourier transformation [START_REF] Liu | Rolling bearing fault diagnosis based on stft-deep learning and sound signals[END_REF], the residual signal is divided into small windows of equal duration, and subsequently fourier transformation is applied on it. Using the short time fourier transformation, The window for short time fourier transformation should be equal to the simulation time used for training the neural network. This assures that the neural network is able to recognise the pattern properly. The window of short time fourier transformation should also be more than the time ∆t used in equation 11 to calculate the energy activity.

Fault Prognosis

Once fault isolation is achieved, the mathematical form of Energy Activity can be used to calculate the real variation in the fault parameter, and furthermore the remaining useful life, assuming that the allowable limits of a parameter are known for failure.

In order to estimate the dynamics of the degradation (i.e. the time variation in the value of parameter), the time derivative of the EA is required.

As the fault (degradation) is due to a modification of the value of the component parameters in time, the Energy Activity for a faulty component can be expressed as a function of the component input signals, which itself is a function of the all the component parameters φ (considered as time varying) and time t. (See equation 13)

EA = f (S(φ), t) (13)
From equation 13 the following can be calculated For an R-element, g is algebraic and from equation 9

EA = θ R b a S 2 dt (17)
Therefore,

∂EA ∂S = θ R b a ∂S 2 ∂S dt = 2θ b a Sdt (18)
For C and I elements, the calculation of ∂EA ∂S is more difficult due to the forms of expressions in equation 10.

Equation 16 can now be integrated as shown in equation 19 to have an estimation of fault parameter from sensor data.

φ = dφ dt dt + φ ideal (19)
The continuous calculation can then be extrapolated according to a know degradation trend. If a degradation trend is unknown, then a polynomial equation can be used to extrapolate the component value. The Remaining Useful Life can be easily calculated if the safe limits of the component values are known beforehand (figure 2). The trend of component degradation is extrapolated to find the point in time when the component value reaches the allowable limit. This point is called the End of Life.

The time difference between present and end of life is the Remaining Useful Life.

APPLICATION

System

In order to check the proposed methodology, a simulation is performed using a simple spring-mass-damper system.

The system is shown in figure 5. The pre-requisites i.e. the dynamic model using bond graph, and the ideal & and safe working limits of components are given by figure 6, and table 1 respectively.

Fault Isolation

During the simulation, a fault condition is indicated by deviation in the component value from ideal. For the purpose of the simulation, the fault magnitude is modelled as change in component value. Fault is introduced in the spring. The variation in spring stiffness is shown in figure 9. A fault is introduced at 5 seconds which continues to decrease the spring stiffness. At 50s from the start of simulation, a corrective action is simulated and spring stiffness starts to increase to recover its initial value at 100s. The Time-Frequency map obtained from the Short Time Fourier Transformation is applied on the obtained residual. The Time-Frequency map is shown in figure 8. The data entries corresponding to the each time interval are given as input to the neural network trained in the previous step. The neural network is able to correctly predict the fault location as spring.

Fault Prognosis

The equation 16 is used for evaluating the spring stiffness change rate. The change in spring stiffness φ introduces a change in the damper input S, which affects the Energy Activity. The calculated stiffness change rate is passed through a median filter in order to remove sharp peaks due to numerical anomalies. The stiffness change rate after filtering is shown in figure 10. This change rate can be integrated to find the actual spring stiffness. The calculated variation of spring stiffness is shown in figure is achieved by using a combination of Neural Network and Short Time Fourier Transformation. Given that the dynamic model of the system is known, the neural network is trained using fault simulations and does not require failure data. The prognosis process is completed using the mathematical nature of Energy Activity for energy dissipators. This can also be a limitation for the proposed process as the prognosis process can not utilize the energy storing elements. The proposed method is simulated for finding the end of life of a spring mass damper system undergoing a fault. The method is able to predict the fault location correctly and recreate the parameter values of the component under fault with good accuracy.

Fig. 1 .

 1 Fig. 1. Prognosis and Health Management process.

Fig. 2 .

 2 Fig. 2. Extrapolation of Component parameter for Remaining Useful Life

S

  is the signal input from the system to the component. According to the causality (derivative or integral) associated with the component, S represents the flow or effort input signal. Function g in equation 7 expresses the relation between the component input and output. It corresponds to the physical law applied by the component. It uses the component parameter value φ. For the bond graph element R, g is algebraic

Fig. 3 .

 3 Fig. 3. Generalised procedure.

  derivative of the Energy Activity calculated from the real system. • ∂EA ∂t represents the variation of the Energy Activity only in time, i.e. due to no change in φ. This term can be calculated as the time derivative of the Energy Activity of the fault free system. • ∂S ∂φ represents the variation of the input signal of the component due to the modification of component parameter value. As the dynamic model of the component is known as a pre-requisite this value can be calculated easily. • ∂EA ∂S depends on the nature of relation g in the equation 7.

Fig. 5 .

 5 Fig. 5. Spring Mass Damper System.

Fig. 6 .

 6 Fig. 6. Bond Graph of Spring Mass Damper System.4.2 Neural Network TrainingThe first step is to pre-train a neural network. The procedure explained in figure4is used to train the neural network. The fault range for training the neural network are those given in table 1. For creating the training dataset, the fault range for every component is divided into 20 equal intervals. The residual defined in equation 1 is calculated at the I element i.e. mass. The frequencyamplitude graph used for training of neural network is shown in figure ??. The x-axis represents the frequency in Hz and the y-axis represents the amplitude. The cross marks of different colors represent the peaks observed at different magnitudes of faults. From the figures it is evident that peaks locations can provide a unique fault signature to the components and can be easily used to train a neural network. A default pattern recognition/classification algorithm provided in MATLAB is used with 5 hidden layers. 70% of the available database is used for training while 15% of database is used for validation and testing each.

Fig. 7 .

 7 Fig. 7. Neural Network Dataset.

Fig. 8 .

 8 Fig. 8. Time-Frequency map of Short Time Fourier Transformation.

Fig. 10 .

 10 Fig. 10. Calculated Parameter variation Rate.11. The error in the calculated spring stiffness is shown in figure12. From the figure it is evident that the spring stiffness is calculated with good accuracy.At any time when the fault is observed the trend of the parameter variation can be extrapolated using a polynomial equation. For the current example a first order polynomial is used. The point of failure a.k.a. End of Life is reached when the extrapolation trend reaches the allowed limit of the component value. The Remaining Useful Life is continuously monitored. Once the corrective action is applied the calculation of Remaining Useful Life is continued. This represents the amount of time for which the corrective action can be applied before the component value overshoots the allowable limits. Calculation of End of Life is shown in figure11.

Fig. 11 .

 11 Fig. 11. Calculation of End of Life.

Fig. 12 .

 12 Fig. 12. Error in calculated spring stiffness.

Table 1 .

 1 Ideal value of system components

	Element	Value in SI Lower limit Upper limit
	Force	10
	Spring stiffness	100	85	105
	Mass	10	9	11
	Damping coefficient	0.5	0.4	0.6

Table 2 .

 2 Calculated End of LifeIn this paper a model based method for Prognosis and Health Management is proposed using Energy Activity. Both the diagnosis and prognosis processes are completed using variants of Energy Activity as a metric. Diagnosis

	Cause of variation	End of Life Time
	Fault	55s
	Corrective action	113s
	5. CONCLUSION


Guy Fayolle
email: guy.fayolle@inria.fr.

Roudol Iasnogorodski
email: iasnogorodski@mail.ru

Conditions for some non stationary random walks in the quarter plane to be singular or of genus 0

Keywords: Algebraic curve, functional equation, generating function, genus, quarterplane, Riemann surface, singular random walk AMS 2000 Subject Classification: Primary 60G50; secondary 30F10, 30D05

, which considers only the special situation where t ∈]0, 1[ is a transcendental number over the field Q(p i,j ). In addition, when p 0,0 = 0, our classification holds for all t ∈]0, +∞].

Historical introduction

Discrete two-dimensional random walks in domains with non-smooth boundaries are of interest to several groups in the mathematical community. In fact, these processes are encountered in pure probabilistic problems, as well as in applications involving queueing theory or enumerative combinatorics. Some four decades ago, original mathematical methods have been proposed to determine the invariant measure of such Markov processes, say (M t , N t ) ∈ Z 2 + at time t, evolving in the positive quadrant, see the book [START_REF] Fayolle | Random Walks in the Quarter Plane: Algebraic Methods, Boundary Value Problems[END_REF]. When the jumps in the the interior of the quadrant are of size one, the bi-variate generating function Q(x, y) of the probability distribution satisfies a functional equation of two complex variables, whose general solution, obtained via the so-called kernel method, amounts to solving a boundary value problem on a Riemann surface of genus 0 or 1, except for some degenerate walks, referred to as singular. The complete classification was obtained in [START_REF] Fayolle | Random Walks in the Quarter Plane: Algebraic Methods, Boundary Value Problems[END_REF]Section 2.3].

It turns out that the method can also be applied to analyze the transient behaviour of these walks, after switching to the Laplace transforms

F (x, y, s) = ∞ 0 e -st E[x Mt y Nt ]dt, (s) ≥ 0,
in which case F (x, y, s) satisfies a functional equation, where s merely plays the role of a parameter.

In this article, we will consider the problem of counting walks in the quarter plane.

The related counting generating function (CGF) is of a transient nature and satisfies the functional equation (2.2), in which the time parameter t is in a certain sense the analogous of s for the Laplace transforms. In this context we will provide the classification of the walks, which, as expected, turns out to be almost identical to that mentioned above.

Model and notation

Enumeration of planar lattice walks has become a classical topic in combinatorics.

For given allowed jumps (or steps), it is a matter of counting the number of paths starting from some point and ending at some arbitrary point in a given time.

Here we assume the random walks evolve in the positive quarter plane, and the jumps have size 1. Hence, the following transitions can take place.

         (m, n) → (m + i, n + j), with probability p i,j , |i|, |j| ≤ 1, ∀(m, n) > (0, 0); (0, n) → (i, n + j), with probability p i,j , 0 ≤ i ≤ 1, |j| ≤ 1, ∀n > 0; (m, 0) → (m + i, j), with probability p i,j , 0 ≤ j ≤ 1, |i| ≤ 1, ∀m > 0; (0, 0) → (i, j), with probability p i,j , 0 ≤ i, j ≤ 1.
Without restricting the generality, we suppose

(i,j)∈S p i,j = 1.
It will be convenient to denote by S the set of the admissible jumps, which are included in the set of the eight nearest neighbors, so that

S ⊂ {-1, 0, 1} 2 \ {(0, 0)}.
Let f (m, n, k) denote the probability for a walk in Z 2 + of reaching the point (m, n) after k steps, starting from (0, 0). Then the corresponding CGF

F (x, y, t) = m,n,k 0 f (m, n, k)x m y n t k (2.1)
satisfies the following functional equation (see [START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF] for the details) K(x, y, t)F (x, y, t) = K(x, 0, t)F (x, 0, t) + K(0, y, t)F (0, y, t) + l(x, y, t), (

where

           K(x, y, t) = xy tS(x, y) -1 , S(x, y) = (i,j)∈S p i,j x i y j , l(x, y, t) = -tp -1,-1 F (0, 0, t) -xy.
(2.3)

F (x,
p ij x i y j -1 (2.5) def = t[a(x)y 2 + b(x, t)y + c(x)] def = t[ a(y)x 2 + b(y, t)x + c(y)], (2.6) with 
                a(x) = 2 i=0 p i-1,1 x i , b(x, t) = p -1,0 + (p 0,0 -1/t)x + p 1,0 x 2 , c(x) = 2 i=0 p i-1,-1 x i , (2.7)
and similar equations for a(y), b(y), c(y).

Let us introduce the discriminant of (2.6) , always viewed as a polynomial in x D(x, t)

def = b 2 (x, t) -4a(x)c(x). (2.8)
Then it will be convenient to write

D(x, t) = 4 i=0 d i (t)x i , (2.9)
where, by using (2.7),

               d 0 (t) = p 2 -1,0 -4p -1,1 p -1,-1 , d 1 (t) = 2[p -1,0 (p 0,0 -1/t) -2(p -1,1 p 0,-1 + p 0,1 p -1,-1 )], d 2 (t) = (p 0,0 -1/t) 2 + 2p 1,0 p -1,0 -4(p 1,1 p -1,-1 + p 1,-1 p -1,1 + p 0,1 p 0,-1 ), d 3 (t) = 2p 1,0 (p 0,0 -1/t) -4(p 1,1 p 0,-1 + p 0,1 p 1,-1 ), d 4 (t) = p 2 1,0 -4p 1,1 p 1,-1 , (2.10)
noting that d 0 and d 4 are in fact independant of t. It is immediate to check that d 2 (t) is always strictly positive. Indeed we have the inequality

d 2 (t) ≥ (1 -1/t) 2 + (p 1,1 -p -1,-1 ) 2 + (p 1,-1 -p -1,1 ) 2 + (p 0,1 -p 0,-1 ) 2 .
Hence, the degree of D(x, t) in x is always 2.

We shall need the following simple notion of positivity for an arbitrary polynomial.

Definition 2.1. A polynomial p(x) is said to be positive (resp. negative) to mean that its coefficients are all real positive (resp. negative) and not simultaneously 0.

3 The roots of the discriminant D(x, t)

The main result of this section resides in he following theorem.

Theorem 3.1. For all t ∈]0, 1], D(x, t) has 4 real roots, exactly two of them being inside the unit disk, and double roots can occur only at x = 0 or x = ∞.

Proof. For all x satisfying a(x)c(x) ≥ 0, we introduce

φ 1 (x, t) def = b(x, t) + 2 a(x)c(x), φ 2 (x, t) def = b(x, t) -2 a(x)c(x), (3.1)
whence, according to (2.8),

D(x, t) = φ 1 (x, t)φ 2 (x, t).
Keeping in mind that a(x) and c(x) are positive (see Definition 2.1), they have no positive real roots. Our goal is to distribute the 4 roots of D(x, t) between φ 1 and φ 2 .

The roots of φ 1 (x, t) and φ 2 (x, t). The next inequalities are immediate

φ 2 (1, t) φ 1 (1, t) 1 -1/t < 0, φ 1 (0, t) 0. (3.2)
Let ξ 0 denote the largest real root of a(x)c(x), when it exists. Then ξ 0 0 and φ 1 (ξ 0 , t) = φ 2 (ξ 0 , t) = b(ξ 0 , t) 0. Consequently, φ 1 (x, t) and φ 2 (x, t) each have a root on ]ξ 0 , 1[. On the other hand, if ξ 0 -1, then a(x)c(x) 0 for x ∈ [-1, 1] and

φ 1 (-1, t) φ 2 (-1, t) > 0,
which, in tandem with (3.2), proves indeed that φ 2 (x, t) has always a root on ] -1, 1[. This conclusion still holds when ξ 0 does not exist, since then a(x)c(x) is positive and φ 1 (x, t), φ 2 (x, t) are real ∀x ∈ R.

So, when degree of D(x, t) = 2, there are two real roots inside the unit disk, harmoniously distributed between φ 1 (x, t) and φ 2 (x, t).

In addition, we note that, when degree of D(x, t) 3,

lim x→+∞ φ 1 (x, t) = +∞, (3.3)
in which case φ 1 (x, t) has a root on ]1, ∞[ and a root on [0, 1[.

Degree of D(x, t) = 3, i.e., d 4 (t) = p 2 1,0 -4p 1,1 p 1,-1 = 0.
Then two roots of D(x, t) come from φ 1 (x, t), and the third one from φ 2 (x, t).

a(x)c(x 0, ∀x ∈] -∞, ξ 1 ], φ 2 (ξ 1 , t) = b(ξ 1 , t) > 0,
and, since lim

x→-∞ φ 2 (x, t) = -∞, it follows that φ 2 (x, t) has a negative real root on ] -∞, ξ 1 [. If ξ 1 does not exist, then φ 2 (x, t) has a root on ] -∞, -1].
In conclusion, the distribution of the roots of D(x, t) between φ 1 (x, t) and φ 2 (x, t) implies that a double root exists if, and only if,

φ 1 (x, t) = φ 2 (x, t) = 0 ⇐⇒ b(x, t) = a(x)c(x) = 0. (3.4)
But the roots of b(x, t) are real positive, while those of a(x) and c(x) have negative real parts. Hence (3.4) can only take place if x = 0 or x = ∞.

The proof of the theorem is terminated.

4 Classification of the singular random walks Definition 4.1. A random walk is called singular (see [START_REF] Fayolle | Random Walks in the Quarter Plane: Algebraic Methods, Boundary Value Problems[END_REF]) if the associated polynomial K(x, y, t) is either reducible or of degree 1 in at least one of the variables.

We establish a useful lemma, which is of an algebraic nature and gives conditions for the factorization of the kernel K(x, y, t).

Lemma 4.2. Let A[x] be the algebra of polynomials in x with coefficients in an arbitrary field A containing the rational numbers Q. For given a, b, c

∈ A[x] satisfying b 2 -4ac = kp 2 , (4.1)
where p is a polynomial and k ∈ A, there exist α, β, f 1 , f 2 ∈ B[x], such that

ay 2 + by + c = (f 1 y -α)(f 2 y -β), (4.2)
where B def = A + √ k denotes the field A to which is added the element √ k. In addition, the relation (4.1) is necessary to have the factorization (4.2).

Proof. Let Z 1 , Z 2 be the roots of

Z 2 + bZ + ac = 0, (4.3)
so that, with an obvious notation,

Z 1,2 = - b 2 ± p 2 √ k.
In B[x], let f 1 be the greatest common divisor (g.c.d) of Z 1 and a, so that

Z 1 = βf 1 and a = f 1 f 2 ,
where β and f 2 are relatively prime. Then, from Z 1 Z 2 = ac, we get βZ 2 = cf 2 .

Hence, the g.c.d. of β and f 2 being the unit element, f 2 divides Z 2 and we shall put

Z 2 = αf 2 .
Setting for a while Z = ay, the announced factorization (4.2) follows directly from the chain of equalities

Z 2 +bZ +ac = (Z -Z 1 )(Z -Z 2 ) = (ay-βf 1 )(ay-αf 2 ) = a(f 2 y-β)(f 1 y-α) (4.4)
together with the identity

Z 2 + bZ + ac = a(ay 2 + by + c).
As for the necessity of (4.1) to have (4.2), we therefore assume

ay 2 + by + c = (f 1 y -α)(f 2 y -β).
Letting Z

1 def = αf 2 , Z 2 def
= βf 1 , one sees that Z 1 , Z 2 are the respective roots of (4.3), which both belong to ∈ B[x] and satisfy

(Z 1 -Z 2 ) 2 = (Z 1 + Z 2 ) 2 -4Z 1 Z 2 = b 2 -4ac. (4.5)
On the other hand, as part of the hypothesis, we can set

Z i def = u i + v i √ k ∈ B[x], for i = 1, 2, where u i , v i ∈ A[x].
Then, using the fact that a, b, c belong to A[x], together with the relations αf 2 + βf 1 = -b, and αβf 1 f 2 = ac, we get immediately

u 1 = u 2 = -b/2, v 1 + v 2 = 0, Hence (Z 1 -Z 2 ) 2 = [2v 1 √ k] 2 = 4kv 2 ,
which, comparing with (4.5), yields exactly (4.1).

The proof of the lemma is terminated.

Proof. Assume k > 0 in Lemma 4.2. Then √ k ∈ R and B = A = R.
The classification of the family of random walks under study is given by the next theorem.

Theorem 4.5. For t ∈]0, 1] and p i,j ∈ [0, 1], p 0,0 = 1, the random walk is singular if, and only if, one of the following conditions holds:

(i) There exists (i, j) ∈ Z 2 , |i| 1, |j| 1, such that only p i,j and p -i,-j are different from 0 (see figure 4.1a and the three cases obtained by rotation); From now on, we shall assume ac ≡ 0.

Case K(x, y, t) = u(x, t)v(x, y, t). Here u (resp. v) stands for a polynomial in x (resp. (x, y). It is immediate to see that the degree of u(x, t) in x is 1. But this means that a, b, c have a common real root in x, which can occur only at x = 0. Indeed, by (2.7), it appears that b(x, t) has always two positive real roots, while the roots of a and c cannot be positive. Thus the only possibility is p -1,1 = p -1,0 = p -1,-1 = 0, or equivalently c ≡ 0, which corresponds to the walks of type b in figure (4.1). Hence, c(y

) ≡ 0 ⇐⇒ K(x, y, t) = tx[x a(y) + b(y, t)].
Similarly, by exchanging the roles of x and y, we obtain the walks of type d in figure (4.1). In particular,

c(x) ≡ 0 ⇐⇒ K(x, y, t) = ty[ya(x) + b(x, t)].
Under Definition 4.1, all singular random walks satisfy (4.2) together with the system

αf 2 + βf 1 = -b, (4.6)
f 1 f 2 = a, (4.7
)

αβ = c, (4.8)
where α, β, f 1 , f 2 are polynomials of degree at most 2 in x, with real coefficients by Corollary 4.4. In the course of the proof, we shall omit the variables (x, t) for the sake of readability.

When a and c are positive (see Definition 2.1), (4.7) and (4.8) imply that (f 1 , f 2 ) are both either positive or negative, and likwise for (α, β). This is true because a and c are of degree at most 2, and we shall choose (f 1 , f 2 ) positive.

Recalling that b ≡ b(x, t) = p -1,0 + (p 0,0 -1/t)x + p 1,0 x 2 , the identity (4.6) leads to consider two different situations.

(i) p 1,0 + p -1,0 = 0. Here it is impossible to satisfy (4.6) for t ∈]0, 1], since the three coefficients of the polynomial b do not have the same sign.

(ii) p 1,0 = p -1,0 = 0. Then (4.6) can hold only if (α, β) are positive and if (αf 2 , βf 1 ) are proportional to x. Upon combining (4.7)-(4.8), we obtain

(αf 2 )(βf 1 ) = ac = Kx 2 ,
where K is a constant. Thus the coefficients of x 4 , x 3 , x, x 0 of the polynomial a(x)c(x) must cancel out, yielding from (2.7) the respective equations

         p 1,1 p 1,-1 = 0, p 1,1 p 0,-1 + p 1,-1 p 0,1 = 0, p 0,1 p -1,-1 + p -1,1 p 0,-1 = 0, p -1,1 p 0,-1 = 0. (4.9)
Now it is easy to check that the puzzle-like system (4.9) generates exactly the 4 following random walks: By exchanging the roles of x and y, we also obtain the singular walk where only p 1,0 and p -1,0 are different from zero.

         p 1,-1 = p 0,-1 = p -1,1 = p 0,1 = 0, figure (
The proof of the theorem is terminated.

Remark 4.6. When p 0,0 = 0, one verifies easily that Theorem 4.5 holds in fact for all t ∈]0, +∞[.

Kernel of genus 0

It is wellknown that, for all non-singular random walks, the Riemann surface S defined by the algebraic curve (2.4) has genus 0 or 1. Morever, the genus of S is 0 if, and only if, the discriminant D(x, t) defined in (2.8) has a multiple zero in x, possibly infinite.

For the sake of completeness, we also introduce the branch points x i (resp. y i ), i = 1, . . . , 4, of the algebraic function Y (x, t) satisfying K(x, Y (x, t), t) = 0 [resp. X(y, t) satisfying K(X(y, t), y, t) = 0]. These branch points depend of course on t.

In the book [3, Section 2.3], the classification in the so-called stationary case t = 1 has been dealt with in great detail. The proof was heavily using the principle of continuity of the roots with respect to the parameters p i,j , thus reducing the analysis to the case of the so-called simple random walk, i.e. when only p 1,0 , p 0,1 , p -1,0 , p 0,-1 are different from 0. Indeed, for t ∈]0, 1[, exactly the same arguments hold. Nonetheless, we present hereafter a self-contained proof. The following theorem holds.

Theorem 5.1. For all non singular random walks and for all t ∈]0, 1[, the algebraic curve defined by K(x, y, t) = 0 has genus 0 if, and only if, D(x, t) has a double root at x = 0 or x = ∞. That corresponds to the following pictures.

1.

x 1 = x 2 = 0 y 3 = y 4 = ∞. Γ 0 Γ y 1 y 2 x 3 x 4 X-plane Y -plane 2. y 1 = y 2 = 0 x 3 = x 4 = ∞. y 4 Γ X-plane Y -plane Γ x 1 x 2 0 y 3 3. x 3 = x 4 = ∞ y 3 = y 4 = ∞. x 2 Γ y 1 y 2 X-plane Y -plane Γ x 1 4. x 1 = x 2 = 0 y 1 = y 2 = 0. x 4 Γ X-plane Y -plane Γ 0 y 3 y 4 0 x 3
In other words, multiple roots of D(x, t) can occur only in either of the two following situations:

(a) x = 0, in which case d 0 (t) = d 1 (t) = 0; (b) x = ∞, in which case d 3 (t) = d 4 (t) = 0.
Proof. Let us recall that, since d 2 (t) > 0, the degree of D(x, t) with respect to x is always 2. From the outset, we shall solve the case where this degree is exactly 2.

Degree of D(x, t) = 2. This means that D(x, t has a double root at x = ∞, whence d 3 (t) = d 4 (t) = 0 and the genus of the kernel is 0 for all non singular random walks. Then, by using (2.10) and remarking that d 3 (t) 0, we get the system p 1,0 = p 1,1 p 0,-1 = p 0,1 p 1,-1 = p 1,1 p 1,-1 = 0.

(5.1)

It is now easy to check that the set of relations (5.1), without taking into account the singular random walk p 1,0 = p 1,1 , = p 1,-1 = 0, is tantamount to the pictures 2. and 3.

Degree of D(x, t) 3. We have proved in Theorem 3.1 that a double root can only occur at x = 0. Then d 0 (t) = d 1 (t) = 0, whence by using (2.10) p -1,0 = p -1,1 p 0,-1 = p 0,1 p -1,-1 = p -1,1 p -1,-1 = 0, which corresponds to the pictures 1. and 4.

The proof of the theorem is terminated.

Remark 5.2. When D(x, t) has a double root both at x = 0 and x = ∞, we know from Lemma 4.2 that a factorization of the form (4.2) holds, and the walk is singular. On the other hand, it is worth noting that, for all t ∈]0, 1[, the discriminant D(x, t) can never have double root at x = 1. This case can occur for t = 1 and corresponds random walks having zero drift vectors (see [START_REF] Fayolle | Random Walks in the Quarter Plane: Algebraic Methods, Boundary Value Problems[END_REF]Section 2.3]).

Remark 5.3. The approach proposed in this paper applies verbatim to functional equations pertaining to the transient distribution of two-dimensional Markov processes, such as those encountered in queueing systems (see e.g., [START_REF] Fayolle | Random Walks in the Quarter Plane: Algebraic Methods, Boundary Value Problems[END_REF]Section 8.4]).

  Degree of D(x, t) = 4, i.e. d 4 (t) = 0. Letting |x| → ∞, there are two possibilities. (i) d 4 (t) > 0. Then lim x→∞ φ 2 (x, t) = +∞, and φ 2 (x, t) has a positive real root on ]1, +∞[. (ii) d 4 (t) < 0. Let ξ 1 denote the smallest finite real root of a(x)c(x), when it exists. Then

Remark 4. 3 .Fig. 4. 1 :

 31 Fig. 4.1: Singular random walks.

(

  ii) There exists i, |i| = 1, such that for any j, |j| 1, p i,j = 0 (see figure 4.1b,c); (iii) There exists j, |j| = 1, such that for any i, |i| 1, p i,j = 0 (see figure 4.1d,e). Proof. Let us first eliminate three simple situations. Case a ≡ 0. Then p 1,1 = p 0,1 = p -1,1 = 0. The kernel is of degree 1 in y and this corresponds to the walks of type e in figure 4.1 and K(x, y, t) = t(by + c). Case c ≡ 0. Then p 1,-1 = p 0,-1 = p -1,-1 = 0, giving the walks of type d in figure 4.1 and K(x, y, t) = ty(ay + b). By exchanging the variables x and y, and writing the kernel as a polynomial in x with coefficients a, b, c (see (2.5)), the cases a ≡ 0 or c ≡ 0 depict the walks of respective types c and b in figure (4.1).

  4.1,a); p 1,1 = p 0,1 = p -1,-1 = p 0,-1 = 0, figure (4.1,a) rotated by 90°; p 1,1 = p 1,-1 = p -1,-1 = p -1,1 = 0, figure (4.1,a) rotated by 45°; p 1,1 = p 0,1 = p -1,1 = p 0,-1 = 0, sub-case of figure (4.1,e); (4.10)


Kévin Barré
email: kevin.barre@edu.mnhn.fr

Christian Kerbiriou

Ros-Kiri Ing

Yves Bas

Clémentine Azam

Isabelle Le Viol

Kamiel Spoelstra

Bats seek refuge in cluttered environment when exposed to white and red lights at night

Keywords: acoustic localization, artificial light, flight behaviour, chiroptera, microphone array, streetlight

Background: Artificial light at night is recognized as an increasing threat to biodiversity. However, information on the way highly mobile taxa such as bats spatially respond to light is limited. Following the hypothesis of a behavioural adaptation to the perceived risks of predation, we hypothesised that bats should avoid lit areas by shifting their flight route to less exposed conditions.

Methods: Using 3D acoustic localization at four experimentally illuminated sites, we studied how the distance to streetlights emitting white and red light affected the Probability of bats Flying Inside the Forest (PFIF) versus along the forest edge.

Results:

We show that open-, edge-, and narrow-space foraging bats strongly change flight patterns by increasing PFIF when getting closer to white and red streetlights placed in the forest edge. These behavioural changes occurred mainly on the streetlight side where light was directed.

Conclusions:

The results show that bats cope with light exposure by actively seeking refuge in cluttered environment, potentially due to involved predation risks. This is a clear indication that bats make use of landscape structures when reacting to light, and shows the potential of vegetation and streetlight orientation in mitigating effects of light. The study nevertheless calls for preserving darkness as the most efficient way.

Background

Artificial light at night (ALAN) is recognized as a prominent and growing threat to global biodiversity [START_REF] Koen | Emerging threat of the 21st century lightscape to global biodiversity[END_REF] and hence there is an urgent need to expand scientific knowledge on its effects on ecosystems, and a demand for efficient solutions to reduce these [START_REF] Gaston | The biological impacts of artificial light at night : The research challenge[END_REF]. ALAN impacts a wide range of taxa, at different spatiotemporal scales [START_REF] Altermatt | Reduced flight-to-light behaviour of moth populations exposed to long-term urban light pollution[END_REF][START_REF] Gaston | Impacts of Artificial Light at Night on Biological Timings[END_REF]. Effects vary from the individual level to the disruption of ecosystem functioning by altering interactions between species and regulatory processes [START_REF] Bennie | Artificial light at night causes top-down and bottom-up trophic effects on invertebrate populations[END_REF][START_REF] Hölker | Light pollution as a biodiversity threat[END_REF][START_REF] Knop | Artificial light at night as a new threat to pollination[END_REF].

ALAN also affects spatial behaviour by disorienting species and forming barriers in the landscape. For instance, artificial light disorients migrating birds [START_REF] Doren | High-intensity urban light installation dramatically alters nocturnal bird migration[END_REF] and obstruct toads [START_REF] Van Grunsven | Behaviour of migrating toads under artificial lights differs from other phases of their life cycle[END_REF] and highly mobile taxa such as bats [START_REF] Hale | The ecological impact of city lighting scenarios: Exploring gap crossing thresholds for urban bats[END_REF]. However, the impact of ALAN on species movement across nightscapes remains poorly documented, in particular underlying mechanisms such as changes in spatial behaviour and movement (e.g. flight speed, flight route) of bats in response to light [START_REF] Polak | Differential effects of artificial lighting on flight and foraging behaviour of two sympatric bat species in a desert[END_REF] which potentially affects energetic cost and fitness of individuals [START_REF] La | Energetic consequences of light speeds of foraging red and hoary bats (Lasiurus borealis and Lasiurus cinereus; Chiroptera: Vespertilionidae)[END_REF].

This topic is all the more important given that bats are mostly nocturnal and well known to be impacted by ALAN in term of activity [START_REF] Stone | Impacts of artificial lighting on bats: A review of challenges and solutions[END_REF]. Depending on species, ALAN positively or negatively impacts bat activity (e.g. Pipistrellus spp and Nyctalus spp., respectively Myotis spp., Plecotus spp. and Rhinolophus spp.) at the streetlight scale [START_REF] Azam | Is part-night lighting an effective measure to limit the impacts of artificial lighting on bats?[END_REF][START_REF] Stone | Street Lighting Disturbs Commuting Bats[END_REF][START_REF] Russo | Adverse effects of artificial illumination on bat drinking activity[END_REF], while evidence is also accumulating that ALAN negatively impacts all bat species at larger scales [START_REF] Azam | Disentangling the relative effect of light pollution, impervious surfaces and intensive agriculture on bat activity with a national-scale monitoring program[END_REF][START_REF] Laforge | Reducing light pollution improves connectivity for bats in urban landscapes[END_REF][START_REF] Pauwels | Accounting for artificial light impact on bat activity for a biodiversity-friendly urban planning[END_REF].

ALAN also affects bat movement, for example by keeping individuals from crossing of lit gaps in wooded corridors [START_REF] Hale | The ecological impact of city lighting scenarios: Exploring gap crossing thresholds for urban bats[END_REF] or lit bridges along waterways [START_REF] Barré | Artificial light may change flight patterns of bats near bridges along urban waterways[END_REF] in urban environments. It was also shown that different spectra reduce commuting activity along hedgerows and that light shy species switch to the unlit side [START_REF] Zeale | Experimentally manipulating light spectra reveals the importance of dark corridors for commuting bats[END_REF].

Among possibilities to reduce these impacts, light spectrum, intensity, directionality, light spill and the duration of lighting are parameters that potentially can be used to reduce negative effects [START_REF] Gaston | Reducing the ecological consequences of night-time light pollution: Options and developments[END_REF]. For instance, light-shy bats such as Myotis and Plecotus spp. appear to be equally active close to red streetlights and at unlit sites [START_REF] Spoelstra | Response of bats to light with different spectra: light-shy and agile bat presence is affected by white and green, but not red light[END_REF]. However, highly light averse species such as Rhinolophus hipposideros have shown to avoid all spectra tested for, including red light [START_REF] Zeale | Experimentally manipulating light spectra reveals the importance of dark corridors for commuting bats[END_REF]. Specific part night lighting schemes, with lights turned off from midnight to 5 am, were not found to substantially mitigate effects of light at night as they were still on during the activity peak of bats [START_REF] Azam | Is part-night lighting an effective measure to limit the impacts of artificial lighting on bats?[END_REF]. Moreover, information on impact distances are also essential to prevent negative effects of lighting setups and allow biodiversity friendly urban planning.

However, how the response of bats relates to the distance to light sourcesand hence light intensityis still relatively unknown. Thus far, the only study available on how the response of bats varies with distance reveals clear species dependent differences between 10 and 50 m from a light source [START_REF] Azam | Evidence for distance and illuminance thresholds in the effects of artificial lighting on bat activity[END_REF]. For Eptesicus serotinus no difference in effects were shown between 0 and 10 m from the streetlight, however strong negative effects were present between 25 and 50 m from the streetlight [START_REF] Azam | Evidence for distance and illuminance thresholds in the effects of artificial lighting on bat activity[END_REF]. These findings suggest the response of bats to light is intensity dependent.

All these effects of ALAN and possible measures for reduction of impact on bats remain so far mostly studied using activity metrics (i.e. indicators of abundance). However, the level of activity close to a light source does not provide information on how light level affects the behaviour of bats. Indeed, using a single microphone only allows the assessment of average bat activity within an acoustic detection range of approximately five to over 100 m, depending on species specific call amplitude [START_REF] Barataud | Acoustic ecology of european bats: species identification, study of their habitats and foraging behaviour[END_REF].

Therefore, to assess respective effects of light intensity and spectrum on bats, there is a need for the assessment of species-specific changes in flight behaviour (e.g. changes in flight paths). Acoustic localization in three dimensions (3D) is an effective tool for the assessment of flight paths [START_REF] Koblitz | Arrayvolution: using microphone arrays to study bats in the field[END_REF]. In a study using this technique, authors found that bats reduce flight height and increase flight speed in presence of artificial light [START_REF] Polak | Differential effects of artificial lighting on flight and foraging behaviour of two sympatric bat species in a desert[END_REF]. The attraction of insects by light [START_REF] Owens | The impact of artificial light at night on nocturnal insects: A review and synthesis[END_REF] creates foraging opportunities for bats [START_REF] Rydell | Exploitation of Insects around Streetlamps by Bats in Sweden[END_REF] which should cause bats to reduce flight speed [START_REF] Grodzinski | Context-dependent flight speed: Evidence for energetically optimal flight speed in the bat Pipistrellus kuhlii?[END_REF]. However, studies found that light increases flight speed, probably due to an increased fear of predation [START_REF] Polak | Differential effects of artificial lighting on flight and foraging behaviour of two sympatric bat species in a desert[END_REF][START_REF] Barré | Artificial light may change flight patterns of bats near bridges along urban waterways[END_REF]. An alternative solution to reduce predation risk is to avoid open spaces when exposed to light. Hence, we hypothesize that bats that have the opportunity to fly closer to vegetation (i.e. when flying close to the forest edge) seek shelter in the vegetation while getting closer to lights. This hypothesis is supported by the fact that typically highly light-averse species [START_REF] Azam | Is part-night lighting an effective measure to limit the impacts of artificial lighting on bats?[END_REF][START_REF] Spoelstra | Response of bats to light with different spectra: light-shy and agile bat presence is affected by white and green, but not red light[END_REF][START_REF] Jones | Foraging strategy and predation risk as factors influencing emergence time in echolocating bats[END_REF] such as Myotis and Plecotus species, that are slowflying and hence more vulnerable to predation by hawking birds, mostly fly in cluttered habitats (hereafter named narrow-space foragers). We also hypothesize this behavioural response (i.e. flight closer to vegetation) to begin at least 10 metres from the light source as shown by Azam et al. [START_REF] Azam | Evidence for distance and illuminance thresholds in the effects of artificial lighting on bat activity[END_REF]. Indeed, distances from the light source at which behavioural responses (avoidance or attraction) are detected vary according to species, and approximatively lie around 50 m for bat species mostly flying in open space (hereafter named open-space foragers) such as Eptesicus serotinus, 10 m for bat species mostly flying at wooded edges (hereafter named edge-space foragers) such as Pipistrellus species and at up to 25 m for narrow-space foragers such as Myotis and Plecotus species [START_REF] Azam | Evidence for distance and illuminance thresholds in the effects of artificial lighting on bat activity[END_REF]. Such distance thresholds correspond to light intensities lower than one lux for narrow-space foragers and between one and five lux for open-space and edge-space foragers [START_REF] Azam | Evidence for distance and illuminance thresholds in the effects of artificial lighting on bat activity[END_REF].

In this study, we hypothesize the distance dependent behavioural response of open, edge and narrow-space foraging bat species to streetlights emitting different spectra. Specifically, in comparison with unlit sites we expect bats to fly closer to the vegetation when getting closer to the light, and as much for spectral composition close to white. Using 3D acoustic localization near experimental light posts in forest edges, we first investigated the probability of bats flying inside the forest versus open habitat in relation with the distance to the light. We studied whether this relationship varies around light posts (e.g. the back-and front side, and above and under the lights) in order to determine the potential of impact reduction by light orientation (i.e. shielding by the light armature).

Methods

(a) Experimental sites

The study was done at four experimentally illuminated sites in The Netherlands, each with four rows (separated by 204±17 m) with five four-meter-tall lampposts (separated by 25 m and the central one at forest edge) placed perpendicular in forest edge habitat (Fig. 1). Each row was randomly assigned to emit white, green or red light (Fortimo white, ClearField red, and ClearSky green light, Philips, Amsterdam, The Netherlands), with one of the rows kept dark (just poles). In this study, we only used the white, red and dark rows. All lights are switched on at sunset, and off at sunrise since spring 2012. All experimental lights emit light in the full spectrum range at low intensity; green lamps have an increased blue and reduced red light emission, and red lamps have an increased red and reduced blue emission (Fig. 1).

All light colours have negligible UV emission (see Spoelstra et al. [START_REF] Spoelstra | Experimental illumination of natural habitat--an experimental set-up to assess the direct and indirect ecological consequences of artificial light of different spectral composition[END_REF] for spectral compositions). The light beam of each light is directed downwards by Philips Residium FGS224 (1xPL-L36WHFP) armatures to project light in preferential directions. The light intensity at ground level is on average 8.7 ± 3.0 lux, which is comparable to the illumination levels of countryside roads (Fig. 1; see Spoelstra et al. [START_REF] Spoelstra | Experimental illumination of natural habitat--an experimental set-up to assess the direct and indirect ecological consequences of artificial light of different spectral composition[END_REF] for a further description of these experimental sites). For more detail about light intensity in relation with the distance to the lamp and the orientation of the lamp, see Additional file 1, Appendix S1.

(b) Sampling design and 3D acoustic localization

Bats were acoustically tracked in three dimensions during the first three hours after sunset for 12 nights between the 10 th and 22 nd of July 2018. As we could deploy only one microphone array at the time, we unable to simultaneously sample different light treatments, so we sampled white light, red light and dark control during separate nights. In order to limit variation in bat behaviour linked with inter-night environmental variations, we always sampled a different light colour between consecutive nights (Additional file 1: Table S1). All nights were sampled under highly favourable and constant temperatures (average=16.2 °C, Standard Deviation=1.7°C) and wind speed conditions (average=1.3 m/s, SD=0.8 m/s) (Additional file 1: Table S1). In total four dark control, four red and four white lights were sampled. White light, red light, and dark control were uniquely sampled in two of the four experimental sites, and combinations of the two spectra (dark control and red, and dark control and white, respectively) were sampled in the two other experimental sites. Within the row of light posts, we always sampled the light post right in the forest edge (i.e. at the border of the forest and the open area; see Fig. 1).

To reconstruct 3D positions of bats, we used a trajectography system (hereafter named microphone array) designed at the Institut Langevin (Paris, France). The system uses echolocation calls recorded in a frequency range from eight to 160 kHz at four microphones (FG 3329, Knowles Acoustics, Itasca, IL USA; see Additional file 1, Appendix S2 for more details about sound recording and triggering of echolocation calls). Microphones were arranged in a horizontal triangle form (i.e. three in the corners and one in the middle; Fig. 1).

The microphone array was set up similarly near each row: the four microphones were placed in a horizontal plane above the ground surface, with the central microphone in the open space at four meters horizontal distance to the forest edge and to the streetlight axis perpendicular to the forest edge. The microphone array was always placed on the same side of the light (Fig. 1). Bat positions in the detection range of the microphones were continuously assessed using the time of arrival difference (TOAD) of bat echolocation calls between microphones in the array, using the call waveform [START_REF] Koblitz | Arrayvolution: using microphone arrays to study bats in the field[END_REF]. Further details about the recording setup, the conversion of bat call arrival times into 3D positions, assigning positions to trajectories and the precision of these can be found in Ing et al. [START_REF] Ing | Echolocation Calls and Flight Behaviour of the Elusive Pied Butterfly Bat ( Glauconycteris superba ), and New Data on Its Morphology and Ecology[END_REF] and Additional file 1, Appendices S2 and S3.

(c) Calculation of the distance from 3D positions to the light and accounting for imprecision

To calculate the distance between each bat position and the light, we used the following equation (D, equation 1):

(Eq. 1)
where x, y and z represent distances to the microphone array for each of the three-dimension axis of a given position i. The microphone height was corrected for by entering the actual height of the microphones as placed in the field (i.e. 0.82 to 1 meter) in the position calculation software. Since lights were located at the forest edges and were at four meters height, and given the microphone array placement, we subtract four meters to x, y and z axis, in order to compute the real distance to the light (Fig. 1).

Since the imprecision was expected to increase with the distance to the microphone array, we calculated the cumulated imprecision of 3D positions (I, equation 2) as follows:

(Eq. 2)
where dx, yx and zx represent the standard deviation of distances to the microphone array estimated for each of the three-dimension axis of a given position i [START_REF] Ing | Echolocation Calls and Flight Behaviour of the Elusive Pied Butterfly Bat ( Glauconycteris superba ), and New Data on Its Morphology and Ecology[END_REF]. We choose to discard positions with a cumulated imprecision of more than one meter (Fig. S1

(d) Assigning species to 3D positions

Individual bat calls used to reconstruct 3D positions were saved by continuously recording sound files. In a second step, sound files were segmented into five-second intervals which is sufficient to cover the average duration of a bat pass [START_REF] Kerbiriou | Potential of bat pass duration measures for studies of bat activity[END_REF]. Each of 25,195 five-second files were then classified to the closest taxonomic level using the Tadarida software [START_REF] Bas | Tadarida : A Toolbox for Animal Detection on Acoustic Recordings[END_REF]. Because In a third step, we linked the 3D positions in each 5 second file to the species group found by Tadarida. In case calls of different species groups were found within the same 5 second file, we were able to assign the correct species group to separate series of calls by making use of sequential 3D positions and the peak frequency.

(e) Statistical analysis

We assessed whether the probability of bats flying inside the forest (PFIF) differed according to the distance to the light, and whether this relationship differed between spectra (i.e. dark control, red and white lights). The relationship between the PFIF and the distance to the light allowed us to define the Flight-Path Switch Distance (FPSD) as the distance at which bats on average flew as much inside as outside the forest. We performed Generalized Linear Mixed Models (GLMM, R package TMB), using the PFIF as a binomial response variable where zero corresponded to positions located in the open habitat, and where one corresponded to positions located inside the forest (Fig. 1). We used as explanatory variables the distance to the light, the spectrum type, and the interaction between them to assess the effect of spectra on FPSD. To account for a part of the pseudo-replication (i.e. an average of 15.4 ± 10.1 positions per trajectory; Additional file 1: Fig. S2), we included a random effect on the trajectory identifier. We also included the date as random term in models to control for potential inter-site (i.e. one site sampled each night) and inter-night variations of bat behaviour in relation with lights. Note that weather conditions were highly favourable to bats and stable throughout the sampling period, and that habitat composition was similar between sites (see Sampling design and 3D acoustic localization section and Additional file 1: Table S1).

Given that imprecisions of positions were slightly positively correlated with the distance to the light for Eptesicus/Nyctalus group (r = 0.05, t = 2.0, df = 1786, p-value < 0.045; Pearson's correlation test), we adapted the weight of response variables to the associated precision of positions (i.e. inverse of the imprecision squared [START_REF] Penone | Urbanisation effect on Orthoptera: which scale matters?[END_REF]) by adding a precision weight term in GLMMs.

Lights were oriented toward the ground and the armature parallel to the forest edge, which results in a heterogeneous distribution of light in horizontal and vertical planes (Fig. 1). Thus, in order to assess the dependence of light effects on bats to their spatial position around a streetlight, we built one model per species group for (i) all positions around the light, for (ii) positions under the light (i.e. z < 4 m), for (iii) positions above the light (z > 4 m), for (iv) positions at the backside of the light (i.e. x > 4 m) and for (v) positions in front of the light (i.e. x < 4 m) (Fig. 1). All GLMMs exhibited much smaller Akaike Information Criteria (AIC) than null models. We assessed for each model the goodness of fit by computing the percentage of variance explained by models using the r2 function (R package sjstats). We checked residual plots of models using the R package DHARMa. All analyses were performed using a significance threshold of 5% in the R statistical software [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF].

Results

(a) 3D acoustic localization

Bat calls within the 25,195 five-second files recorded allowed the assessment of a total of 28,646 (3D) positions, with an imprecision of less than one meter. Of all positions, 91.3 % were assigned to the Pipistrellus group, 6.2 % to the Eptesicus/Nyctalus group and 2.4 % to the Myotis/Plecotus group (Additional file 1: Table S2). The number of locations was higher at the white-light poles, followed by red and then dark control poles for Eptesicus/Nyctalus and Pipistrellus groups, but higher around red-light poles followed by white and then dark control poles for Myotis/Plecotus group (Additional file 1: Table S2). Overall, the cumulative imprecision for each location varied between 0.10 and 0.39 m on average, and was dependent on species group, but similar between the three light treatments although slightly lower for all groups in red sites (Additional file 1: Table S2). More than 70 % of positions had a cumulative imprecision lower than 0.2 meters (Fig. S1).

(b) Effect of spectrum on the flight behaviour

Overall, the average probability of bats flying inside the forest (PFIF) was significantly higher near red and white light posts compared to dark control poles for Pipistrellus group, while only significantly higher near white light posts for Myotis/Plecotus and Eptesicus/Nyctalus groups (Table 1). Then, all bat groups were found to have a greater PFIF when getting closer to the light. The increase in PFIF when getting closer to the light was stronger for red and white light posts compared to dark control poles for Pipistrellus and Eptesicus/Nyctalus groups, and only stronger for white light posts for the Myotis/Plecotus group. The increase in PFIF when getting closer to the light was even stronger for red compared to white light posts for the Pipistrellus group, and even stronger for white compared to red light posts for Eptesicus/Nyctalus (Table 1; Fig. 2). Irrespective of these differences, white lighting increased the Flight Path Switch Distance (FPSD) for Eptesicus/Nyctalus group (i.e. 6.1 m), red lighting generated a greater FPSD for Pipistrellus group (i.e. 2.0 m), and white light generated a FPSD of 5.5 m for Myotis/Plecotus group (Additional file 1: Table S3; Fig. 2).

At three to five metres from the light, the PFIF for Myotis/Plecotus and Eptesicus/Nyctalus groups even reached 100 % for white lighting treatment, while the usual PFIF at such distance in unlit sites was under 1 % (Fig. 2). Similarly, the PFIF for Pipistrellus group reached more than 85 % and 50 % at one meter from red and white lights respectively, while close to 0 % in unlit conditions whatever the distance (Fig. 2).

(c) Variation of responses according to location around lights

All species groups increasingly flew inside the forest when getting closer to the light. For both spectra, this effect was only present at the front side of the light, except for Eptesicus/Nyctalus group around red lights (Table 1; Fig. 3). For this group, the response was furthermore limited for bats flying above light posts of both spectra. For the Pipistrellus group, this response was unrelated to flight height (Table 1; Fig. 3).

Concerning the distance dependency, Flight Path Switch Distance (FPSD) was greater for

Eptesicus/Nyctalus flying above and at the front side of white compared to red light poles (5.1 m and 7.0 m versus 4.0 m and 0.2 m, respectively), and was greater for Myotis/Plecotus flying at the front side of white light poles (6.7 m versus no prediction possibility due to insufficient PFIF, respectively) (Additional file 1: Table S3; Fig. 3). The FPSD was also higher under and at the front side of red-light poles compared to white light poles for Pipistrellus group (3.2 m and 2.4 m versus no prediction possibility due to insufficient PFIF, respectively) (Additional file 1: Table S3; Fig. 3).

Finally, fixed effects of models overall almost always explained a large part of the variance (0.11-0.90 r squared; Additional file 1: Table S3).

Discussion

We show that artificial light located at forest edges significantly increases the Probability of For open-and edge-space foraging bats that take advantage of around light accumulated insects, the presence of cluttered habitat (i.e. forest in our case) could further facilitate foraging around streetlights by providing shelter against predators. This result is also consistent with the antagonist effects of ALAN at different spatial scales for open-and edgespace foragers. At the streetlight scale these groups can appear light-opportunistic [START_REF] Stone | Impacts of artificial lighting on bats: A review of challenges and solutions[END_REF][START_REF] Zeale | Experimentally manipulating light spectra reveals the importance of dark corridors for commuting bats[END_REF][START_REF] Spoelstra | Response of bats to light with different spectra: light-shy and agile bat presence is affected by white and green, but not red light[END_REF][START_REF] Azam | Evidence for distance and illuminance thresholds in the effects of artificial lighting on bat activity[END_REF][START_REF] Azam | Is part-night lighting an effective measure to limit the impacts of artificial lighting on bats?[END_REF], however, at a larger scale these species are negatively impacted by ALAN [START_REF] Azam | Disentangling the relative effect of light pollution, impervious surfaces and intensive agriculture on bat activity with a national-scale monitoring program[END_REF][START_REF] Laforge | Reducing light pollution improves connectivity for bats in urban landscapes[END_REF][START_REF] Pauwels | Accounting for artificial light impact on bat activity for a biodiversity-friendly urban planning[END_REF]. The observation that light-opportunistic open-and edge-space foragers seek refuge in cluttered environment near light sources may explain this negative impact, especially in areas with little vegetation around light.

Comparable behavioural changes were expected for Myotis/Plectus species as they are adapted to fly in cluttered environments and are known to be light-averse [START_REF] Azam | Is part-night lighting an effective measure to limit the impacts of artificial lighting on bats?[END_REF][START_REF] Spoelstra | Response of bats to light with different spectra: light-shy and agile bat presence is affected by white and green, but not red light[END_REF][START_REF] Jones | Foraging strategy and predation risk as factors influencing emergence time in echolocating bats[END_REF]. Both red and white light increased the PFIF for Myotis/Plecotus group compared to dark control sites.

Although the effect of red lights was much less important compared to white lights, this finding is important as red light has been reported to have limited to absent effects on the activity of these species [START_REF] Zeale | Experimentally manipulating light spectra reveals the importance of dark corridors for commuting bats[END_REF][START_REF] Spoelstra | Response of bats to light with different spectra: light-shy and agile bat presence is affected by white and green, but not red light[END_REF]. Our results thus suggest that red light may actually not be entirely effective in avoiding behavioural changes of narrow-space foragers, and even less for open-and edge-space foragers.

Overall, Flight Path Switch Distances (FPSD) were mostly longer around white lights, which is likely due to the fact that bats may perceive white light as more intense compared to red light due their spectral sensitivity [START_REF] Gorresen | Ultraviolet Vision May be Widespread in Bats[END_REF][START_REF] Winter | Ultraviolet vision in a bat[END_REF]. Such differences in FPSD could also be linked with differences in light intensity at equal distance, higher for white than red lamps we studied (4.83 more lux in average in a 5 m radius around red lights; see Additional file 1: Appendix S1 for graphical representation of light intensity in relation to the distance to white and red lights), which is known to be one of light parameters driving impacts on bats [START_REF] Lacoeuilhe | The influence of low intensities of light pollution on bat communities in a semi-natural context[END_REF]. When considering the vertical location of bat positions, we found that compared to white light poles, the FPSD was higher for individuals located under (i.e. for Pipistrellus group) and above (i.e.

for Pipistrellus and Eptesicus/Nyctalus groups) red light poles. This is likely directly related to the distribution of the light around the streetlights (supplementary material S1), and aligns with the distance relationand hence intensity dependenceof activity reported by [START_REF] Azam | Evidence for distance and illuminance thresholds in the effects of artificial lighting on bat activity[END_REF].

However, it should be noted that Eptesicus/Nyctalus group mainly flew above lights (Additional file 1: Table S2) which likely explain the absence of response under lights.

We also found changes in bat behaviour in front of light posts for all groups but not at the backside. Thus, the directionality of the light post matters and can be used to reduce the adverse impacts of artificial lighting on bats. In our study, we had a sharper cut-off in light at the backside of the light posts, and hence the effects there disappeared at shorter distances.

Individuals could forage at the backside of streetlights where the predation risk is low, and hence not seek refuge inside forest when getting closer to the light. Further investigations are needed to understand mechanisms involved. Concerning the overall higher effect of red light compared to white light on Pipistrellus group, further studies would be needed to understand why, and if they possibly turn back when getting closer to light instead of seeking refuge inside the forest. However, open-space foragers do not show the same pattern and react similarly to red and white light. We could hypothesise that their higher flight height allows for flying in or above the canopy (i.e. in a potentially more open space than for Pipistrellus group) while increasing their flight speed in response to light, which could explain their different response than edge-space foragers. Further investigations are also needed to address these aspects.

Finally, depending on bat location around streetlights, Flight Path Switch Distances (FPSD) in front of streetlights overall started from 7 m and 4 m for white and red lights, respectively. These distances of impact correspond to light intensities around 6 lux for both white and red lights. However, we defined FPSD as the distance at which bats on average flew as much inside as outside the forest, but impacts likely start before this arbitrarily chosen threshold.

When we look at the beginning of behavioural perturbation, i.e. when the PFIF previously close to zero increases towards positive PFIFs, corresponding FPSD would be around 15 m and 9 m for white and red lights, respectively. Such distances correspond to light intensities around one lux for both white and red lights. These thresholds seem to be consistent with a study which looked at thresholds in light intensity affecting bat activity [START_REF] Azam | Evidence for distance and illuminance thresholds in the effects of artificial lighting on bat activity[END_REF]. However, it is important to be cautious about the definition of safety thresholds for bats and further studies should confirm these results by testing wider distance and intensity ranges around streetlights, by sampling all spectra simultaneously, and by studying more sites and nights per spectrum. We were not able to measure forest height, thus further studies could accurately account for bat position in relation to forest canopy as open space foragers such as Nyctalus spp. can fly above (average flight height 9±4 m; Additional file 1: Fig. S3). However, considerably lower flight heights recorded for Pipistrellus (5.8±2.8 m) and Myotis/Plecotus groups (5.2±2.5 m), and light effects generalized to under and above light positions for Pipistrellus group support these findings (Additional file 1: Fig. S3).

Conclusion

Our study demonstrates that spectrum type, intensity and directionality of streetlights has an effect on the flight behaviour of all bats, including light-opportunistic species, highlighting the need to consider simultaneously all these characteristics when studying ALAN impact on bats. In contrast to the absence of changes in bat activity in response to red light reported earlier, we here show that bats can have a comparable change in flight behaviour in response to red and white light. This finding first shows that bats actively seek refuge in cluttered environment when getting closer to light sources. This is a clear indication that bats make use of landscape structures when dealing with light, and shows the potential of vegetation in mitigating negative impacts of artificial light at night, but calls for preserving darkness as the most efficient way.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Table 1. Estimates, standard errors and p-values of the effect of the distance to the light, the spectrum and the mutual interaction on the probability of bats flying inside the forest when unlit control (A) and white spectrum (B) were used as intercept (***p < .001, **p < .01, *p < .05). Results are presented for all positions, and positions above, under, behind and in front of the light (see Fig. 1 for placement definitions), and derived from generalized linear mixed models.

Eptesicus/Nyctalus

Myotis/Plecotus Pipistrellus

All positions

N=1,788 N=692 N=26,166

  ) and those not included in any bat individual trajectory (i.e. composed of several positions; for details on trajectory reconstruction see Additional file 1: Appendix S3), which led us to keep 28,646 positions on the 35,067 recorded.

  the identification by echolocation to the species level is difficult, we limited identification to following species groups: the Eptesicus/Nyctalus group including Eptesicus sp. and Nyctalus sp., the Myotis/Plecotus group including Myotis sp. and Plecotus sp., and the Pipistrellus group including Pipistrellus sp. These three groups respond differently to light: Eptesicus/Nyctalus group (i.e. open space forager) are usually considered as light opportunistic or light shy (context dependent), species in the Myotis/Plecotus group (i.e. narrow space forager) are light shy, and species in the Pipistrellus group (i.e. edge space forager) are light opportunistic.

Flying

  Inside the Forest (PFIF) for open-, edge-and narrow-space foragers (Eptesicus/Nyctalus, Pipistrellus and Myotis/Plecotus groups, respectively) regardless of light spectrum.

Figure 1 .

 1 Figure 1. Location of study sites plotted on a nightly light emission map (A), schematic

Figure 2 .

 2 Figure 2. Predicted curves of the probability of bats flying inside the forest (left y axis) in

Figure 3 .

 3 Figure 3. Predicted curves of the probability of bats flying inside the forest (left y axis) in

Acknowledgments

We thank Staatsbosbeheer, Natuurmonumenten, the Dutch Ministry of Defence, Het Drentse Landschap and the municipality of Ede for the use of their terrain. We are grateful towards Julie Pauwels and Charlotte Roemer which recently help to improve R script allowing to reconstruct 3D bat positions. We also thank Ray van Dijk for his help with the nightly field work.

Funding

This work was supported by the Infrastructures de Transports Terrestres, Ecosystèmes et Paysages (ITTECOP) program, including funding sources from the Agence de l'Environnement et de la Maîtrise de l'Énergie (ADEME) and the Fondation pour la Recherche sur la Biodiversité (FRB). The setup and maintenance of the experimental facilities of "Licht Op Natuur" is financed by the Dutch Technology Foundation STW, part of the Netherlands Organization for Scientific Research (NWO).

Availability of data and materials

R scripts and data used for analyses are available in the Zenodo repository, https://doi.org/10.5281/zenodo.4036279.

Competing interests

We declare we have no competing interest.

Author's contribution K.B., I.L., C.K, and K.S. conceived the ideas and designed the methodology; R.K.I. built the microphone array; K.B. and K.S. collected the data, K.S. set up experimentally illuminated sites; K.B. and Y.B. analysed the data and wrote R scripts; K.B, I.L, C.K. and K.S interpreted results; all authors led the writing of the manuscript. All authors contributed to the drafts and gave their final approval for publication.