{'fraction_non_alphanumeric': 0.027895944517492974, 'fraction_numerical': 0.00781189573826282, 'mean_word_length': 4.744668246445498, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Detection and characterization of instrumental transients in LISA Pathfinder and their projection to LISA

The LISA Pathfinder (LPF) mission succeeded outstandingly in demonstrating key technological aspects of future space-borne gravitational-wave detectors, such as the Laser Interferometer Space Antenna (LISA). Specifically, LPF demonstrated with unprecedented sensitivity the measurement of the relative acceleration of two free-falling cubic test masses. Although most disruptive non-gravitational forces have been identified and their effects mitigated through a series of calibration processes, some faint transient signals of yet unexplained origin remain in the measurements. If they appear in the LISA data, these perturbations (also called glitches) could skew the characterization of gravitational-wave sources or even be confused with gravitational-wave bursts. For the first time, we provide a comprehensive census of LPF transient events. Our analysis is based on a phenomenological shapelet model allowing us to derive simple statistics about the physical features of the glitch population. We then implement a generator of synthetic glitches designed to be used for subsequent LISA studies, and perform a preliminary evaluation of the effect of the glitches on future LISA data analyses.

The LISA Pathfinder (LPF) mission succeeded outstandingly in demonstrating key technological aspects of future space-borne gravitational-wave detectors, such as the Laser Interferometer Space Antenna (LISA). Specifically, LPF demonstrated with unprecedented sensitivity the measurement of the relative acceleration of two free-falling cubic test masses. Although most disruptive nongravitational forces have been identified and their effects mitigated through a series of calibration processes, some faint transient signals of yet unexplained origin remain in the measurements. If they appear in the LISA data, these perturbations (also called glitches) could skew the characterization of gravitational-wave sources or even be confused with gravitational-wave bursts. For the first time, we provide a comprehensive census of LPF transient events. Our analysis is based on a phenomenological shapelet model allowing us to derive simple statistics about the physical features of the glitch population. We then implement a generator of synthetic glitches designed to be used for subsequent LISA studies, and perform a preliminary evaluation of the effect of the glitches on future LISA data analyses.

I. INTRODUCTION
LISA Pathfinder (LPF) was launched on December 3 rd , 2015 as a technological demonstrator for the Laser Interferometer Space Antenna (LISA), a future spacebased gravitational-wave observatory [1]. The LPF mission measured the relative acceleration between two cubic free-falling test masses (TMs) housed within the spacecraft (SC), up to a precision of 30 fm s −2 Hz −1/2 for mHz frequencies, as depicted by the orange spectrum in Fig. 1.
LPF demonstrated the feasibility of LISA free-fall precision requirement via the assessment of parasitic differential acceleration ∆g(t) between the two free-floating TMs, which constitutes the main science measurement of the LPF mission. This performance reached the required free-fall levels for LISA [2,3], and was far more sensitive than LPF's precision target. Several incremental corrections and calibrations were necessary to reach this better-than-expected goal [4], along with the characterization of various subsystems. The satellite was successfully operated until its planned decommissioning on July 18 th , 2017.
The TMs differential acceleration data ∆g(t), once cleaned of the inertial effects, were found to be affected by spurious non-gravitational forces of instrumental origin, which degraded LPF noise performance and limited its sensitivity. These perturbations, often referred to as glitches, were transient in nature. To produce the results in Ref. [3], they were subtracted from ∆g(t) data in order to recover the full LISA-like sensitivity, allowing to go from the green to the orange periodogram in Fig. 1. The removal of these transients was particularly critical to obtain the full sensitivity at the lowest frequencies. Should they arise in LISA data, their presence could impact the characterization of gravitational-wave sources. Adequately accounting for them in the data modeling is therefore crucial both for noise and signal estimation.
The physical nature of these spurious events is still unknown. They manifested as uncorrelated non-periodic events during ordinary runs, occurring with varying intensity and duration. The glitch phenomenology presented here appears to be connected to a series of physical counterparts that are currently under investigation and are outside the aim of this work. Hypotheses regarding their physical origin are under analysis and require further simulation and eventually on-ground testing at torsion pendulum facilities. Note that these unmodeled perturbations in the differential acceleration measurements are different from micrometeorids impacts considered by the authors of Ref. [5], which are detected with satellite position telemetry.
Despite the lack of a physical understanding of their origin, glitches can be modeled mathematically. In this work, we adopt an approach to allow for an automated and systematic glitch detection. We use generic basis functions called shapelets [6] to form a template bank and detect events through matched filtering. Here, we only consider TM differential acceleration measurements, and we detect events without any assumption on their physical nature. The objective is two-fold: i) provide a arXiv:2112.07490v1 [gr-qc] 14 Dec 2021 FIG. 1. Amplitude spectral density ( √ PSD) of parasitic differential acceleration of LPF TMs as a function of frequency for the best performance noise run, with (green) and without (orange) spurious signals.
mathematical and comprehensive study of glitch population that can be informative for LISA and ii) obtain a distribution of event parameters serving as an input for LISA simulation studies, as well as data challenges [7].
The article is outlined as follows. In Section II, we start by describing the measurements used in this analysis. We then introduce the shapelet model and the detection method in Section III. In Section IV we present the results of the glitch search and evaluate its performance, before analyzing the statistics of glitch parameters in Section V. We show in Section IV D that these parameters can also be estimated through Bayesian analysis from matched filtering outputs, for an optimal modeling and removal. Finally, based on the statistics we obtain, we describe the design of a glitch generator readily usable for LISA studies in Section VI.

II. ANALYZED MEASUREMENTS
A. Description of the data The differential acceleration ∆g(t) is derived from the most sensitive measurement of the TMs position, performed via the differential interferometer o 12 [8,9], as where the ω 2 i terms are the spring-like constants of noisy forces acting on the TMs, due to electrostatic, magnetic, and local gravity effects, while g c (t) are the commanded forces, known up to a multiplicative amplitude calibration constant.
Equation (1) describes the dynamics along the x-axis, but does not include any inertial effects caused by cross-coupling with other degrees-of-freedom. Thus the understanding and accurate calibration of ∆g(t) required a series of dedicated experiments, designed to identify the dynamical relations between the three-body system constituted by the TMs and SC, and the various environmental noises contributing to the overall noise budget [4,10].
After subtraction of the inertial effects, ∆g(t) includes the sensitive x-axis dynamics of LPF TMs free of the inertial forces acting on the system, allowing for evaluation of the LPF noise performance on noise-only measurements. Noise-only measurements were acquired during mission operations, in-between the different calibration and system identification experiments, lasting up to two weeks at a time, in order to benchmark LPF performance.
The analysis in this paper only examines the noise-only ∆g measurements containing spurious signals, which are publicly available on the LPF legacy data archive [11]. We select continuous time series segments where no signal injections took place, and we use the developed detection method to identify the transient signals in the noise (see Section III below). The actual time segments are listed in Appendix A.

B. Target instrumental perturbations
The data exhibit spurious signals of varying duration and amplitudes. These artefacts are distinguishable from stationary Gaussian noise, and can be described by a transient process with short time duration compared to the overall measurement.
We observe two types of glitches. The majority include a sharp rise followed by a slow exponential decay back to the previous acceleration trend. The initial slope can be a positive or negative acceleration. Adopting the gammaray burst terminology [12], we will refer to these glitches as fast rise and an exponential decay (FRED) events. An example is shown in Fig. II B. The other common type of glitches, which occur less frequently, have a sine-Gaussian shape characterized by two sequential excursions with opposite signs, as shown in the example of Fig. II B). Although this is beyond the scope of this work, ongoing analysis suggest that the two glitch types are connected to different physical phenomena [13,14].
During nominal mission operations, glitches appeared stochastically with a characteristic time interval. However, during the course of LPF mission, the system was cooled down in order to further investigate the effects contributing to the overall noise budget. Unexpectedly, this increased the glitch event rate substantially. We will henceforth refer to these measurements as cold runs, and we will separately investigate them in addition to the ordinary noise runs.

III. DETECTION METHOD
A. Penalized likelihood model In order to make an inventory of transient signals in LPF ∆g data, we need a flexible framework to detect the corresponding excess of power. This can be done by decomposing glitch signals into a dictionary of functions where they are sparse, i.e. only a reduced number of nonzero coefficients is necessary to accurately represent the signal.
The general model of the measured data can be written as follows: where y is the observed time series of size N sampled at frequency f s , s is the glitch signal, and n is the realization of a stochastic random variable assumed to follow a zeromean, stationary Gaussian distribution. Formally, we design the detection pipeline by finding the maximum of a penalized likelihood function whose logarithm can be written as where p θ (y|θ) is the unpenalized likelihood function of the model parameters θ, equal to the Gaussian probability density function in this problem. The second term pen λ (θ) is the penalty function, which controls the thresholding of the function amplitudes. The vector θ is the vector of parameters to be estimated, and λ is the regularization parameter, controlling the model quality in representing the data and the sparsity level. For a Gaussian zero-mean distribution, the unpenalized likelihood writes where Σ is the covariance matrix of the noise vector n, x † refers to the Hermitian conjugate of any vector or matrix x and |X| is the determinant of any matrix X. Note that in this formulation, the model parameter vector includes both signal parameters θ s and noise parameters θ n , and can therefore be written as θ ≡ (θ s , θ n ). We further assume that the noise is wide-sense stationary, with power spectral density (PSD) S n (θ n , f ), so that its covariance can be approximated by a circulant matrix that is diagonalizable in the Fourier basis: where F is the discrete Fourier transform matrix with elements F kn = e −2jπkn/N and Λ is a diagonal matrix whose diagonal elements are proportional to the PSD Λ kk = N f s S n (θ n , f k ) In the following, we choose a penalty function equal to the l 0 pseudo norm of the coefficient vector, which imposes the sparsest constraint on the representation, and is well adapted to sharp features: The l 0 -norm of a vector gives its number of non-zero coefficients, thus minimizing it enforces sparsity.

B. Modeling transient signals using shapelets
As we do not postulate any physical process underlying the glitch phenomena, we choose to adopt a representative phenomenological model. Given previous observations, this model should have the following characteristics: it should allow for i) a sharp rise, ii) an exponential decay and iii) one or a few oscillations. Most of wavelet functions, including Gaussian-shape sinusoids, are not well suited to that kind of transients because of their symmetric shape. Asymmetric functions defined for positive times only are preferable. In addition, the model should be general enough to describe the various types of glitches that have been cataloged.
A good option that meets these requirements is the exponential shapelet function described in Ref. [6]. This is a collection of exponentially damped functions that are eigen-wavefunctions of the normalized 1D-hydrogen atom. They can be defined in the time domain as follows: where c n = (−1) n−1 n − 3 2 is a normalizing factor ensuring that the quadratic sum of the components is equal to 1, and h + (t) is the Heaviside function which is equal to 1 for t ≥ 0, 0 otherwise. A glitch perturbation is then modeled by a finite linear combination of shapelets: The parameters to be estimated for each shapelet component are therefore the scale parameter β, which acts as a characteristic damping time and hence is related to the glitch duration, the arrival time τ and the amplitude α. We denote as θ s ≡ (β, τ, α) the vector gathering the shapelet parameters.

C. Matching pursuit algorithm
A common approach for signal detection in gravitational-wave data analysis is to perform a matched filtering of the data with a template bank, correlating the observed signal with a parametric waveform model. This method is equivalent to calculating the ratio between two likelihoods: the likelihood in the numerator corresponds to the hypothesis that a signal is present, and the one in the denominator corresponds to the null hypothesis (absence of signal).
However, when the signal is assumed to be the superimposition of several non-orthogonal waveform atoms, direct matched filtering can be cumbersome. Instead, here we use an iterative version called matching-pursuit [15], where matched filtering is applied at each iteration to the residuals of the fit at the previous iteration. The final estimation is given by the sum of the signals estimated individually. It can be shown that this algorithm approximately solves the maximum penalized likelihood problem where L (α) is the penalized log-likelihood given in Eq.
(3) with fixed covariance parameter θ n (assumed to be known) and with penalty function given by Eq. (6).
Let us consider iteration i, and let us label r (i) the residuals from the previous iteration. We assume that we already found i − 1 components fitting the signal, and we now want to find the next significant element, if any, to add to the sum in Eq. (8). We are looking for the value of θ s that best fits the signal left in r (i) . This is done by matched filtering the residuals with the waveform model h (t, θ) given by i.e., finding the set of parameters τ, β that maximize the output signal-to-noise ratio (SNR) where we defined the scalar products ·|· as The symbolx represents the discrete Fourier transform (DFT) of any vector x, defined as where the normalization depends on the time series size N and the sampling frequency f s , and is chosen such that the square modulus ofx k is homogeneous to a power spectral density. If we denoteψ n (f ) the continuous Fourier transform of the shapelet function in Eq. (7), we can writeh We now fix the shapelet order to n = 1. This choice greatly simplifies the model and hereby lightens the computational cost of matched filtering. The counterpart is that we may need more shapelet components to fit the rare events with multiple oscillations. We approximate the DFT of the shapelet function by its continous Fourier transform evaluated at discrete Fourier frequencies f k , as given in Ref. [6]: Plugging Eq. (14) into Eq. (12) yields Note that we choose the grid of arrival times τ m = m/f s to be the same as the time series sampling cadence, then Eq. (16) can be efficiently computed with an inverse Fast Fourier transform (FFT) algorithm.

D. Detection and false alarm probability
Here we derive the detection threshold used in the matched filtering process. In gravitational-wave data analysis, a common way to do that is through the Fstatistics, which is the likelihood in Eq. (4) maximized with respect to the parameters on which the model depends linearly, called extrinsic parameters [16]. In our problem, the extrinsic parameters are the amplitudes of the shapelet functions α i . Maximizing the logarithm of Eq. (4) with respect to α yields where θ intr ≡ (β, τ ) denotes the intrinsic parameters, i.e. all other parameters different than the extrinsic ones. The scalar product h|y involves a weighted sum of the data DFT componentsỹ k which, in the absence of signal, follows a zero-mean Gaussian distribution. Therefore the scalar product follows a zero-mean Gaussian distribution with covariance h|h −1 , and 2F y is a chi-squared distribution with 2 degrees of freedom [17]. The cumulative probability density of F y is where Γ r is the regularized gamma function, given by where Γ is the gamma function. In practice, we will calculate the F-statistics for a given grid of values of the arrival time τ and the damping parameter β. Assuming that the different variables F i ≡ F (τ i , β i ) are uncorrelated, the probability that neither of them will exceed the threshold F 0 is where q is the number of explored parameters. The total false alarm probability is therefore If we impose a total false alarm probability α, then the threshold must be such that For large q, we can approximate the last equality by We impose to have a α = 0.01% false alarm probability, meaning that there is 0.01% chance that the detection could be triggered by Gaussian random noise only. The corresponding threshold for a 2.5-day segment sampled every 10 seconds with a grid of 80 values for β is nearly SNR = √ F 0 ∼ 5. Adopting this value, we chose to focus on the robustness of the detection over its completeness, in particular to be conservative with respect to power spectral density estimation errors.

IV. APPLICATION TO LPF NOISE RUNS
The matching pursuit method described in the previous section is applied to 56 LPF ∆g noise-only measurements acquired between March 2016 and July 2017, lasting 1 to 10 days. Some of them contain one or several transient events. In this section, we describe the tuning of the search and its outputs.

A. Pre-processing
Prior to running the algorithm, we divide each noise run in segments of T seg = 2.5 days (or smaller, when the noise measurement run itself is shorter). This segmentation allows us to bound the memory cost of the matched filtering, while encompassing the longest events in the search window. After computing the SNR as a function of the arrival time, we crop the SNR time series by T seg /8 = 7.5 hours at the start and at the end, to avoid edge effects (see, for example, [18]). As a result, there is an overlap of 2/8 = 25% between two consecutive segments.
Afterwards, we perform a raw estimate of the noise PSD, by fitting a spline to the log-periodogram. This fit is made by taking the median of the PSDs estimated on 8 subsegments, improving the robustness against glitches. This provides a first estimate for the matching pursuit run, during which the PSD estimate will be updated using the residuals of the fit at each iteration. The parametrization of the search is described below.

B. Parametrization of the search
The parameter grid that we explore is evenly sampled in logarithmic space of β, with 80 values in the interval [0.1, 50000] s. The lower bound corresponds to the sampling time, and the upper bound corresponds to about twice the duration of the longest glitch. In addition, we evenly sample arrival times at the same rate as the input time series (0.1 second). The order of the shapelet model is fixed to n = 1 to mitigate the computational cost of the search. Since a single shapelet of order 1 may not be sufficient to account for one glitch shape, we allow for several shapelet elements to describe a single event. For example, sine-Gaussian events are usually modelled with two elements of the basis. Note that a new shapelet component is added to the model only if the associated SNR in the residuals of the previous iteration exceeds the detection threshold, thereby preventing overfitting.
The iterations of the matching pursuit algorithm stop when there are no more signal in the residuals that exceeds the SNR threshold. The estimated glitch signal is then the sum of all detected components, as given by Eq. (8).

C. Search results
At each iteration, the matching pursuit algorithm computes the SNR time series from the residuals of the previous iteration, using the current estimate of the noise PSD. In Fig. 3, we present the SNR values, maximized over β, computed at the first iteration of the third search segments in the noise run covering the period between February 13 th 2017 and March 2 nd 2017 (or noise run number 66). We choose this excerpt because it corresponds to the best noise performance achieved [3], and we zoom in the interval where two events of different amplitudes are clearly visible. For these events, the matching pursuit detections coincide with events recorded in the legacy catalog made to produce the data in Ref. [3]. In this earlier analysis, glitch detection was performed visually using low-pass filtered data. Then, glitches were subtracted using ad hoc models based on exponential functions or filtered impulses [13]. This comparison confirms our ability to automatically detect transients that were manually highlighted.
We plot in Fig. 4 two examples of spotted events. The first one is a FRED event which corresponds to the SNR peak exhibited in the upper panel of Fig. 3. We also plot a sine-Gaussian event occurring earlier in the run, corresponding to the peak shown in the lower panel of Fig. 3.
In addition, to cross-check the arrival times with the events spotted in the best performance noise run of February 2017, in Fig. 5 we compare the shapelet fit residuals with the de-glitched data that are shown in Ref. [3]. We plot the matching pursuit shapelet residuals (red) against Ref. [3]'s data (orange). We observe a slight excess of power at low frequency, below 0.2 mHz. This difference is due to the grid of damping parameter β used in the matching pursuit search which is not fine enough to provide the best fit of the glitch signal, although it is sufficient to yield a reliable detection. This is particularly true for long-lived glitches where the maximum likelihood value lies between two grid elements. To improve the fit, we can refine the parameters estimation by sampling their posterior distribution using a Markov chain Monte-Carlo (MCMC) technique, as described in Section IV D. Starting with the matching pursuit outputs, we can quickly converge toward a more optimal value. By performing such a refinement on the longest events, we obtain the black residuals in Fig. 5 which are close to Ref. [3]'s data within error bars. Note that this comparison only provides a consistency check between two deglitching approaches, since we do not have access to pure noise residuals.

D. Refined parameter estimation
In this section, we check the results of the matched filtering process on selected events by performing a refined parameter estimation. To this aim, we use a Bayesian approach by sampling the posterior distribution of the glitch parameters: where p(x|θ) is the likelihood function, p(x) is the evidence and p(θ) is the prior distribution of the model parameters. We estimate the posterior distribution for glitch parameters θ s = (α, β, τ ). The prior distribution for all glitch parameters are chosen to be uniform around the values output by the matched filtering. To probe the posterior, we use ptemcee [19,20], a stochastic sampling method based on MCMC with parallel tempering and linear transformations among chain states as proposals. We choose the longest glitch present in the noise run 66 as an example (see table I), as this event has a significant impact on the low-frequency part of the residual spectrum shown in Fig. 5. After sampling, the resulting posteriors for this investigation are shown in Fig. 6. Their marginal distribution follows closely the Gaussian probability density, and becomes stationary around the input values provided by the matching pursuit algorithms, which were θ s = 5.6 × 10 −13 ms −2 , 1.096 × 10 3 s, 3.6 × 10 5 s .
Not surprisingly, the arrival time is the parameter best  spotted by the matching pursuit detection, because the grid is as fine as the data sampling rate. The damping time grid being coarser, the Bayesian estimation allows one to better fit β to an intermediate value. This has also an impact on the amplitude estimate α, since it is correlated with β. The related improvement is visible on the residual plot shown in Fig. 5, where the noise level is lower at frequencies below 0.2 mHz. This result advocates for including a glitch model as part of LISA global fit, i.e. the full characterization of all detectable sources in the data [22]. A glitch detection and subtraction as a pre-processing phase might not be enough to mitigate their impact on the noise PSD estimate and on the bias of gravitational waves (GW) source parameters.

V. ANALYSIS OF GLITCH PARAMETER STATISTICS
A. Description of the observed population In this section, we analyse the distribution of glitch parameters among all detected events. Since some glitches are described by more than one basis component, we assume that flagged arrival times occurring within less than 5 seconds belong to the same event. We represent such multi-component glitches by an equivalent triad (α eq , β eq , τ eq ) corresponding to the component of maximum SNR among the sum of components describing the event. Note that due to the normalization of the shapelet function in Eq. (7), the effective amplitude of a glitch is proportional to α/β. Therefore the parameter α eq corresponds to an equivalent transferred impulse α eq = ∆v eq expressed in units of velocity (m s −1 ). This quantity is important as it is directly related to the event SNR.
We gather in Fig. 7 the histograms of the parameters from all detected events. We separate the cold runs statistics (which represent 20% of the number of events) from the ordinary runs statistics. The upper left panel shows the histogram of time intervals between glitches, which follows an exponential distribution. We estimate a glitch rate of λ ≈ 5 · 10 −5 s −1 (about 4 events per day) for ordinary runs, and a rate about ten times larger for cold runs.
The histograms in the upper right panel exhibit transferred impulses mostly ranging between 10 −15 and 10 −13 m s −1 . Note that the lower bound is likely to be due to the observation bias, i.e. the impossibility to detect events below the noise level. For example, for the estimated PSD in noise run 66, for β ≈ 1 minute and a SNR detection threshold of 5, events with amplitudes below 2 × 10 −16 m s −1 are not detected.
The lower left panel represents the damping times distributions, which show that most events last less than one minute, with rarer events lasting a few hours. The peaks visible at the leftmost edge of the damping time histograms are due to events whose timescales are either close or below the sampling time, resulting in Dirac-like signals.
We also plot the statistics of the number of shapelet components needed to describe single events. In most cases, glitches have FRED shapes and one component is enough to describe them. However, some glitches require two or more components for their fit, like the sine-Gaussian signal shown in the lower panel of Fig. 4.
To examine correlations among the glitch waveform parameters, we plot in Fig. 8 the joint distribution of β eq and ∆v eq in the form of a corner plot. We can crudely distinguish two populations of events, especially when considering the ordinary runs (gray). The first category gathers short-duration events with damping times below one second (lower contours in the figure), while the second category includes longer events, with a peak around one minute (upper contours). The impulses of short events are on average slightly larger than the long ones.

B. Sampling from the glitch parameter distribution
For future instrument and data analysis characterization efforts, it will be important to be able to generate as many glitch events as needed for LISA simulations. To this aim, we developed a method to sample the 2-dimensional distribution of the glitch amplitudes and damping times, to enable the corresponding LISA waveform generation.
Drawing new samples representative of the LPF glitch population requires sampling the complicated twodimensional distribution of the glitch parameters (α, β) presented in Fig. 8, and sampling the times of arrival τ . The approach that we adopt to draw samples from the desired distribution follows the idea that one can define an invertible map from the simple normal distribution N (0, I) to the desired distribution. This can be formal- ized using the following change of variable equation: where the parameters of the glitch are defined as θ = (α, β), p(θ) is the distribution in question and z(.) is the map that we want to estimate. The map can be estimated using machine learning techniques called normalising flows. We used a particular implementation called neural spline flows [23]. We sample the glitches arrival times using the exponential distribution of the intervals between ordinary run events that we fitted from the upper left histogram in Fig. 7.

VI. EXTRAPOLATING THE GLITCH DISTRIBUTION FROM LPF TO LISA
Our analysis of LPF glitches can serve to inform a tentative assessment of the impact of spurious instrumental transients on LISA data processing. Another objective is to develop the tools needed to test data analysis techniques which will be robust to their presence. To do this assessment, we need a process to translate LPF measurement into equivalent LISA simulations.
In LPF, the presence of glitch signals could not be detected in any measurement other than the most sensitive channel that was the differential interferometer measuring the position of one TM relative to the other. If we consider the same scenario in LISA, we have to assume that glitches may only be visible in the most sensitive channels as well, which are time-delay interferometry (TDI) observables [24,25]. TDI variables allow for the cancellation of laser frequency noise and are obtained from a post-process performed on ground. Besides, data analyses are generally performed in units of relative laser frequency deviations, also referred to as fractional frequency. Converting LPF data into equivalent relative laser frequency deviations requires integration with respect to time of the analytical shapelet model. Then, it needs to be projected to LISA arms and transformed through the TDI algorithm.
This approach was implemented in LISA Glitch [26], a Python module easily interfaced with the rest of the simulation tool chain. LISA Glitch is able to generate glitch signals of various shapes, including the shapelet model introduced in Sec. III B. These glitch signals are associated with one injection point and written to a glitch file. This glitch file is then read by the instrument simulator and glitch signals are injected into their respective injection points during the simulation.

A. Derivation of relative frequency deviations induced by glitches
To adjust the glitch shapes to the commonly used units we first integrate the acceleration perturbation ∆a into a fractional frequency signal dν/ν 0 as Plugging the shapelet model into the above equation yields (see Appendix B) Graphs of the original glitch in differential acceleration and the integrated glitch in fractional frequency are shown in Fig. 9.

B. Time-delay interferometry response to glitches
A thorough investigation of glitch propagation would require introducing them at different simulated injection points, propagating the signals through the LISA dynamics and observation of the output in different interferometers. We postpone such an in-depth study for the future, as the development of LISA dynamics simulations is still in progress at the time of the writing of this paper. In the absence of dynamics, we inject glitches directly on the test-mass interferometer or on the inter-spacecraft (science) interferometer. We look at how the resulting output in the TDI channels will differ based on the place of injection. In reality the appearance of the glitches in one of the interferometers will depend on their frequency content and on the Drag Free and Attitude Controllers (DFACS) [27]  We consider the triple quasi-orthogonal combinations A, E and T which are typically used in the literature [28]. The advantage of these combinations is that they exhibit quasi-orthogonal noise and T is a "null" combination which suppresses most of the GW signal at low frequency. To propagate the glitches through LISA data streams, we use the pyTDI package [29] which computes TDI combinations from raw interferometer measurements. We show an example in Fig. 10. We see that the smooth step function shape at the interferometer level in Fig. 9 is turned into sharper features in the TDI combinations, which act as a differentiation on a duration of a few inter-spacecraft delays (about 8 seconds). We inject the glitch signal into one single science interferometer out of six. Note that in reality, the reaction of the spacecraft drag-free system control loop may result in occurrences of the glitch signal in more than one interferometer. More accurate and realistic propagation of glitches in LISA is the subject of an on-going study.
C. Glitch SNR in LISA Using the simulation described in the previous section, we can compute the glitch TDI response, and hence the SNR in simulated LISA data. Although this computation relies on the assumptions we made regarding the injection point, it can provide us with a good order of magnitude estimate. In the following, we choose to treat glitches as a perturbation on one of the TMs.
We have to assume a noise model for the secondgeneration pseudo-orthogonal TDI channels, reflecting the main pathlength noises in the measurement. We use the model provided by the LDC software [7], based on the LISA science requirement document [30].
First, we compare the glitch power in the TDI combinations A, E, and T with the theoretical PSD of the stochastic instrumental noise. We choose an observation time of 10 5 s, which encompasses the duration of the longest glitch, and compute the empirical PSDs of all detected events. We plot the median and the standard deviation of this PSD distribution in Fig. 11, along with the noise PSD.
We see that the median of the glitch power lies below the mean noise level in the LISA frequency band. However, the spreading of the power across a large bandwidth implies large SNRs. Integrating the distribution over all frequencies yields the events SNRs in TDI channels: wheres i (f ) is the Fourier transform of the glitch response in TDI i, and S i (f ) is the instrumental noise one-sided PSD in channel i. We restrict the integration to frequencies greater than f min = 10 −5 Hz and lower than f max = 1 Hz, which are the boundaries of the observatory sensitivity goal, as stated in the LISA Science Requirement Document [30]. Fig. 12 shows the distribution of TDI SNRs derived from LPF events detected in ordinary (gray) and cold (blue) runs. It exhibits values that range between 10 −2 to 10 4 , with 50% of events having SNRs larger than 10.
The maximal values are large with respect to typical GW sources, which can reach a few hundreds for galactic binaries to a few thousands for supermassive black holes. However, glitch signals are very short in time so they can easily be isolated in time-frequency representations. This very preliminary analysis suggests the need for an adequate processing to mitigate the impact of glitches on science performance.

VII. CONCLUSION
In this paper, we investigated the presence of transient noise artefacts in the LPF noise runs. We developed a method to detect these glitches based on the matched filtering of shapelets, a family of adapted mathematical functions. Each shapelet component is parametrized by an amplitude and a decay time, which can be interpreted as an equivalent transferred impulse amplitude and duration, respectively.
We characterize the glitch population by applying the detection scheme to 56 differential acceleration measurements, from which we obtain a distribution of amplitudes, damping times, and arrival times. Then, we develop a machine-learning method that allows one to draw an arbitrary number of events from this distribution. We explore possible consequences for LISA by making assumptions on glitch injection points, i.e. on the system components from which they may originate. Based on these assumptions, we sketch a preliminary process to compute equivalent perturbations in LISA TDI data. We find that half of the detected events have a significant SNR when extrapolated to LISA measurements. This result makes it all the more necessary to study in depth the impact of transient disturbances on LISA performance.
The tool that we develop constitutes a first building block to inform studies assessing the impact of glitches on LISA scientific performance, and develop adapted data analysis methods to mitigate it. The outputs of this study are currently used in the LISA Data Challenge dataset codenamed Spritz [7], where we simulate LISA data corrupted by LPF-like glitches. Investigating the effect of glitches on the detection and characterization of GW sources is paramount for the LISA mission development. While some studies already tackled the problem of distinguishing instrumental glitches from GW bursts [31], more realistic simulation campaigns including LPF-like signals and dynamical instrumental response will be carried out in the near future.