L Leblanc
email: leblanc@isir.upmc.fr

E Saghbini

J Da Silva

A Harlé

S Vafadar

T Chandanson

R Vialle

G Morel

B Tamadazte

Automatic Spinal Canal Breach Detection During Pedicle Screw Placement

Precise pedicle screw placement is imperative for a range of spinal procedures, demanding exact geometric alignment while also carrying inherent risks. The literature reports a complication rate in the case of screw mispositioning of up to 18%. To increase the accuracy of pedicle screw placement, we developed a robotic setup and a breach detection algorithm that could detect a possible perforation of the spinal canal. The robotic setup includes a robotic arm, a drilling system, and specific sensors, e.g., electrical conductivity at the drill bit's tip. The breach detection algorithm consists of a Bayesian-based method providing online and real-time analysis of the electrical conductivity signal to predict a breach. The robotic setup and the perforation detection algorithm were assessed in two ex-vivo experiments. First, data collection was performed by drilling 80 fresh pig vertebrae pedicles, followed by precise data labelling by a surgeon. The evaluation of the proposed algorithm was conducted numerically. Finally, the assessment was performed online by automatically drilling into pedicles in conditions similar to the operating room. The results demonstrated that the algorithm could predict perforations and prevent the robotic setup from going through in 100% of 24 drilled vertebrae. Results demonstrate that using electrical conductivity combined with a robotic setup allows the detection of imminent perforations of the spinal canal during pedicle drilling. This was the first study to evaluate spinal canal perforation detection during ex-vivo pig pedicle drilling.

I. INTRODUCTION

Pedicle screw (PS) placement in spine surgery has been considered a standard surgical procedure for the past forty years. It is widely used for various spinal procedures, especially in the case of spinal deformities, e.g., scoliosis [START_REF] Hwang | Posterior correction of adolescent idiopathic scoliosis with high-Fig. 10. Sample of 5 (among a total of 24) vertebrae drilled experimentally using the developed robotic set-up and the proposed perforation detection algorithm. a) shows a global view of 5 vertebrae, b) the corresponding zoom on the drilled holes, and c) depicts the CT-scan images of the drilled holes. As can be seen (green circle) the drilling was stopped before the spinal canal, i.e. at the moment of deformation of periostal tissue. density pedicle screw-only constructs: 5 years of follow-up[END_REF]. Different forms of scoliosis can occur in both young and adult patients. Clinical grading of scoliosis is based on the severity of the angular deformity, measured by the Cobb angle [START_REF] Langensiepen | Measuring procedures to determine the cobb angle in idiopathic scoliosis: a systematic review[END_REF], defined as the angle between the endplates of the most tilted vertebrae of the curve [START_REF] Cobb | Outline for the Study of Scoliosis[END_REF]. Scoliosis can be idiopathic, creating pathologies like neuromuscular disorders, malformations, and other syndromes [START_REF] Aebi | The adult scoliosis[END_REF]. It can cause muscular pain associated with a continuous effort to maintain the spine's balance or even neurological pain associated with narrowing the spinal canal. From 70 • not 75 • angulation onwards, there are repercussions on respiratory function [START_REF] Koumbourlis | Scoliosis and the respiratory system[END_REF]. Scoliosis treatment depends on the angular deformity's aetiology and severity. In young patients, the risk of evolution and angulation determines the treatment. For instance, when the deformation is below 20 • , the patient will typically require only regular monitoring. For curves between 20 • and 50 • , the patient is treated with a brace to limit the progression of the disease. However, surgical intervention is usually indicated for more severe deformities, i.e., those greater than 50 • . The surgical treatment, also known as arthrodesis, consists of spinal fusion by implantation of PS and rods [START_REF] Cuartas | Use of all-pedicle-screw constructs in the treatment of adolescent idiopathic scoliosis[END_REF]. The PS should be positioned correctly to obtain maximum stability of the bone-screw interface while avoiding a perforation in the spinal canal during insertion (Fig. 1). This highlights that the success of spinal surgery depends substantially on the ability to position the PS safely and accurately [START_REF] Weinstein | Anatomic and technical considerations of pedicle screw fixation[END_REF].

Drilling the pedicles for screw positioning is a challenging procedure. Failure can lead to injuries to major vessels (aorta, vena cava) or the spinal cord [START_REF] Faraj | Early complications of spinal pedicle screw[END_REF]. The literature reports a complication rate in the case of screw mispositioning of up to 18% [START_REF] Tang | Position and complications of pedicle screw insertion with or without image-navigation techniques in the thoracolumbar spine: a meta-analysis of comparative studies[END_REF]. The screw insertion task is performed free-hand, with intraoperative X-ray control or navigation guidance. Patient-specific guides are used to increase accuracy [START_REF] Farshad | Accuracy of patient-specific templateguided vs. free-hand fluoroscopically controlled pedicle screw placement in the thoracic and lumbar spine: a randomized cadaveric study[END_REF]. In some cases, a passive robot is used to immobilize the screw trajectory at the entry point while allowing the surgeon to perform the drilling task [START_REF] Ortmaier | A hands-on-robot for accurate placement of pedicle screws[END_REF]. However, the most established techniques are navigation systems combined with a preoperative registration phase and robotic solutions. A wellknown technical limitation of these solutions is the loss of accuracy due to imaging [START_REF] Joseph | Current applications of robotics in spine surgery: a systematic review of the literature[END_REF] and registration errors [START_REF] Park | Radiographic and clinical outcomes of robot-assisted posterior pedicle screw fixation: two-year results from a randomized controlled trial[END_REF], [START_REF] Pijpker | Accuracy of patientspecific 3d-printed drill guides for pedicle and lateral mass screw insertion: an analysis of 76 cervical and thoracic screw trajectories[END_REF] in addition to the risk of radiation exposure resulting from X-ray imaging.

II. MOTIVATION AND CONTRIBUTIONS

Although navigation approaches have contributed to improving spine surgical procedures, they present several limitations, particularly in terms of usability and accuracy. Accuracy is reduced because of the patient's breathing movement. Furthermore, a slight change in the patient's position during the pedicle preparation is sufficient to introduce registration errors [START_REF] Joseph | Current applications of robotics in spine surgery: a systematic review of the literature[END_REF], [START_REF] Shuman | Intraoperative navigation in spine surgery: Effects on complications and reoperations[END_REF].

One of the advantages of a robotic approach in PS placement, beyond the functional accuracy it provides, is safety. Indeed, the use of a robotic system during the drilling phase allows instantaneous reaction (e.g., stopping the drilling) in tens of milliseconds, compared with manual drilling performed by a surgeon who may be slower to react. This is particularly true in certain procedures (e.g., on adult patients), where the surgeon has to apply relatively great force to drill into the bone. Additionally, drilling parameters (rotation speed, applied force and torque, etc.) can be precisely defined and kept constant (or adapted to suit the patient) throughout the procedure contrary to a manual procedure.

In the remainder of the paper, we will demonstrate that the use of the robot combined with an automatic breach detection method enables drilling to be stopped very early, i.e., within the first milliseconds of a change in the electrical conductivity signal used in our algorithm. In this spirit, we have proposed a robotic platform for PS insertion equipped with specific sensors, including the electrical conductivity sensor (hereinafter referred to as DSG sensor) developed by SpineGuard SA [START_REF] Suess | Control of pedicle screw placement with an electrical conductivity measurement device: initial evaluation in the thoracic and lumbar spine[END_REF], [START_REF] Yurube | Improved accuracy and safety of pedicle screw placement by using a probe with an electrical conductivitymeasuring device during severe syndromic and neuromuscular scoliosis spine surgery[END_REF] and many other sensors (force, position, vibrations, etc.). We then investigate the use of signals recorded by the mentioned sensors to detect perforations early during pedicle drilling. Real-time perforation detection is performed using a Bayesian-based method to detect changes in the behavior of the DSG signal. The developed robotic setup and the perforation detection algorithm are extensively evaluated, numerically and experimentally, based on the Gertzbein-Robbins classification of PS placement errors [START_REF] Gertzbein | Accuracy of Pedicular Screw Placement In Vivo[END_REF]. According to this classification, grade "A" indicates no spinal canal perforation, grade "B" indicates less than 2 mm of perforation in the spinal canal, grade "C" indicates 2-4 mm spinal canal perforation and grade "D" indicates more than 4 mm of perforation. Higher grades are associated with a higher risk of neurological complications. Grades "A" and "B" are considered clinically satisfactory. This paper extrapolates the Gertzbein and Robbins classification scale to pedicle drilling since no screw insertion is performed.

The numerical validation is achieved using a database created by performing 80 drillings. The results obtained are promising. Regarding the detection rate, an imminent perforation point, which corresponds to a deformation of the periosteal tissue that occurs just before the actual perforation, was detected in 25% of the drilling signals. Perforation was detected less than 2 mm after the actual perforation point in the 75% remaining signals. In total, 100% of perforation alerts were raised in a zone corresponding to the grades "A" and "B" for pedicle drilling. Regarding accuracy, the perforation point was detected with an average error of 0.50 mm (standard deviation = 0.71 mm) over the 80 drilled lumbar vertebrae.

An experimental validation run completes these first results. Twenty-four [START_REF] Sun | State recognition of bone drilling with audio signal in Robotic Orthopedics Surgery System[END_REF] vertebrae are drilled. For all vertebrae, the drilling process was automatically stopped just before perforation, exhibiting a 100% success rate. The experimental results reinforce those obtained in simulation and fully meet the requirements of such a surgical procedure.

III. RELATED WORKS AND BACKGROUND

Bone drilling is a routine procedure in various surgical disciplines (orthopaedics, neurosurgery, traumatology). As discussed, the drilling task is a delicate and challenging procedure that requires accuracy and safety. Detecting a potential breach when drilling a bone has consequently become an active research field. Various sensors are used to measure the drilling progress (force, speed, position, etc.). Therefore, identifying a possible breach is equivalent to detecting an unexpected change in the measured signal. Promising methods have been proposed in the literature. For instance, Brett et al. [START_REF] Brett | An Automatic Technique for Micro-Drilling a Stapedotomy in the Flexible Stapes Footplate[END_REF] proposed an algorithm for stapedotomy based on detecting the simultaneous increase of torque and the decrease of the feed force during the drilling procedure. Another similar work suggested tracking the gradient of the force signal (after filtering) to detect unexpected variations in the signal, meaning that the drill is breaching a femur bone [START_REF] Allotta | Study on a mechatronic tool for drilling in the osteosynthesis of long bones: Tool/bone interaction, modeling and experiments[END_REF]. Lee et al. [START_REF] Lee | Control and breakthrough detection of a three-axis robotic bone drilling system[END_REF] proposed using the trend and a threshold on both the torque and feed rate on various bone types (e.g., scapula, skulls, and femurs). Although these methods are promising, qualitative and quantitative evaluations are not provided. More recently, Hu et al. [START_REF] Hu | State Recognition of Pedicle Drilling With Force Sensing in a Robotic Spinal Surgical System[END_REF] suggested using feature functions based on the measured force for different drill shapes. For this purpose, cattle vertebrae drilling signals are used to identify threshold values which are then used to stop the drilling when the force reaches the transitional zone from the cancellous tissue into the cortical tissue. This study performed six experimental drillings; however, only one drilling hole was evaluated, demonstrating that the drilling stopped automatically at 0.97 mm before the breach. Besides, monitoring the electrical current through the drill motor [START_REF] Hsu | A modular mechatronic system for automatic bone drilling[END_REF] was used to detect unexpected changes in the signal to identify a potential breach within the safety zone. However, no metric evaluation is reported on sample size, success rate, or detection accuracy. The use of audio signal has also been investigated for breach detection purposes [START_REF] Sun | State recognition of bone drilling with audio signal in Robotic Orthopedics Surgery System[END_REF] demonstrating that processing the audio signal allowed discriminating cortical and cancellous bone tissue by analyzing the power spectral density. It is reported that using audio signals provides an 88% success rate in finding potential breaches over 35 drillings performed in pig vertebrae. However, no information on accuracy is provided in this study.

Preliminary work on breach detection [START_REF] Silva | Automatic bone breach detection for spine surgery based on bio-electrical conductivity sensing: Ex-vivo experimental validation[END_REF] shows that the DSG signal can provide very interesting results. The reported breach detection algorithm was tested on 104 lamb vertebrae, giving a 100% success rate. However, the drilling was carried out in the spinous processes emerging perpendicularly into the canal, thereby avoiding the tangential breaches that can occur when the pedicles are drilled. Furthermore, the proposed breach detection algorithm includes signal filtering and calibration phases. This method uses several parameters which must be defined beforehand such as absolute and relative conductivity thresholds and the definition of the expected local conductivity gradient (thresholds). Furthermore, it also includes an online filtering phase (moving averaging) of the DSG signal, which may introduce a delay in breach detection.

The following section describes the proposed method, which consists of Bayesian distribution for online inference. This method does not require a filtering phase or calibration parameters.

IV. BREACH DETECTION METHODOLOGY

A suitable and effective perforation detection algorithm should identify meaningful changes in the trend of signals of interest, operate in real-time, be accurate, robust to measurement noise (regardless of the sensor type), and have minimal parameter tuning. The latter indicates that the surgeon should not need to perform complex parameterization before or during the surgical procedure, e.g., depending on the vertebral bone quality.

Probabilistic and recursive methods are suitable candidates for the accurate detection of an abrupt change in a time series. Such change occurs at so-called changepoints (CPs), i.e., time instants when the probability distribution of a time series changes. The Bayesian Online Change Point Detection (BOCPD) algorithm [START_REF] Adams | Bayesian Online Changepoint Detection[END_REF] is one of the most promising CP detection methods and has been demonstrated in different fields (e.g., finance and environment). It has been found to be one of the best-performing methods for offline and online CP detection [START_REF] Burg | An Evaluation of Change Point Detection Algorithms[END_REF], [START_REF] Aminikhanghahi | A survey of methods for time series change point detection[END_REF].

A. Bayesian Online Changepoint Detection

Let x = (x 1 , x 2 , ...., x n ) ⊤ denote the data samples observed over time. In the present work, x is univariate, but the method can be extended to multivariate cases. Let x 1:t denote the data samples observed between the initial time instant 1 and time instant t. CPs occur at unknown time instants between the initial time instant t = 1 and final time instant t = n. These CPs divide x into production partitions [START_REF] Adams | Bayesian Online Changepoint Detection[END_REF]. For each partition ρ, it is assumed that the data samples within that partition are independent and identically distributed and sampled from some probability distribution P (x t |η ρ ), where the parameters η ρ are also considered independent and identically distributed.

BOCPD relies on computing the probability distribution of an intermediate variable, the "run length", i.e., the time since the most recent CP occurred [START_REF] Adams | Bayesian Online Changepoint Detection[END_REF]. The run length at time instant t is denoted as r t . Observing r t = 0 indicates that a CP occurred at time instant t. Arbitrarily, the first data sample observed is considered an additional CP occurring at initial time instant t = 1, implying that r 1 = 0. Let x r t denote the subset of observations associated with the run r t . Between two consecutive time instants, the run length can either be reset to 0 if a CP occurred or be incremented by 1 if no CP occurred:

r t = 0 if a CP occurred r t-1 + 1 otherwise (1)
An illustration is provided in Fig. 2. A hypothetical univariate signal x is plotted in Fig. 2.a. As long as no CP occurs, r t is incremented as shown in Fig. 2.b where solid lines indicate that the run length grows at the next time instant. Significant changes in the distribution of x occur at time instant 4 and time instant 9. At these time instants, r t resets to 0, as illustrated in Fig. 2.b, where the dotted lines indicate that the current run is truncated and the run length drops to 0.

Once a new observation x t is observed, we determine whether the run length r t is incremented or reset to 0. To this end, we compute the posterior probability distribution P (r t |x 1:t ). In other words, BOCPD aims at finding the run length r t at every time instant that best matches the data x 1:t observed so far. In Bayesian terms, for every time instant t, one wants to find r t = arg max r P (r t = r|x 1:t ) where:

P (r t |x 1:t ) = P (r t , x 1:t ) P (x 1:t ) (2)
Since in (2) only the joint distribution over the run length and observed data P (r t , x 1:t ) depend on r t , one can focus on this distribution and rewrite it recursively [START_REF] Adams | Bayesian Online Changepoint Detection[END_REF]:

P (r t , x 1:t ) =
rt-1 P (r t |r t-1 ) P (x t |r t-1 , x r t ) P (r t-1 , x 1:t-1 )

(3) Equation (3) exhibits that P (r t , x 1:t ) can be expressed as a function of P (r t-1 , x 1:t-1 ). This provides a recursive message-passing algorithm for the joint distribution over the current run length and the data, assuming two other factors are computed: the CP prior and the posterior predictive. Fig. 2.c illustrates this message-passing algorithm. In this plot, the circles represent run-length hypotheses. The lines between the circles show the recursive transfer of probability mass between consecutive time instants. A probability distribution is associated with each circle, i.e., run-length hypothesis. For each hypothesis, the probability distribution parameters are obtained by updating the initial parameters prior using every data sample observed since the beginning of the run.

B. Changepoint Prior

P (r t |r t-1 ) is the prior over r t given r t-1 . This prior is non-zero at only two outcomes, since r t only changes to 2 possible values, as shown in [START_REF] Hwang | Posterior correction of adolescent idiopathic scoliosis with high-Fig. 10. Sample of 5 (among a total of 24) vertebrae drilled experimentally using the developed robotic set-up and the proposed perforation detection algorithm. a) shows a global view of 5 vertebrae, b) the corresponding zoom on the drilled holes, and c) depicts the CT-scan images of the drilled holes. As can be seen (green circle) the drilling was stopped before the spinal canal, i.e. at the moment of deformation of periostal tissue. density pedicle screw-only constructs: 5 years of follow-up[END_REF]. This provides computational efficiency to the method and the prior simplifies to [START_REF] Aebi | The adult scoliosis[END_REF], where H is referred to as a hazard function. b) The most probable run length rt drops to 0 when a changepoint is met. c) The trellis on which the passing algorithm lives: at each time instant, the probability of every run length is computed based on the data observed so far during the said run.

P (r t |r t-1 ) =    H(r t-1 + 1) if r t = 0 1 -H(r t-1 + 1) if r t = r t-1 + 1 0 otherwise (4)
C. Posterior Predictive P (x t |r t-1 , x r t ) is the posterior predictive, i.e., the probability that the most recent datum belongs to a given run. Computing this posterior predictive can be challenging. Using conjugate exponential models can allow for more computational efficiency since they provide a closed-form expression of the posterior predictive. For instance, exponential families, which include many of the most common distributions, are good candidates and allow sequentially updating the distribution parameters as new data are observed [START_REF] Adams | Bayesian Online Changepoint Detection[END_REF].

D. Changepoint Detection

The original literature does not explicitly specify a determination criterion to declare a CP at a specific time instant [START_REF] Adams | Bayesian Online Changepoint Detection[END_REF]. Identifying a CP after a single sample of a new distribution may be challenging in certain cases, making it necessary to "wait" for N samples and evaluate the probability of a change happening N samples prior. Therefore, in terms of run length, detecting that r t drops to 0 after a new observation is not an efficient criterion for accurate CP detection. Instead, we detect that r t drops in a [0, N ] range.

E. Bayesian Online Perforation Detector

Building on the BOCPD approach, a Bayesian Online Perforation Detector is proposed. The overall functioning of the said detector is detailed in Algorithm 1. Regarding practical implementation, the following hypotheses are made:

1) A CP is declared if the run length drops in a [0, [START_REF] Brett | An Automatic Technique for Micro-Drilling a Stapedotomy in the Flexible Stapes Footplate[END_REF] range. 2) P (r t |r t-1 ) is assumed to be constant both for simplicity and since no other prior information is available. Therefore, in Equation 4, H(r t ) = λ -1 , where λ is a parameter to adjust the algorithm sensitivity. A value of λ = 250 is empirically established. 3) As a generic assumption, normal distribution unknown mean and variance is assumed as a likelihood. The corresponding conjugate prior is a Normal-Inverse-Gamma distribution. The resulting posterior distribution is a generalized Student's T distribution with center µ t , precision Λ = αtκt βt(κt+1) , and degree of freedom 2α t [START_REF] Murphy | Conjugate Bayesian analysis of the Gaussian distribution[END_REF]:

P (x t |r t-1 , x r t ) = t 2αt (x t |µ t , β t (κ t + 1) κ t α t ) (5)
Parameters µ t , κ t , α t , and β t are initialized with preperforation distribution priors µ 0 , κ 0 , α 0 , and β 0 at time instant t = 0. Then, at each time instant t, we want to model a new possible run of length 0, corresponding to the possibility of a CP. In this case, posterior predictive is computed using post-perforation distribution priors µ prior , κ prior , α prior , and β prior . These priors may be different from the initial values µ 0 , κ 0 , α 0 , and β 0 to encompass the knowledge that a strong rise in the conductivity signal is expected upon perforation. The values of the distribution parameters priors used in the Bayesian Online Perforation Detector are provided in Table I. At each time instant t, we also model the growth for already existing runs. This is achieved by computing the posterior predictive using parameters µ t , κ t , α t and β t for each run length value before updating these values with the latest datum observed:

TABLE I DISTRIBUTION PARAMETERS PRIORS USED IN

µ t+1 = κ t µ t + x t+1 κ t+1 (6)
κ t+1 = κ t + 1 (7) α t+1 = α t + 0.5 (8)
β t+1 = β t + κ t (x t -µ t ) 2 2(κ t + 1) (9)

V. DATA ACQUISITION SETUP AND PROTOCOL

To assess the proposed perforation detection algorithm, we conducted an ex-vivo experiment for data acquisition (to create a comprehensive database) and evaluated the proposed method. The experiments are done in conditions close to those of an operating room by drilling fresh pedicles and mimicking the physiological behaviour of the cerebrospinal fluid filling the spinal canal. The used drill bit is derived from instruments developed by SpineGuard SA (Fig. 3.a. It consists of a 3 mm diameter drill bit for cutting purposes with a pyramidal tip embedding a conductivity sensor and a threaded shaft (Fig. 3.b). The threaded shape of the drill continuously maintains the contact between the drill and the bone and avoids unwanted motion along the z-axis (i.e., the drilling axis). The drilling is carried out by a custom-made power drilling unit delivering 1.5 Nm nominal torque and a maximum speed of 922 rpm (revolutions per minute), able to perform drilling and screwing tasks. Finally, the drilling system is mounted on the robot's end-effector with an angle of 30 • concerning the last robot axis (Fig. 3.c).

A. Drill bit and drilling system

B. Robotic Setup and Embedded Sensors

The drilling system is rigidly fixed on a robotic arm, with an angle of 30 • , as can be seen in Fig. 4. This robotic arm consists of a KUKA LBR Med 7 R800 robot adapted to medical requirements. It is equipped with built-in joint Calculate posterior predictive:

P (xt|r t-1 , x r t ) 4:
Calculate growth probability:

P (rt = r t-1 +1, x 1:t ) = (1-λ -1
prior )P (xt|r t-1 , x r t )P (r t-1 , x 1:t-1 ) 5:

Calculate changepoint probability:

P (rt = 0, x 1:t ) = r t-1 λ -1
prior P (xt|r t-1 , x r t )P (r t-1 , x 1:t-1 ) 6:

Calculate marginal probability: P (x 1:t ) = r t P (rt, x 1:t ) 7:

Calculate run length distribution: P (rt|x 1:t ) = P (rt, x 1:t )/P (x 1:t ) 8:

Find rt = arg max r P (rt = r|x 1:t )

9:

Check that run length has exceeded the detection threshold once: 10:

if rt > threshold then 11:

DetectionF lag ← T rue 12:

end if 13:

Evaluate whether a change occurred: 14:

if DetectionF lag & rt < threshold then 15:

Alert ← T rue 16:

end if 17:

Update distribution parameters 18: end for torque sensors allowing it to operate in a collaborative mode. It has a payload of 7 kg, featuring position accuracy of ±0.15 mm and joint redundancy (7 degrees of freedom). The DSG sensor of the drilling system allows for distinguishing the different layers the drill bit passes through during the drilling process. The conductivity signal is low in the cortical layers compared to the one detected in the cancellous layers. However, the signal increases significantly once it intersects body fluid. Fig. 5 depicts an example of a typical DSG (proportional to the electrical conductivity) signal acquired by drilling a fresh pig vertebra. These tissue discrimination characteristics are suitable for detecting perforations (lateral or medial). Although only the conductivity signal is used in the perforation detection algorithm, different signals issued from various embedded sensors (torque, position, velocity, time, video, etc.) are recorded and added to the database.

C. Drilling Workflow

We designed a drilling box where the vertebra (Fig. 6.a) is rigidly fixed in the drilling box. To simulate experimental conditions close to actual pedicle drilling in the operating room, the vertebra is submerged in a saline solution that approximately reproduces the electrical conductivity of cerebrospinal liquid contained in the spinal canal (Fig. 6.b). The saline solution is obtained by mixing 7g of salt per litre of tap water at a temperature in a [19.7, 24.6] • C range. A camera is placed in front of the spinal canal to record the whole drilling procedure (Fig. 6.c).

D. Vertebrae Collection and Data Acquisition

To assess the developed methods and material, we collected fresh pig lumbar vertebrae from the butcher (Fig. 6.a). The fresh vertebrae were collected one day before the lab experiments and placed only in a refrigerator (since a freezer could have changed the properties of the bones or soft tissues).

For each drill, the surgeon collaboratively manipulated the robot to place the drill bit at the desired entry point and voluntarily aim for the spinal canal. They then adapted the contact force and drilling rotational speed according to the bone quality of the vertebrae. The force controller was developed earlier in another study [START_REF] Silva | Force control of the KUKA LBR Med without external force sensor[END_REF]. Principally, the required force for drilling was empirically determined and set to f d = 25 N. However, based on the bone's mechanical properties, the force f d could reach up to 100 N. The drilling rotational speed was set at 60 rpm for all the performed drillings.

We recorded several additional signals such as torque, position (depth penetration of the drill bit), velocity, and time. Moreover, the camera filmed the spinal canal and captured the exact instant when the drill perforated into the canal.

A posteriori, three people independently watched the recorded video to identify the perforation time to avoid potential bias. The surgeon then labelled the exact time instant the perforation started. This allowed us to label all the other signals to build the database (see range highlighted in green in Fig. 7). This provides the ground truth for further assessment of the proposed perforation detection algorithms.

In total, 80 lumbar vertebrae were drilled.

VI. NUMERICAL AND EXPERIMENTAL VALIDATIONS

The first evaluation of the proposed perforation detection algorithm was carried out on the collected data. Hence, we evaluated offline the electrical conductivity signals recorded for each drilling (80 in total). Note that for each signal, the surgeon defined the exact time instant of perforation by a posteriori visualizing the videos recorded during the drilling. This task allows grading perforation detection according to our derivation of the Gertzbein-Robbins classification [START_REF] Gertzbein | Accuracy of Pedicular Screw Placement In Vivo[END_REF]. Detecting an imminent perforation before the actual point of perforation corresponds to grade "A," while declaring a perforation less than 2 mm after the actual perforation point corresponds to grade "B." In other words, the detection is considered acceptable if a perforation alert is raised before the right-most limit blue line of the green area in Fig. 7. For the 80 drilled vertebrae, the success rate of the perforation detection is 100% meaning that the perforation was detected before it actually occurred or just at the instant of deformation of the periostal tissue (inner surface covering the spinal cord). From a clinical point of view, this corresponds to grades "A" and "B" according to the Gertzbein-Robbins classification [START_REF] Gertzbein | Accuracy of Pedicular Screw Placement In Vivo[END_REF]. In detail, 25% of the drillings are classified as "A" (imminent perforation is detected before the actual perforation point) and 75% as "B".

A. Numerical Validation

In terms of accuracy, our algorithm allows stopping the drilling with an average error (i.e., distance from the perforation point) of 0.50 mm with a standard deviation (STD) of 0.71 mm.

Figure 7 shows some examples of electrical conductivity signals for which the perforation is detected. As can be seen, the detection of the perforation (red lines shown in Fig. 7) corresponding to the abrupt change in the signal can occur before or slightly after the perforation point. Variability in the electrical conductivity signals can lead to detection after the perforation point, but always less than 2 mm after said point. The perforation detection over the whole set of the 80 performed drilling is summarized in Fig. 8, where Fig. 8.a depicts the distance of the automatic stop to the perforation point for each drilled vertebra, and Fig. 8.b shows the histogram of all the computed distances.

B. Experimental Results

We experimentally assessed the robotic setup and the perforation detection algorithm based on the numerical eval- uation and obtained promising results regarding accuracy and safety during the drilling procedure.

The experimental procedure is identical to the one used to build our database. For this purpose, we used another group of fresh vertebrae. The surgeon moved the drill bit to the entry point of the pedicle (by manipulating the robotic arm), and the drilling was automatically performed. The robot drilled 24 fresh lumbar pig vertebrae. The perforation detection algorithm stopped the drilling online automatically (without any calibration or pre-setting procedure compared to the numerical validation procedure (Section VI-A)). Once the vertebrae were drilled, the surgeon cut the pedicles along the drilled hole to verify where the drilling stopped. After investigation, it was shown that for all 24 drilled vertebrae, the drilling process was automatically stopped just before the perforation, i.e., exhibiting a 100% success rate. Furthermore, the robotic platform and the breach detection method were evaluated in another context, i.e. drilling a fresh sheep skull. It was shown that breaches were avoided (several drillings were carried out) without any change to the algorithm parameters or any further calibration phase being required.

VII. CONCLUSION

Spinal canal perforation can be associated with screw instability and neurological complications during PS placement. To avoid such perforations, we developed a robotic setup equipped with a multitude of sensors as well as a perforation detection algorithm using the electrical conductivity signal. This method was based on a Bayesian approach, allowing for accurate real-time perforation detection without filtering the signal or prior calibration.

To evaluate the proposed materials and methods, we created a database including 80 drillings in fresh pig lumbar vertebrae. This database contains different types of signals (electrical conductivity, torque, velocity, and position). Numerical and experimental validation on perforation detection was convincing (perforations were avoided in 100% of cases). The results show that electrical conductivity allows the detection of imminent perforations of the spinal canal during actual ex-vivo pedicle drilling. Combined with a robotic drilling setup, this permitted reliably stopping the drilling upon detection in real-time.

Future work will evaluate the method's robustness through in-vivo tests. These will involve a pig whose spine is subjected to physiological motions such as breathing and deformations resulting from the pressure of the drill on the pedicle. The other future task will be to extend the drilling database and make it open-source by including several other signals (position, speed, force, time, video, etc.).

Fig. 1 .

 1 Fig. 1. Different scenarios for screw insertion into the vertebrae: a) leading to a medial perforation (i.e. spinal canal perforation), b) leading to a lateral perforation and c) well-placed screw.

Fig. 2 .

 2 Fig. 2. Illustration of the BOCPD method for a) hypothetical univariate signal with distribution changes occurring at time instants t = 4 and t = 9.b) The most probable run length rt drops to 0 when a changepoint is met. c) The trellis on which the passing algorithm lives: at each time instant, the probability of every run length is computed based on the data observed so far during the said run.

Fig. 3 .

 3 Fig. 3. a) Used drill bit equipped with an electrical conductivity sensor designed by SpineGuard SA, b) zoom on the "threaded" part of the drill and c) represent the drill bit mounted on the robotic arm.

Algorithm 1

 1 Bayesian Online Perforation Detector using the electrical conductivity signal Input: xt, DSG (electrical conductivity) signal in mV Output: Alert, flag used to stop the drilling 1: Initialization: P (r 1 = 0) ← 1, DetectionF lag ← F alse, 2: for xt do 3:

Fig. 4 .

 4 Fig. 4. Photography of the developed robotic set-up for spine surgery.

Fig. 5 .

 5 Fig. 5. An example of DSG signal acquired during a vertebra drilling process.

Fig. 6 .

 6 Fig. 6. Illustration of the data acquisition workflow: a) fixing the vertebra in salty water for drilling, and b) view of the frontal camera monitoring the drilling progress.

Fig. 7 .

 7 Fig. 7. Four examples of detected perforations (see vertical red lines) on the recorded DSG signals. A ±2 mm green zone (delimited by the vertical blue lines) centred on the perforation instant (dotted lines) is plotted for visualization. Detection must occur before the upper limit of the green zone to be considered satisfactory (i.e. grade "A" or "B" perforations).

Fig. 8 .

 8 Fig. 8. Qualitative analysis over the 80 performed drillings. a) distance of the automatic stop to the bone interface for each drilled vertebra according to our derivation of the Gertzbein-Robbins classification (the boundaries of the "A" and "B" grades are represented by a red line) and b) the histogram of the resulting computed distances over the 80 drilled vertebrae.

Figure 9 .

 9 a shows an example of an automatic stop of a drilling task before the perforation on the spinal canal with a safety margin of approximately 1 mm (Fig. 9.b).

Fig. 9 .

 9 Fig. 9. Illustration of an automatic stop of a drilling process before perforation.

Figure 10

 10 Figure 10 depicts a sample of 5 (among a total of 24) drilled vertebrae after cutting the pedicles along the drilled holes. Fig. 10.a shows a global view of the 5 vertebrae and Fig. 10.b shows the corresponding zoom on the drilled holes, while Fig. 10.c depicts the CT-scan images realized retrospectively. It can be seen on the CT-scan images that the proposed algorithm, coupled with the robotic platform, makes it possible to stop just at the moment of deformation of the periostal tissue covering the inside of the spinal canal. This means that although the detection error of the perforation is

This work has been supported by EU's H2020 research and innovation program under grant agreement No. 101016985 (FAROS. project)


F Chapoton

On a q-analogue of the Zeta polynomial of posets

research documents, whether they are published or not. The documents may come

In the study of finite partially ordered sets (posets), one uses frequently polynomial invariants. They are useful to distinguish the posets or recognize them under different disguises, but also for testing the solidity of our understanding. Computing these invariants may require a good handle on the combinatorial mechanism behind the scene. This can lead to structural results. For example, sometimes these polynomials factor nicely and one would like to understand the reason for this.

One of the most classical such polynomial is the Zeta polynomial, whose values at positive integers count chains of elements. This polynomial is related quite closely to the order polynomial, which describes integer points in a polytope naturally attached to the poset. Yet another polynomial is the characteristic polynomial, recording the values of the Möbius function. An important but less well-known polynomial is the Coxeter polynomial, which contains information about the derived category of modules over the incidence algebra.

The aim of the present article is to introduce and study a q-analogue of the Zeta polynomial, defined for posets endowed with an height function. This is a polynomial in the variable x, whose values at q-integers count chains according to the sum of heights of their elements. We will show that it shares many properties of the Zeta polynomial, to which it reduces when setting q = 1.

The definition of the q-Zeta polynomial is motivated by the introduction in [START_REF] Chapoton | q-analogues of Ehrhart polynomials[END_REF] of a q-analogue of the Ehrhart polynomial for polytopes. A special case of this q-Ehrhart theory provides a q-analogue of the order polynomial. By following the known relationship between Zeta-polynomials of distributive lattices and order polynomials of their poset of join-irreducibles, this leads to the present construction.

The reader should be warned that, in some sense, a large part of the surrounding context and proofs are not new and indeed very classical since the pioneering works of Stanley on P -partitions. The precise relation of q-Ehrhart polynomials with the theory of P -partitions has been described in [9, §4]. A nice reference on the history of the subject of P -partitions can be found in [START_REF] Gessel | A historical survey of P -partitions[END_REF],

Introduction

Let us first recall briefly the classical Zeta polynomial of posets.

Let P be a finite poset. For every integer n ≥ 2, consider the set C n (P ) of chains e 1 ≤ e 2 ≤ • • • ≤ e n-1 in P . Then the cardinality of C n (P ) is given by the value at n of a polynomial Z P (x), called the Zeta polynomial of P .

One can easily prove this fact by the following computation. For k ≥ 1, let Ch k (P ) be the set of strict chains with k elements in P , i.e. sequences c 1 < • • • < c k .

By gathering chains according to their underlying strict chain and multiplicities, one finds that

Z P (n) = e1≤•••≤en-1 1 = k≥1     m 1 ,...,m k ≥1 i mi=n-1 1     # Ch k . (1)
But the inner sum is just the binomial coefficient n-2 k-1 , so that

Z P (n) = k≥1 n -2 k -1 # Ch k , (2)
which is obviously a polynomial evaluated at n. The Zeta polynomial is an useful invariant of posets, with nice general properties. By the formula [START_REF] Baumann | Chains in shard intersection lattices and parabolic support posets[END_REF], it belongs to the ring of polynomials with rational coefficients taking integer values on N. The Zeta polynomial is multiplicative with respect to the Cartesian product of posets, additive with respect to disjoint union and invariant by duality. For more on this classical subject, see for instance [19, §3.11].

q-analogue

Let us now turn to our proposal for a q-analogue of the Zeta polynomial.

The letter q stands for the indeterminate in Q(q). We use the standard notation [n] q = (q n -1)/(q -1) for the q-analogue of the integer n ∈ Z, which is a polynomial in q if n ≥ 0 and a Laurent polynomial otherwise. We will denote by [n]! q the q-factorial of n when n ≥ 0. Note that [0] q = 0 and [1] q = 1.

Let P be a finite poset. We need the additional data of an height function h : P → N such that h(x) < h(y) for every cover relation x < y in P .

A poset P is graded if there exists an height function that increase by 1 along every cover relation. Then there is a preferred choice for such an height function h, by assuming further that it has minimal value 0 on every connected component of P . This specific height function will be denoted by rk.

Every poset P can be endowed with an height function by choosing an arbitrary linear extension and using it as an height function.

Note that some posets that are not graded are nevertheless naturally endowed with natural height functions, for example the Tamari lattices using their description as posets of tilting modules.

For a strict chain c = (c 1 , . . . , c k ) in Ch k (P ), let h(c) denote the sequence of heights (h(c 1 ), . . . , h(c k )) and let h(c) denote the sum of this sequence. The hypothesis on h ensures that the sequence h(c) is strictly increasing.

For every integer n ≥ 2, let us define

Z P,h ([n] q ) = e1≤•••≤en-1 q j h(ej ) = k≥1 c∈Ch k (P ) m 1 ,...,m k ≥1 i mi=n-1 q i mih(ci) . ( 3
)
This can be rewritten as

Z P,h ([n] q ) = k≥1 c∈Ch k (P ) q h(c) m ′ 1 ,...,m ′ k ≥0 i m ′ i =n-k-1 q i m ′ i h(ci) . (4)
At this point, one can recognize the innermost sum as being essentially the q-Ehrhart polynomial of a simplex. Before proceeding further, let us recall the theory of q-Ehrhart polynomials, whose details can be found in [START_REF] Chapoton | q-analogues of Ehrhart polynomials[END_REF].

q-Ehrhart polynomials and reciprocity

Fix a lattice polytope Q and an integral linear form ℓ on the ambient lattice. Assume that ℓ is not constant on any edge of Q and takes values in N on Q.

Then there exists a unique polynomial

Ehr Q,ℓ (x) ∈ Q(q)[x] such that Ehr Q,ℓ ([n] q ) = z∈nQ q ℓ(z) , (5)
for every integer n ≥ 0. In words, this is counting lattice points in the dilates of Q according to the value of ℓ on each lattice point. By the formula above, the polynomial Ehr Q,ℓ belongs to the ring of polynomials whose value at every q-integer [n] q with n ≥ 0 is a polynomial in q with positive integer coefficients. This is called the q-Ehrhart polynomial of the lattice polytope Q with respect to the linear form ℓ. Note that setting q = 1 recovers the classical Ehrhart polynomial of lattice polytopes, not depending on the linear form ℓ. The degree of Ehr Q,ℓ is the maximum value of ℓ on Q.

Let us also recall Ehrhart reciprocity in this setting. Let n ≥ 1 be an integer. Then the evaluation of Ehr Q,ℓ at [-n] q is given up to sign by the similar sum over interior points in the dilates of Q:

Ehr Q,ℓ ([-n] q ) = (-1) d z∈Int(nQ) q -ℓ(z) , ( 6
)
where d is the dimension of the polytope Q and Int denotes the interior of a polytope. Note that the interior of a 0-dimensional polytope (a point) is just itself.

For k ≥ 1, the lattice polytope in N k whose vertices are the basis vectors in N k will be called the standard basic simplex.

For every k-tuple a of distinct elements of N, let us denote by E a the q-Ehrhart polynomial of the standard basic simplex in N k with respect to the linear form given by the standard scalar product with a. The degree of E a is the maximal element of a. At q = 1, E a becomes the Ehrhart polynomial of the standard basic simplex, namely x+k-1 k-1 , no longer depending on a. For example, for the tuple a = (1, 2, 3), one gets E a = ((q -1) x + 1) • (qx + 1) • (q 2 x + q + 1) q + 1 .

Definition of q-Zeta polynomial

Let us go back to our proposed definition (4) for Z P,h . One therefore finds that

Z P,h ([n] q ) = k≥1 c∈Ch k (P ) q h(c) E h(c) ([n -k -1] q ). (7)
Here one has to be a little cautious about the substitution of the inner summation in (4) by the q-Ehrhart polynomial E h(c) . This is a priori allowed only if n ≥ k + 1. On the one hand, if 2 ≤ n ≤ k, then the inner summation in (4) vanishes because it runs over an empty set. On the other hand, by Ehrhart reciprocity, the value E h(c) ([-d] q ) vanishes when 1 ≤ d ≤ k -1, because in this case the d-th dilate of the standard basic simplex with k vertices has empty interior. By the appropriate shift of variables relating [n] q and [n -k -1] q , one finally reaches the following definition, consistent with all previous formulas. Definition 1. The q-Zeta polynomial Z P,h of a finite poset P with respect to the height function h is given by

Z P,h (x) = k≥1 c∈Ch k (P ) q h(c) E h(c) x -[k + 1] q q k+1 . ( 8
)
This polynomial is an element of the ring Q(q)[x]. Its degree is the maximal value of the height function h on P . It belongs to the sub-ring of polynomials whose values at every q-integer [n] q with n ≥ 2 is a polynomial in q with positive integer coefficients.

The first and second properties are clear from the formula (8) and the general properties of the q-Ehrhart polynomials. The third property follows directly from (4), which holds by construction, as well as [START_REF] Bergeron | Posets related to the connectivity set of Coxeter groups[END_REF].

In other words, by (3), one has the following statement:

Lemma 1.1. The values of Z P,h at q-integers [n] q for n ≥ 2 are q-analogues of the numbers of chains e 1 ≤ • • • ≤ e n-1 in P , where the power of q is the sum of the heights of elements in the chain.

As expected, the polynomial Z P,h deserves the name of q-analogue of the Zeta polynomial.

Lemma 1.2. For any height function h, the specialisation of Z P,h at q = 1 is the usual Zeta polynomial Z P of the poset P .

Proof. This follows by comparing ( 1) and (3).

Examples

Let us now give a few examples, using the height function rk coming from the grading of the posets, unless stated otherwise.

Example 1.3. For the unique poset • with one element with respect to rk, the q-Zeta polynomial is 1. More generally, for the height function on • with value H ∈ N, the value of the q-Zeta polynomial at [n] q is q (n-1)H and therefore

Z •,H (x) = 1 + (q -1)x q H . (9)
Example 1.4. For the total order on d ≥ 2 elements, one finds

d-2 j=0 ([j] q + q j x) [d -1]! q . ( 10
)
This follows from the identification of the value at [n] q , using the formula (3) as a weighted sum over all chains, with the standard q-binomial coefficient counting lattice paths in a (n -1) × (d -1) rectangle according to the area below.

Example 1.5. For the graded poset on 5 elements with one minimum 0, one maximum 1 and 3 pairwise incomparable elements in between, one finds

x • ((q + 2) x -1) q + 1 . ( 11
)
This was found using a computer. By hand, it can be computed by interpolation, as the degree is known to be 2. This is made easier by the statements about the values at [0] q and [1] q given below in lemma 2.3 and lemma 2.2.

Example 1.6. Let us also consider the poset on {a, b, c, d} where a and b are both smaller than both c and d. One obtains

2(q + 1) q (x -1), (12)
whose value at [START_REF] Athanasiadis | Shellability of noncrossing partition lattices[END_REF] q is 0.

Example 1.7. For the poset on {a, b, c} where a is smaller than b and c, one

obtains 2x -1, (13)
whereas for the dual poset one gets

((q + 1) x -1) q . ( 14
) Example 1.8. For the poset on {a, b, c, d, e} where a ≤ b, a ≤ c, b ≤ d, c ≤ e, one obtains 2qx 2 + 2x -q -1 q + 1 . ( 15
)

Basic properties

Let us present in this section a few basic properties of Z P,h . Let P and Q be two posets with height functions g and h. Consider P × Q with the height function g + h sending (a, b) to g(a) + h(b). Consider also P ⊔ Q with the height function g ⊔ h defined by g on P and h on Q.

Lemma 2.1. With the notations above, Z P ×Q,g+h = Z P,g Z Q,h and Z P ⊔Q,g⊔h = Z

P,h + Z Q,h .
Proof. This is most easily seen using formula [START_REF] Bergeron | Posets related to the connectivity set of Coxeter groups[END_REF].

A chain e in the Cartesian product P × Q is the same as a pair of chains (e ′ , e ′′ ) in P and Q. The sum over e of the heights in P × Q is that of e ′ in P plus that of e ′′ in Q.

For the disjoint union P ⊔ Q, chains are either entirely in P or entirely in Q and the result follows.

For a poset P , let χ P denote the Euler characteristic of the order complex of P , which is the simplicial complex made of strict chains. Recall that if P has a unique minimum or maximum, the order complex is contractible.

Lemma 2.2. The value of Z P,h at [START_REF] Athanasiadis | Shellability of noncrossing partition lattices[END_REF] q is χ P . In particular, if P has a unique minimum or maximum, this is 1.

Proof. Let us first compute E h(c) ([-k] q ) for k ≥ 1 and c ∈ Ch k (P ). By Ehrhart reciprocity, E h(c) ([-k] q ) = (-1) k-1 q -h(c)
because the only interior point in the k-th dilate of the standard basic simplex is (1, . . . , 1). Then using formula [START_REF] Brady | K(π, 1)'s for Artin groups of finite type[END_REF], one finds

Z P,h ([1] q ) = k≥1 (-1) k-1 # Ch k (P ).
This is exactly the expected Euler characteristic.

Lemma 2.3. Assuming that P has a unique minimum 0, the value of Z P,h at [0] q is q -h( 0) (1 -χ P \ 0). In particular, if P also has a unique maximum 1 distinct from 0, this is 0.

Proof. The special case when P has only one element is clear, see example 1.3. Assume now that P has at least 2 elements. First one can show that

E h(c) ([-k -1] q ) = (-1) k-1 q -h(c) i q -h(ci) ,
using Ehrhart reciprocity, for k ≥ 1 and c ∈ Ch k (P ). Therefore one finds

Z P,h ([0] q ) = k≥1 (-1) k-1 c∈Ch k (P ) i q -h(ci) .
In this sum, the strict chain reduced to 0 contributes the term q -h( 0) . On the remaining chains, removing or adding 0 define bijections exchanging strict chains with k + 1 elements containing 0 and strict chains with k elements not containing 0. The contributions of such a pair of chains to the previous sum almost cancel mutually, except one summand q -h( 0) in the inner sum, contributing

k≥1 (-1) k c∈Ch k (P \ 0) q -h( 0) , in which one can recognize -q -h( 0) χ P \ 0.
Recall that a poset is bounded if it has unique minimum 0 and unique maximum 1.

Lemma 2.4. If the poset P is bounded, then Z P,h ([-1] q ) = q -h( 0)-h ( 1) µ P ( 0, 1), [START_REF] Stanley | Ordered structures and partitions[END_REF] where µ P is the usual Möbius function of the poset.

Proof. The proof is very similar to that of the previous lemmas.

The special case when P has only one element is clear, see example 1.3. One can therefore assume that 0 = 1.

Using Ehrhart reciprocity to evaluate

E h(c) at [-k -2] q for c ∈ Ch k (P ), one gets the formula k≥1 (-1) k-1 c∈Ch k (P ) i   q -2h(ci) + i =j q -h(ci)-h(cj)   .
Then one first separates strict chains according to whether they start by 0 or not. Contributions of pairs of chains almost cancel. In the remaining sum, one separates strict chains according to whether they end by 1 or not. Once again, there are cancellations by pairs. There remains only

q -h( 0)-h( 1)   -1 + k≥1 (-1) k-1 # Ch k (P ′ )   , where # Ch k (P ′ ) is the number of strict chains c 1 < c 2 < • • • < c k in the poset P ′ = P \ { 0, 1}
. By a classical result, the reduced Euler characteristic between the parentheses is the Möbius number µ P ( 0, 1).

Let us consider a poset P and its dual poset P . To any height function h on P and for any integer H at least equal to the maximal value of h, the function H -h is an height function on P . In that case, the q-Zeta polynomial of P w.r.t. h and that of P w.r.t. H -h are not directly related in a simple way, but their values are related as follows.

Lemma 2.5. One has the following relation:

Z P ,H-h ([n] q ) = q (n-1)H (Z P,h ([n] q )) | q=1/q , ( 17
)
for all n ∈ Z.

Proof. For n ≥ 2, this follows directly from (3). The case n = 1 also follows from lemma 2.2. It is therefore enough to check that the right hand side is the value of a polynomial at the q-integer [n] q . Indeed, let us write

Z P,h (x) = H j=0 κ j x j ,
as the degree of this polynomial is at most H. Then the right hand side in the expected relation is

H j=0 κ j | q=1/q (q n-1 ) H-j [n] j q ,
which is the value at [n] q of the polynomial

H j=0 κ j | q=1/q 1 + (q -1)x q H-j x j .
See example 1.7 for a simple case of this relationship between q-Zeta polynomials of dual posets.

Lemma 2.6. Let P be a poset with height function h. Let Z P,h (x) be the corresponding q-Zeta polynomial. Then the q-Zeta polynomial Z P,h+1 for P with the shifted height function h + 1 is given by

Z P,h+1 (x) = 1 q (1 + (q -1)x) • Z P,h (x). ( 18
)
Proof. Using the formula (3) for the height function h + 1 gives directly

Z P,h+1 ([n] q ) = q n-1 Z P,h ([n] q ),
from which the result follows.

Lemma 2.7. Let P be a poset with height function h. Let Z P,h (x) be the corresponding q-Zeta polynomial. Let D ≥ 1 be an integer. Then the q-Zeta polynomial Z P,Dh for P with the scaled height function Dh is given by

Z P,Dh (x) = Z P,h | q=q D (1 + (q -1)x) D -1 q D -1 . ( 19
)
Proof. From (3), it follows that Z P,Dh

([n] q ) = Z P,h | q=q D ([n] q D ) for n ≥ 2.
Then expressing the argument [n] q D using [n] q gives the formula.

This property can be seen in example 1.3.

Recall that the coefficients of the flag f -vector of a graded poset P are the numbers of strict chains c with a fixed sequence of ranks rk(c). For more on this, see [START_REF] Billera | Linear inequalities for flags in graded partially ordered sets[END_REF] and [19, §3.12].

Lemma 2.8. Let P be a graded poset and consider the height function rk. The polynomial Z P,rk is entirely determined by the flag f -vector of P .

Proof. It follows from formula (7) that Z P,rk is a linear combination of q-Ehrhart polynomials E a whose coefficients are exactly coefficients of the flag f -vector.

Note: one can wonder about the converse implication. Most probably, the flag f -vector should be a finer invariant than the q-Zeta polynomial.

3 Relation with the q-order polynomial

The order polytope Q P of a poset P [START_REF] Stanley | Two poset polytopes[END_REF] is the lattice polytope in N P defined by inequalities 0 ≤ z p ≤ 1 for all p ∈ P and z p ≤ z q if p ≤ q in P . The application (z p ) p → (1 -z p ) p is a bijection between Q P and Q P for the dual poset P .

Let us recall the q-order polynomial L P of a poset P as introduced in [START_REF] Chapoton | q-analogues of Ehrhart polynomials[END_REF]. This is the q-Ehrhart polynomial of the order polytope Q P , with respect to the linear form (1, 1, . . . , 1), i.e. the sum of coordinates. The value of L P at [n] q for n ≥ 0 is therefore

L P ([n] q ) = z∈nQP q p zp . ( 20
)
The degree of L P is #P . Let P be a finite poset. Let J(P ) be the distributive lattice of lower ideals in P under the order of containment. This lattice is graded by the cardinality of the lower ideal. Fix n ≥ 2, and consider a chain e 1 ≤ . . . ≤ e n-1 of elements of J(P ). This is an increasing chain of lower ideals of P . As such, it is characterised by the following data : for each element p of P , let z p be the smallest integer between 0 and n -2 such that p ∈ e zp+1 if it exists and n -1 otherwise. Then these vectors (z p ) p∈P satisfy that z p ≤ z q if p ≤ q in P . This map defines a bijection between chains e 1 ≤ . . . ≤ e n-1 in J(P ) and lattice points in the (n -1)-dilate of the order polytope Q P . Proposition 3.1. There holds Z J(P ),rk (1 + qx) = L P (x).

Proof. Let us compare the values at x = [n -1] q for n ≥ 2. Let us compute the monomial in q attached in Z J(P ),rk (1 + qx) = Z J(P ),rk ([n] q ) to one chain e 1 ≤ . . . ≤ e n-1 in terms of the corresponding lattice point (z p ) p∈P in (n -1)Q P . For the chain, the monomial is q to the power i #e i = i p∈ei 1. By the bijection between chains and lattice points, the exponent of q becomes p∈P (n -1 -z p ) where z ∈ (n -1)Q P . Using the bijection (z p ) p → (1 -z p ) p between Q P and Q P , the exponent of q becomes p∈P z p where z ∈ (n -1)Q P . Summing these monomials over z gives exactly the value at [n-1] q of the q-order polynomial of P .

For example, let P be the poset with three elements a, b and c, such that a is less than b and c. Then the q-Ehrhart polynomial of P is 1 [START_REF] Baumann | Chains in shard intersection lattices and parabolic support posets[END_REF] q [3] q

• (qx + 1) • (q 2 x + q + 1) • ( q 3 + q 2 x + q 2 + q + 1), and the q-Zeta polynomial of J(P ) is

1 [2] q [3] q • x • (qx + 1) • ( q 2 + q x + 1).
4 Ehrhart series and volumes 4.1 P -partitions and q-Ehrhart series for posets

Let us first recall the following classical setting, part of the more general theory of P -partitions, due to Stanley [16, §8].

The q-Ehrhart series of any poset P is the formal power series defined by

L P = n≥0 L P ([n] q )t n , (21)
and can be expressed as a rational fraction

L P = H P (q, t) (1 -t)(1 -qt) . . . (1 -q #P t) , (22)
where H P is a polynomial in q and t with non-negative integer coefficients. This polynomial has a known combinatorial interpretation, using descents for t and major indices for q, as a sum over all linear extensions of the poset P . Let us also introduce the q-volume of a poset, as defined in [9, §4.2].

Definition 2. The q-volume of a poset P is the leading coefficient of the q-Ehrhart polynomial L P times the q-factorial [#P ]! q .

Proposition 4.9 in [START_REF] Chapoton | q-analogues of Ehrhart polynomials[END_REF] gives the following relationship between the polynomial H P and the q-volume of P . Proposition 4.1. The q-volume of P is equal to q ( #P +1

2

) times H P (1/q, 1).

q-Ehrhart series for q-Zeta polynomials

Inspired by proposition 3.1 which identifies the q-Zeta polynomials of distributive lattices with q-order polynomials of their posets of join-irreducibles, it is natural to extend the constructions of the previous paragraph to the general case, for the q-Zeta polynomials of arbitrary posets.

Let (P, h) be a poset endowed with an height function h. In this section, H will denote the maximal value of h.

Let us consider the generating series of values of Z P,h :

Z P,h = n≥0 Z P,h ([n + 1] q )t n , (23)
similar to the q-Ehrhart series (21).

Proposition 4.2. This series can be expressed as a rational fraction:

Z P,h = H P,h (q, t) (1 -t)(1 -qt) . . . (1 -q H t) , (24)
where H P is a polynomial in q, q -1 and t with integer coefficients. The degree of H P,h with respect to t is at most H.

Proof. Let us write Z for Z P,h . Consider the polynomial Z(1 + qx). By eq. ( 3) and lemma 2.2, its value at every [n] q for n ≥ 0 is a polynomial in q with integer coefficients. By proposition B.2, it can therefore be expressed as a sum

Z(1 + qx) = 0≤j≤H c j B j (x), (25)
where c j ∈ Z[q, q -1 ] and B j are polynomials defined in Appendix B by (33). Then one concludes using the generating series from lemma B.5.

By analogy with definition 2, let us introduce the volume associated with the q-Zeta polynomial. Definition 3. The q-Zeta volume of a poset P with respect to the height function h is the leading coefficient of the q-Zeta polynomial Z P,h times the q-factorial [H]! q .

In this context, one has the following general property, for arbitrary posets. Let us choose H -h as the height function on the dual poset P . Proposition 4.3. The q-Zeta volume of P w.r.t. h is q ( H 2 ) times H P ,H-h (1/q, 1).

Proof. The proof is similar to the proof of [START_REF] Chapoton | q-analogues of Ehrhart polynomials[END_REF]Prop. 4.9]. Let us write

Z P,h = H j=0 h j t j H ℓ=0 1 -q ℓ t , (26)
for some coefficients h j in Z[q, q -1 ]. By lemma 2.5 translated into generating series using (23), one has the equality H P ,H-h (q, t) = H P,h (1/q, q H t).

After replacing q by 1/q in the evaluation at t = 1, one gets

H P ,H-h (1/q, 1) = H j=0 h j q -jH .
On the other hand, one can compute directly the coefficient of x H in Z P,h (1+ qx) from (26). For the term j = 0, one can use lemma B.5 and (33) that gives the corresponding coefficient of x H explicitly as

1 [H]! q H j=1 q j = q ( H+1 2 ) [H]! q .
For the terms of index j ≥ 1, one notes that the product by t in the series amounts to replace x by (x-1)/q in the polynomial coefficient, which multiplies the leading coefficient of x H by q -H . The total leading coefficient of

Z P,h (1 + x) is therefore given by   H j=0 h j q -jH   q ( H+1 2 ) [H]! q .
For the leading coefficient of Z P,h (x), one has moreover to divide by q H . One concludes by comparing the obtained expressions.

For example, for the self-dual poset of example 1.5 with H = 2, the series is

Z P,rk = 1 + 2qt (1 -t)(1 -qt)(1 -q 2 t)
and the q-volume is q + 2.

Positivity properties

By proposition 3.1 and the discussion in §4.1, when P is a distributive lattice, then H P,rk is a polynomial in q and t with non-negative integer coefficients.

This positivity property of H P,rk is not true for all posets. Small counterexamples are the posets of example 1.7 and example 1.8. In the case of example 1.8 with H = 2, one gets the numerator -q 3 t 2 + 2 + q + q 2 t -1.

In this section, we will give a sufficient criterion for positivity of H P,rk , in terms of the existence of an R-labelling.

R-labellings

Let L be a set of labels, endowed with an arbitrary relation denoted by

♥ ↔. For p < q in P , a maximal chain p = e 0 < e 1 < • • • < e n+1 = q is increasing 1 if λ(e i-1 , e i ) ♥ ↔ λ(e i , e i+1 ) for all 1 ≤ i ≤ n.
An R-labelling of a poset P by (L, ♥ ↔) is an application λ from the set of edges of the Hasse diagram of P to L such that:

• for every pair of comparable elements p ≤ q in P , there is exactly one maximal chain from p to q.

1 This terminology comes from the case where ♥ ↔ is a partial order relation, which is not assumed here. One may say concatenable of friendly for a better intuition.

For example, the weak order on the symmetric group S 3 has no R-labelling, as both maximal chains from the minimum to the maximum are necessarily increasing.

Remark 5.1. Every EL-labelling, as defined in [START_REF] Björner | Shellable and Cohen-Macaulay partially ordered sets[END_REF], is also an R-labelling, as it satisfies a stronger condition. Supersolvable lattices and upper-semimodular lattices always have an R-labelling. This is proved in [19, examples 3.13.4 and 3.13.5].

Positivity criterion

Proposition 5.2. If P is a bounded and graded poset that admits an R-labelling, then H P,rk is a polynomial with positive integer coefficients.

Proof. The proofs follows from theorem A.1 that gives positivity of the flag h-vector under the given hypothesis and the fact that the coefficient of H P,rk are non-negative linear combinations of the flag h-vector elements, as proved in theorem A.2. Lemma 5.3. If P is a bounded and graded poset that admits an R-labelling, then the q-Zeta volume of Z P,rk is a non-negative q-analogue of the number of maximal chains in P .

Proof. This follows from proposition 5.2 and proposition 4.3. The relationship between the q-Zeta volume and H P,rk involves the dual poset P , but the set of maximal chains is preserved by duality. Proposition 5.2 above applies to several classical posets attached to finite Coxeter groups. The first ones are non-crossing partitions lattices issued from the absolute order [START_REF] Brady | Non-crossing partition lattices in finite real reflection groups[END_REF][START_REF] Brady | K(π, 1)'s for Artin groups of finite type[END_REF][START_REF] Bessis | The dual braid monoid[END_REF][START_REF] Mccammond | Noncrossing partitions in surprising locations[END_REF] and shard-intersection orders [START_REF] Reading | Noncrossing partitions and the shard intersection order[END_REF][START_REF] Reading | Noncrossing partitions and the shard intersection order[END_REF]. These two families are known to be EL-shellable by results of [START_REF] Athanasiadis | Shellability of noncrossing partition lattices[END_REF][START_REF] Petersen | On the shard intersection order of a Coxeter group[END_REF].

The intersection lattice of an essential central hyperplane arrangement is a geometric lattice and is therefore EL-shellable [START_REF] Björner | Shellable and Cohen-Macaulay partially ordered sets[END_REF]Th. 3.1]. This applies in particular to posets of generalized set partitions defined as intersection lattices of reflection hyperplane arrangements.

Hence, all these posets have an R-labelling and non-negative H P,rk .

There is a more confidential family of posets attached to finite Coxeter groups, the parabolic-support posets [START_REF] Bergeron | Posets related to the connectivity set of Coxeter groups[END_REF][START_REF] Baumann | Chains in shard intersection lattices and parabolic support posets[END_REF]. Every interval in this family is shellable and upper-semimodular, hence has an R-labelling [START_REF] Bergeron | Posets related to the connectivity set of Coxeter groups[END_REF]. But these posets are not bounded, hence the criterion does not apply to the full posets. Positivity of H P,rk seems nevertheless to hold for the parabolic-support posets. This remains to be proved and explained.

Another interesting example is given by the root posets of type B, which are not bounded posets either. One can endow these posets with an R-labelling according to the direction (North-East or North-West) of the cover relations. It seems that the coefficients of powers t k in H P,rk are non-negative q-analogues of n k 2 , the type B Narayana numbers.

One can wonder what would be a necessary and sufficient condition on posets for positivity of the coefficients of H P,rk .

6 Value at q = 0 and characteristic polynomial Assume in this section that P is bounded and graded with a unique minimal element 0 and a unique maximal element 1. Let H be the maximal value of the rank function rk on P .

The characteristic polynomial of P is defined as

X P (y) = p∈P µ( 0, p) y H-rk(p) . (27)
This can be expressed using the flag h-vector as follows.

Lemma 6.1. The polynomial

y H X P (1/y) = 1 + H j=1 (β P ([1, . . . , j -1]) + β P ([1, . . . , j])) (-y) j . ( 28
)
Proof. The polynomial on the left is

1 + H j=1     p∈P rk(p)=j µ( 0, p)     y j .
Then, for j ≥ 1, the coefficient of y j is the difference Theorem 6.2. The q-Zeta polynomial Z P,rk has no pole at q = 0 and its value at q = 0 is related to the characteristic polynomial by

Z P,rk | q=0 (1 -y) = y H X P (1/y). ( 29
)
Remark 6.3. It can well happen, when the poset P has no unique minimum, that Z P,h has poles at q = 0, for instance for the poset of example 1.6 and the dual poset in example 1.7.

Proof. Let us denote Z P,rk by Z. Let us start with theorem A.2. One deduces that tH P,rk (q, t) = S⊆{1,...,H-1}

β P (S)t #S+1 q ΣS .
As P is bounded, Z(0) = 0 by lemma 2.3. Therefore the previous expression is the numerator of the generating series of values of Z for n ≥ 0. By the correspondence stated in appendix B, this means that Z P,rk (x) = S⊆{1,...,H-1} β P (S)q ΣS H-#S-1;x H q .

One can then use lemma B.8 to let q be 0 and obtain

Z| q=0 (x) = H-1 j=0 β P ([1, . . . , j])x(x -1) j .
Indeed, all sets S that are not formed of consecutive integers starting at 1 appear with a strictly positive power of q, hence vanish when q = 0. One deduces

Z| q=0 (1 -y) = H-1 j=0 β P ([1, . . . , j])(1 -y)(-y) j .
One can then conclude by an easy comparison with lemma 6.1.

A Flag f -vectors and h-vectors

Let us recall the standard definitions of flag f -vectors and flag h-vectors and then state theorems about their relationship with maximal chains.

Let P be a graded and bounded poset with unique minimum 0 and unique maximum 1. Let rk be the rank function on P , with minimal value 0 on 0 and maximal value H on 1.

For a subset S ⊆ {1, 2, . . . , H -1} of cardinality k, let α P (S) be the number of chains 0 < p 1 < p 2 < • • • < p k < 1 in P such that S = {rk(p 1 ), . . . , rk(p k )}.

The flag f -vector of P is the collection of numbers α P (S), indexed by subsets S of {1, 2, . . . , H -1}.

The flag h-vector of P is the collection of numbers β P (S), also indexed by subsets of {1, 2, . . . , H -1}, defined by Möbius inversion as follows: ↔ is a partial order relation.

= p 0 < p 1 < p 2 < • • • < p H = 1 in P , one can associate its descent set 2 D(M ) = {i ∈ {1, . . . , H -1} | not λ(p i-1 , p i ) ♥ ↔ λ(p 1 , p i+1 )}. ( 31
The action of ∆ q can be translated into an action on sequences of Laurent polynomials as follows.

Abusing notation, let ∆ q be the linear operator acting on sequences (a n ) n∈Z of Laurent polynomials by the formula (∆ q (a)) n = (a n -a n-1 )/q n .

(36) This is the same as (34) when a n = P ([n] q ). Lemma B.6. A sequence of Laurent polynomials (a n ) n≥0 is the sequence of values of an element of A q at q-integers [n] q with n ≥ 0 if and only if it is annihilated by some power of ∆ q .

Proof. Assume first that a n = P ([n] q ) for some element P of A q and for all n ≥ 0. Then because the operator ∆ q acts by sending the basis element B k to B k-1 , iterating ∆ q sufficiently gives the constant sequence 0. Conversely, one proceeds by induction on d where the sequence a is annihilated by the power d + 1 of ∆ q . If the sequence a is identically zero, the statement is obvious.

Otherwise, the sequence ∆ d q (a) is a constant Laurent polynomial c. The sequence a ′ n = a n -cB d ([n] q ) is then annihilated by ∆ d q , hence by induction a ′ (n) = P ([n] q ) for some element P of A q . The statement follows.

For integers a and b, let us introduce the following polynomials

a;x b q = a j=a-b+1 [j] q + q j x [b] q ! . (37)
Their evaluation at x = [n] q is the standard q-binomial coefficient a+n b q . Note that B k (x) = k;x k q for k ≥ 0. These polynomials are useful in the following correspondence.

Proposition B.7. Consider the generating series n≥0

Z n t n = d k=0 h k t k d ℓ=0 (1 -q ℓ t) , (38)
where h k are arbitrary coefficients. Then the Z n are the values at [n] q of the polynomial

d k=0 h k d-k;x d q . ( 39
)
Proof. Starting from lemma B.5 corresponding to k = 0, one can obtain by induction on k the equality

n≥0 d-k+n d q t n = t k d ℓ=0 (1 -q ℓ t)
.

The statement readily follows.

Thee polynomials a;x b q also have a very simple limit at q = 0, once normalized appropriately.

Lemma B.8. For 0 ≤ i ≤ d, the polynomial d-i;x d q q ( i 2 ) has a limit when q = 0 and this limit is 1 if i = 0 and x(x -1) i otherwise.

Proof. One uses the explicit product formula (37) for these polynomials. The statement is clear if i = 0. When i > 0, this is a simple computation term-byterm in the product, by distributing the factor q ( i 2 ) in the appropriate factors to compensate for poles in q.

C q-analogue of incidence algebra

The aim of this section is to explain an interpretation of the values of q-Zeta polynomials as corner coefficients of powers of matrix-like objects, similar to the classical case taking place in the incidence algebra.

The idea of a possible q-analogue of the incidence algebra has been alluded to in the last paragraph of [START_REF] Stanley | Subdivisions and local h-vectors[END_REF].

Let us consider a poset P endowed with the height function h.

Let Inc q (P, h) be the vector space of square matrices A with rows and columns indexed by P , with coefficients in the ring of Laurent polynomials in q, and such that A x,y = 0 if x is not smaller than or equal to y in P . Let D h be the diagonal matrix with coefficients q h(x) for x ∈ P .

The space Inc q (P, h) is endowed with the following bilinear product

(A × q B) x,z = AD h B. (40)
The product × q is associative, with the diagonal matrix D -1 h as unit. This ring is a twisted version of the usual incidence algebra, to which it reduces when q = 1. It can be identified with a sub-ring of the usual ring of square matrices by multiplying all matrices on the left by D h .

Let Z ∈ Inc q (P, h) be the zeta matrix of P , defined by Z x,y = 1 for all relations x ≤ y in P and 0 otherwise. The matrix Z is invertible in Inc q (P, h), because it is upper triangular with invertible diagonal.

Let now P be a bounded poset. Let H be the maximal value of h on P . For n ∈ Z, let Z ×qn be the n-th power of Z in Inc q (P, h).

We will use the ring A q and the operator ∆ q as defined in §B.

Proposition C.1. 7 The sequence of matrices (Z ×qn ) n∈Z is annihilated by ∆ H+1 q . For all n ∈ Z, the value of the q-Zeta polynomial Z P,h ([n] q ) is the coefficient of index ( 0, 1) in the matrix Z ×qn . Proof. For n ≥ 2, the definition of the matrix Z ×qn as a product ZD h ZD h • • • Z implies directly that its corner coefficient is the weighted sum over all chains 0 ≤ e 1 ≤ • • • ≤ e n-1 ≤ 1, where the weight is q j h(ej ) . This is exactly eq. ( 3). This implies that the sequence (Z ×qn 0, 1 ) n≥2 is annihilated by ∆ H+1 q . The same proof applies, for every relation x ≤ y, to the sequence of coefficients (Z ×qn x,y ) n≥2 , which is therefore also annihilated by ∆ H+1 q . The sequence of matrices (Z ×qn ) n≥2 is annihilated by ∆ H+1 q acting by (36). Now consider the sequence ∆ H+1 (Z ×qn ) n≥N for some N ∈ Z. This is ∆ H+1 (Z ×qn ) n≥N = Z ×q(N -2) × q ∆ H+1 (Z ×qn ) n≥2 = 0.

(41)

It follows that the whole sequence (Z ×qn ) n∈Z is annihilated by ∆ H+1 q . Hence every sequence of coefficients of fixed index (x, y) is a polynomial of degree at most H evaluated at [n] q , and in particular the corner coefficient coincides with the q-Zeta polynomial of P .

Note that this statement also implies lemma 2.4.

  p∈P rk(p)≤j µ( 0, p) -p∈P rk(p)≤j-1 µ( 0, p), in which each sum is minus the Möbius number µ S of a rank-selected subposet, for the sets S = [1, . . . , j] and S = [1, . . . , j -1]. Using the relationship µ S = (-1) #S-1 β P (S) ([19, 3.12]), one can then conclude.

  #S-#T α P (T ) and conversely α P (S) = T ⊆S β P (T ). (30) Note that the numbers β P (S) are not obviously non-negative, as their definition by Möbius inversion involves signs.Assume now that P admits an R-labelling λ for the relation ♥ ↔, as defined in §5.1. To each maximal chain M : 0

)

  Theorem A.1(Björner and Stanley). Let P be a bounded and graded poset with an R-labelling. Let S ⊆ {1, 2, . . . , H -1} be any subset the set of heights. The number α P (S) counts maximal chains in P with descent set contained in S. The number β P (S) counts maximal chains in P with descent set S and is therefore non-negative.

	For the proof, see [19, §3.13] or [6, Th. 2.7]. One can check that the proof
	of [19, Th. 3.13.2] works verbatim without the hypothesis that	♥

Again, the terminology comes from the case where ♥ ↔ is a partial order.

Richard Stanley, with my sincere admiration.

The author has made use of SageMath [10] in the preparation of this article.

The author has been supported by the ANR Combiné (ANR-19-CE48-0011) and ANR Charms (ANR-19-CE40-0017).

Some interesting information about this statement can be found in [11, §8.1].

There is a simple relationship between the flag h-vector and the numerator of the q-Ehrhart series as defined in 4.2.

Theorem A.2. The polynomial H P,rk is the sum S⊆{1,...,H-1} β P (S)t #S q (S) (32)

and has therefore non-negative coefficients.

Proof. Introduce formal variables u 0 , u 1 , . . . , u H . For a chain K in P \ { 0, 1}, let u K be the product of u rk(p) over elements p in the chain. For a subset S ⊆ {1, 2, . . . , H -1}, let u S be the product of u i over elements of S.

According to [START_REF] Stanley | Enumerative combinatorics[END_REF]Ex. 3.67(b)], one has the formula

where the sum over K is running over all chains in P \ { 0, 1}. By extending the sum over K to all chains in P , one gets

with just two additional factors in the denominator accounting for the multiplicities of 0 and 1 in the chain. By specializing every u ℓ to q ℓ t, one gets the equality

Comparing with the definition of H P,rk as a numerator in proposition 4.2, one obtains the expected formula.

B Rings of polynomials with integer Laurent values

Let A q be the sub-ring of Q(q)[x] made of polynomials P such that P ([n]) q ∈ Z[q, q -1 ] for all n ≥ 0. This is an analogue of the ring of integer-valued polynomials.

For k ≥ 0, let B k be the polynomial defined by

Because the values B k ([n] q ) for n ≥ 0 are standard q-binomial coefficients, every B k belongs to A q . We will use the convention that B -1 = 0.

Let ∆ q be the operator acting on polynomials in x defined by ∆ q (P )(x) = P (x) -P ((x -1)/q) 1 + (q -1)x .

The right hand side is well defined as a polynomial because the numerator has a root at x = 1/(1 -q). Assuming moreover that P ∈ A q , this formula implies that the values ∆ q (P )([n] q ) for n ≥ 1 are Laurent polynomials in q with integer coefficients.

Lemma B.1. The operator ∆ q maps B k to B k-1 for all k ≥ 0.

Proof. It is enough to prove that it holds for the value at every q-integer [n] q for n ≥ 1. This in turn follows from the classical formula

for q-binomial coefficients.

Proposition B.2. The polynomials B k form a basis over Z[q, q -1 ] of the ring A q .

Proof. These polynomials are linearly independent, as B k has degree k. It remains to prove that they span A q over Z[q, q -1 ]. The proof is by induction on the degree. This is clear in degree 0 as B 0 = 1. So let P be a polynomial of degree d > 0 in A q . Then Q = ∆ q (P ) is a polynomial with values in Z[q, q -1 ] for n ≥ 1. Moreover, if the leading coefficient of

Iterating d times the operator ∆ q , one reaches a constant polynomial, with values in Z[q, q -1 ] for n ≥ d and with leading term c. Therefore c ∈ Z[q, q -1 ] and the polynomial P -cB d belongs to A q , with degree at most d -1. By induction, this difference is a linear combination over Z[q, q -1 ] of B k for 0 ≤ k ≤ d -1. Hence P is a linear combination over Z[q, q -1 ] of B k for 0 ≤ k ≤ d.

Corollary B.3. Let P be an element of A q . Then P ([n] q ) ∈ Z[q, q -1 ] for all n ∈ Z. The operator ∆ q acts on A q .

Proof. The first statement holds because this is true for all basis elements B k by an easy computation. The second one follows from the first and (34).

Proof. The special case k = 0 is immediate. It is enough to prove this equality when evaluated at [n] q for n ≥ 1. This reduces to a standard property of q-binomial coefficients.

The following statement can be found for instance in [START_REF] Chapoton | q-analogues of Ehrhart polynomials[END_REF]Lem. 4.7].