# Multi-MeV electrons from above-threshold ionization of the neon K-shell

A. Yandow  Center for High Energy Density Science, The University of Texas at
Austin, 2515 Speeday Stop C1600, Austin, TX 78712  Lawrence Livermore National
Laboratory, Livermore, CA 94551<EMAIL_ADDRESS>T. N. Ha  C. Aniculaesei  H.
L. Smith  C. G. Richmond  M. M. Spinks  H. J. Quevedo  S. Bruce  M. Darilek
C. Chang  D. A. Garcia  E. Gaul  M. E. Donovan  B. M. Hegelich  T. Ditmire
Center for High Energy Density Science, The University of Texas at Austin,
2515 Speedway Stop C1600, Austin, TX 78712

###### Abstract

We present measurements of integrated electron energies produced by above-
threshold ionization (ATI) of neon in a laser field with intensity exceeding
$10^{20}$ W/cm2. We observe electrons with energy exceeding 10 MeV ejected in
the laser forward direction above a threshold intensity of 2 $\times 10^{20}$
W/cm2. We compare to ATI models using both tunneling (ADK-PPT) and barrier
suppression ionization and observe the onset of ATI at a higher threshold
intensity than predicted by these models.

Suggested keywords

††preprint: AAPM/123-QED

Above-threshold ionization (ATI) is a fundamental response of an atomic system
to an intense flux of photons. The first experimental evidence of ATI showed
the absorption of seven photons in a six-photon multiphoton ionization pathway
of xenon [1]. As near-infrared laser intensity increases beyond $10^{14}$
W/cm2, ATI can be well-described by a quasi-classical two-step model. The
ionization process can be described by the Ammosov-Krainov-Delone and
Perelomov-Popov-Terent’ev (ADK-PPT) model [2][3], where the photons sum
coherently to create a strong electric field which liberates electrons through
a quasi-static tunneling process. The ATI electron is then “born” into the
laser field with initial conditions consistent with the ADK-PPT tunneling
model, and the absorbed laser energy can be found by integrating the classical
Lorentz force equations. The ADK-PPT ionization rate model has been well-
validated by measurements of ion charge states produced in the laser focus at
nonrelativistic intensities [4]. Extension of these experiments to intensity
above $10^{20}$ W/cm2 has been elusive even with advances in laser technology
due to the limited repetition rate of high-energy ultrafast lasers and the gas
density limits imposed by conventional time-of-flight experiment designs [5].

Yamakawa et al. and Chowdhury et al. performed precision measurements of
highly-charged noble gas ions at relativistic laser intensities exceeding
$2\times 10^{19}$ W/cm2[6][7]. Both authors exploited the strongly nonlinear
dependence of tunneling ionization probability on laser intensity to calculate
a model intensity from the relative charge state yields. The intensity
computed by Chowdhury et al. was at the estimated lower bound of their
experimental intensity. Yamakawa et al. found unexpectedly that the laser
intensity calculated from the ADK-PPT tunneling model depended on the noble
gas ion species when the laser parameters were held constant, with the
intensity calculated from the model decreasing systematically with increasing
atomic number. Recent modeling suggests the ionization of helium-like ions is
more robust to uncertainties in the ionization modeling [8][9], motivating
investigation into ionization of the neon K-shell.

In this Letter we explore the observation of K-shell neon charge states
produced in a laser focus by detecting the high-energy ATI electrons produced
by the laser-ion interaction. We select high-energy ATI electrons as our
experimental observable because direct laser acceleration of the highly-
charged ions will severely degrade the resolution of time-of-flight
spectrometers as intensity approaches $10^{21}$ W/cm2[10]. The behavior of
relativistic ATI electrons is well-characterized by experiment [11][12][13],
with higher-energy electrons confined to a smaller forward cone in the laser
forward direction at an angle

$tan(\theta)=\sqrt{\frac{2}{\gamma-1}}$ (1)

from the laser propagation direction[13]. Modulation of the ATI electron
energy spectrum and spatial distributions induced by the gaps in appearance
intensity between different atomic shells has also been confirmed
experimentally in argon and xenon [12]. The large gap in ionization potential
between the L-shell ($<$ 239 eV) and K-shell electrons (1196 eV and 1362 eV)
of neon result in the K-shell electrons being “born” into a field nearly two
orders of magnitude higher in intensity, and the K-shell ATI electrons that
interact with the laser field at peak strength will be ejected in a narrower
cone centered on the laser forward direction than the L-shell electrons that
will be ponderomotively expelled by the leading edge of the laser pulse and
will therefore attain lower energies.

Figure 1 shows a simplified diagram of our experimental setup, which is
described in greater detail elsewhere [14]. A multi-Joule, linearly polarized,
ultrafast laser pulse produced by the rod amplifier of the Texas Petawatt
laser is focused with an $f$/1.4 off-axis parabolic mirror, focusing to a spot
with central maximum full-width at half-maximum measured to be $2.6\pm 0.2$
$\mu$m with optimal wavefront. The intensity at the focal plane in the target
chamber is estimated using indirect measurements of the wavefront, an
estimated pulse duration deconvolved from a second-order autocorrelation
measurement assuming a Gaussian pulse shape, and a calibrated energy
measurement. A peak intensity of $5\times 10^{20}$ W/cm2 is attained in this
configuration. Intensity was scanned by decreasing the laser energy with the
laser wavefront and pulse duration remaining optimized.

Figure 1:  A conceptual diagram of the experimental setup. Two detectors
discussed in this Letter are oriented at angles of $0^{\circ}$ and
$30^{\circ}$ from the laser forward direction and in the laser polarization
plane. The dashed lobes demonstrate angular separation of the L-shell and
K-shell electrons. This drawing is not to scale. Color figures available
online.

The focused laser pulse interacts with a low-density plume of neon gas near
the target chamber center that is introduced by a flow-calibrated orifice with
a diameter of 100 $\mu$m backed with 60 torr of ultra high-purity neon gas. We
estimate the maximum density of the gas in the interaction volume to be about
$3\times 10^{14}$ cm-3 using Ansys Fluent[15] simulation of the steady-state
gas flow into vacuum, below the threshold density for collective plasma
effects. The electrons produced by the laser-ion interactions are detected by
scintillating calorimeter detectors placed around the laser focus. The three
detectors discussed in this Letter were oriented along the polarization plane
at 30∘ from the laser forward direction, along the laser forward direction,
and at a control position 110∘ from the laser forward direction and out of the
polarization plane.

(a)

(b)

Figure 2:  Electron energy deposited in the scintillator at the
(a)$30^{\circ}$ position and (b) $0^{\circ}$ position (b). Shield thicknesses
of 1 mm and 2.6 mm have $50\%$ efficiency cutoffs at 1.4 MeV and 2.8 MeV,
respectively. Predictions of the ADK-PPT model (solid curve) and Augst BSI
model (dashed curve) are presented alongside, with simulation intensities at
the marker points. Color figures available online. The floors mark the maximum
values and intensity ranges for measurements falling below the lowest detector
charge threshold in each shielding configuration.

Each detector consisted of a 50 mm diameter, 40 mm long cylinder of long-
lifetime (285 ns) scintillating plastic (Eljen Technologies EJ240) coupled to
a photomultiplier tube with a tapered voltage divider for optimal pulse
linearity. The scintillator plastic and photomultipier tubes (PMT) were
encased in a vacuum-compatible PTFE housing that was made light-tight with
colloidal graphite and aluminum foil. The relatively large solid angle ($\sim$
0.03 steradians) subtended by the detectors captured several hundred ATI
electrons at each detector, enabling accurate calorimeter energy measurements
with only a few shots at each laser intensity. The output current pulse from
each photomultiplier tube was recorded on a Tektronix TDS5054 oscilloscope and
digitally filtered to eliminate ringing from on-shot electromagnetic noise.
The upper and lower charge thresholds at each voltage were estimated from the
breakdown of the linear relationship between current amplitude and integrated
charge due to pulse saturation and residual ringing, respectively.

We performed a series of control experiments to confirm our signal originated
from K-shell ATI electrons. We verified the signal was not due to
electromagnetic pulse effects by verifying it disappeared when the gas flow
was turned off. We verified the signal from forward-directed radiation was at
least 500 times greater than backward-detected radiation at the control
detector position. We swapped the detectors located at the $30^{\circ}$ and
control positions to confirm the observed signal was due to radiation detected
in the laser forward direction and not a detector-specific artifact. A series
of control shots was also taken using helium gas as a target to simulate the
laser interaction with the L-shell electrons. We observed no repeatable signal
that would correspond to L-shell electrons with energy greater than 1.4 MeV at
the 0∘ position and 2.8 MeV at the 30∘ position.

Figures 2a and 2b show our integrated ATI electron energy yields at the
$30^{\circ}$ and $0^{\circ}$ detector positions. Aluminum filters with
thicknesses of 1 mm and 2.6 mm are inserted in front of the detectors to stop
electrons with energy below 1.4 MeV and 2.8 MeV, respectively, to eliminate
low-energy L-shell electrons observed scattered toward the laser forward
direction [14]. We observe a similar intensity threshold effect at both
positions, where the energy carried forward by ATI electrons becomes
measurable above $2\times 10^{20}$ W/cm2 and increases rapidly with intensity.
The minimum signal level, below which detector ringing renders our estimated
uncertainty unreliable, is marked for both shielding configurations presented
in Figure 2b, showing that the ATI electron signal drops off more than an
order of magnitude at the K-shell ionization threshold intensity. A scaling
transition characteristic of tunneling ionization [4][16] is visible at an
intensity around $3\times 10^{20}$ W/cm2 in Figure 2a, where the measured ATI
electron energy yield transitions to a power-law intensity dependence
dominated by the volume of the focal region exceeding the ionization threshold
intensity. Below this intensity, the integrated ATI electron energy dependence
is dominated by the probability of K-shell ionization in the most intense
region of the laser focus, which is a highly nonlinear function of intensity.

We observe the K-shell ionization threshold intensity is nearly double the
barrier suppression intensity, and that both the ADK-PPT and Augst barrier
suppression models described in detail elsewhere [14] underestimate the
K-shell ionization intensity threshold. The primary sources of error originate
from the detector energy calibration ($\sim 20\%$) and our method of
calculating intensity ($\sim 15-20\%$), and are not sufficient to explain the
discrepancy between the observed ATI electron energy yields and the Monte
Carlo models.

(a)

(b)

Figure 3:  a) Detector efficiencies at different aluminum thicknesses (right
axis) and a simulated K-shell ATI electron spectrum using the ADK-PPT model at
intensity of $1.06\times 10^{20}$ W/cm2 (left axis) shown for comparison. b)
Integrated ATI electron energy at the two positions at two different average
intensities, with ADK-PPT simulation predictions (open markers) for
comparison.

The limited number of laser shots and the low density of target gas prevented
measurement of the electron spectrum using a magnetic spectrometer. We placed
a series of aluminum filters of different thicknesses in front of the
scintillating plastic and took repeated measurements to gain information about
the energy spectrum and the maximum electron energy. Figure 3a shows a series
of detector efficiency curves calculated using G4beamline [17] for different
aluminum filter thicknesses, with thicker filters shielding the scintillating
plastic from lower-energy electrons. Although the lack of a sharp cutoff in
the energy efficiency curves makes it impossible to invert our integrated
measurements to obtain a unique electron energy spectrum, we can gain
qualitative information on the shape of the spectrum and estimate a maximum
ATI electron energy range by comparing the ratio of the measured integrated
electron energy for the two thickest shields and to the ratio of the
respective efficiency curves.

Figure 4:  Comparison between the peak ATI electron energies predicted by
different models and our experimental results. Curves are analytic models (see
Refs. [18] and [19]) and the markers represent the average of the top 10$\%$
most energetic electrons in the ADK-PPT Monte-Carlo simulations. The shaded
red regions represent peak ATI electron consistent with the data in Fig. 3b

Figure 3b shows the measured energy yields at the $30^{\circ}$ and on-axis
detectors, where the modeling predicts the highest number of K-shell ATI
electrons. At an average intensity of $4.1\pm 0.4\times 10^{20}$ W/cm2, we
found the most favorable comparison was to the ADK-PPT ATI electron model with
a model intensity of $1.06\times 10^{20}$ W/cm2. The Augst BSI model predicted
a significantly higher ratio of energy yield in the laser forward direction
than the ADK-PPT model at every intensity, and a model intensity range
consistent with the energy yields observed at both detector positions could
not be found for the Augst BSI model [14]. From the experimental measurements
at $30^{\circ}$ at two intensities ($(4.1\pm 0.4)\times 10^{20}$ W/cm2 and
$(2.2\pm 0.4)\times 10^{20}$ W/cm2), the ADK-PPT model falls off quicker with
increasing shield thickness than the measured yields because the model
significantly underestimates the proportion of electrons with energy $>$ 6 MeV
at both model intensities. We cannot make a conclusive observation about the
K-shell electrons with energy $<3$ MeV because our helium control shots
indicate that lower-energy L-shell electrons are scattered further into the
laser forward direction than predicted by modeling [14].

The maximum ATI electron energy ranges consistent with our measured integrated
energies at these two intensities are 5.6-13 MeV and 10-16 MeV, respectively.
Figure 4 compares these ranges to different analytic models of peak ATI
electron energy: the ponderomotive, relativistic ponderomotive, and
superponderomotive “wave-particle resonance” [18][19] models. The
experimentally determined ranges fall between the relativistic and
nonrelativistic ponderomotive models, with the Monte Carlo model
overestimating the maximum ATI electron energy by a factor or 2-3. The Monte-
Carlo modeling likely overestimates the integrated ATI electron energies
because a Gaussian laser focus is assumed to make the model computationally
tractable. The predicted electron energy at the $0^{\circ}$ position will be
overestimated because higher-order spatial modes will dephase an electron co-
propagating with the laser field more quickly due to stronger Gouy phase shift
associated with higher-order modes, decreasing the maximum electron energy
along the laser direction[19]. Our experimental finding that the maximum ATI
electron energy falls between the ponderomotive and relativistic ponderomotive
models raises an important theoretical question about whether the
superponderomotive scaling of the maximum ATI electron energy at the onset of
“wave-particle resonance” predicted by D. F. Gordon et al. [18] and
demonstrated by our Monte-Carlo modeling in Figure 4 would be a feature of a
non-Gaussian laser focus.

We find that neither the ADK-PPT Monte-Carlo model nor the Augst BSI model
provide a fully consistent quantitative description of integrated ATI electron
energy measurements. While a significant threshold intensity shift is
unexpected for the ionization of a simple helium-like atom, it is not
inconsistent with the other ionization experiments at relativistic
intensities. Yamakawa et al. compared direct measurements of ion charge states
to the ADK-PPT model and found that the ADK-PPT model predicted a laser
intensity that was a factor of 2-8 lower than the indirectly calculated laser
field intensity of $2.6\times 10^{19}$ W/cm2, with heavier noble gas ions
consistent with lower field intensities [7]. While laser intensity has been
claimed to exceed the barrier suppression intensity of the neon K-shell states
since at least 2006 and charge states as high as Kr24+ have been observed
[20], our work represents the first inference of neon K-shell ionization that
has been reported in the literature to our knowledge and the first detection
of ATI electrons exceeding 10 MeV, corresponding to the absorption of $10^{7}$
photons during the ionization process.

The detection of ATI electrons from K-shell states of argon ($\sim 3\times
10^{21}$ W/cm2) and krypton ($\sim 10^{23}$) W/cm2 are relatively
straightforward extensions of our work made possible by development of
repetition-rated multi-petawatt laser systems [21]. The interaction of
Petawatt laser pulses with highly-charged ions can generate ATI electrons with
energy exceeding 100 MeV[10][19], would exit the focus in a small forward cone
that could be monitored with a magnetic spectrometer or a scintillating
detector array. Extension of this basic experimental design and the
development of computational techniques to model the interaction of ultra-
relativistic electrons with realistic spatio-temporal laser fields are active
areas of research necessary to further develop a direct intensity measurement
technique using ATI electrons from highly-charged ions in an intense laser
field.

A. Y. acknowledges helpful conversations with E. Chowdhury regarding the
design of this experiment. This work was supported by the DOE, Office of
Science, Fusion Energy Sciences under Contract No. DE-SC0021125: LaserNetUS: A
Proposal to Advance North America’s First High Intensity Laser Research
Network, the Air Force Office of Scientific Research through Awards No.
FA9550-14-1-0045 and No. FA9550-17-1-0264, the National Nuclear Security
Agency (NNSA) through Award No. NA0002008, and was prepared by LLNL under
Contract DE-AC52-07NA27344. A. Y. gratefully acknowledges the generous support
of the Jane and Michael Downer Fellowship in Laser Physics in Memory of Glenn
Bryan Focht.

## References

  * Agostini  _et al._ [1979] P. Agostini, F. Fabre, G. Mainfray, G. Petite, and N. K. Rahman, Free-Free Transitions Following Six-Photon Ionization of Xenon Atoms, Physical Review Letters 42, 1127 (1979).
  * Ammosov  _et al._ [1986] M. V. Ammosov, N. B. Delone, and V. P. Krainov, Tunnel ionization of complex atoms and of atomic ions in an alternating electromagnetic field, Journal of Experimental and Theoretical Physics 64, 1191 (1986).
  * Perelomov  _et al._ [1966] A. Perelomov, V. Popov, and M. Terent’ev, Ionization of Atoms in an Alternating Electric Field, Soviet Journal of Experimental and Theoretical Physics 23, 924 (1966).
  * Augst  _et al._ [1989] S. Augst, D. Strickland, D. D. Meyerhofer, S. L. Chin, and J. H. Eberly, Tunneling ionization of noble gases in a high-intensity laser field, Physical Review Letters 63, 2212 (1989).
  * Link  _et al._ [2006] A. Link, E. A. Chowdhury, J. T. Morrison, V. M. Ovchinnikov, D. Offermann, L. Van Woerkom, R. R. Freeman, J. Pasley, E. Shipton, F. Beg, P. Rambo, J. Schwarz, M. Geissel, A. Edens, and J. L. Porter, Development of an in situ peak intensity measurement method for ultraintense single shot laser-plasma experiments at the Sandia Z petawatt facility, Review of Scientific Instruments 77, 10E723 (2006).
  * Chowdhury  _et al._ [2001] E. A. Chowdhury, C. P. J. Barty, and B. C. Walker, “Nonrelativistic” ionization of the L -shell states in argon by a “relativistic” $10^{19}$ W/cm2 laser field, Physical Review A 63, 042712 (2001).
  * Yamakawa  _et al._ [2003] K. Yamakawa, Y. Akahane, Y. Fukuda, M. Aoyama, N. Inoue, and H. Ueda, Ionization of many-electron atoms by ultrafast laser pulses with peak intensities greater than $10^{19}$ W/cm2, Physical Review A 68, 065403 (2003).
  * Ciappina and Popruzhenko [2020] M. F. Ciappina and S. V. Popruzhenko, Diagnostics of ultra-intense laser pulses using tunneling ionization, Laser Physics Letters 17, 025301 (2020), arXiv:1911.11233 .
  * Lötstedt  _et al._ [2020] E. Lötstedt, M. F. Ciappina, and K. Yamanouchi, Static-field ionization model of He-like ions for diagnostics of light-field intensity, Physical Review A 102, 013112 (2020).
  * Yandow  _et al._ [2019] A. Yandow, T. Toncian, and T. Ditmire, Direct laser ion acceleration and above-threshold ionization at intensities from $5\times 10^{21}$ W/cm2 to $3\times 10^{23}$ W/cm2, Physical Review A 100, 053406 (2019), arXiv:arXiv:1909.02158v2 .
  * DiChiara  _et al._ [2008] A. D. DiChiara, I. Ghebregziabher, R. Sauer, J. Waesche, S. Palaniyappan, B. L. Wen, and B. C. Walker, Relativistic MeV Photoelectrons from the Single Atom Response of Argon to a $10^{19}$ W/cm2 Laser Field, Physical Review Letters 101, 173002 (2008).
  * Ekanayake  _et al._ [2013] N. Ekanayake, S. Luo, P. D. Grugan, W. B. Crosby, A. D. Camilo, C. V. McCowan, R. Scalzi, A. Tramontozzi, L. E. Howard, S. J. Wells, C. Mancuso, T. Stanev, M. F. Decamp, and B. C. Walker, Electron Shell Ionization of Atoms with Classical, Relativistic Scattering, Physical Review Letters 110, 203003 (2013).
  * Moore  _et al._ [1995] C. I. Moore, J. P. Knauer, and D. D. Meyerhofer, Observation of the Transition from Thomson to Compton Scattering in Multiphoton Interactions with Low-Energy Electrons, Physical Review Letters 74, 2439 (1995).
  * Yandow  _et al._ [2023] A. Yandow, T. N. Ha, C. Aniculaesei, H. L. Smith, C. G. Richmond, M. M. Spinks, H. J. Quevedo, S. Bruce, D. A. Garcia, M. Darilek, C. Chang, E. Gaul, M. E. Donovan, B. M. Hegelich, and T. Ditmire, Above threshold ionization at laser intensity greater than $10^{20}$ W/cm2 (2023).
  * [15] Ansys Fluent, Release 12.0.
  * Augst  _et al._ [1991] S. Augst, D. D. Meyerhofer, D. Strickland, and S. L. Chin, Laser ionization of noble gases by Coulomb-barrier suppression, Journal of the Optical Society of America B 8, 858 (1991).
  * Roberts and Kaplan [2007] T. J. Roberts and D. M. Kaplan, G4beamline simulation program for matter-dominated beamlines, in  _2007 IEEE Particle Accelerator Conference (PAC)_ (IEEE, 2007) pp. 3468–3470.
  * Gordon  _et al._ [2017] D. F. Gordon, J. P. Palastro, and B. Hafizi, Superponderomotive regime of tunneling ionization, Physical Review A 95, 033403 (2017), arXiv:1511.03915 .
  * Maltsev and Ditmire [2003] A. Maltsev and T. Ditmire, Above Threshold Ionization in Tightly Focused, Strongly Relativistic Laser Fields, Physical Review Letters 90, 053002 (2003).
  * Akahane  _et al._ [2006] Y. Akahane, J. Ma, Y. Fukuda, M. Aoyoma, H. Kiriyama, J. V. Sheldakova, A. V. Kudryashov, and K. Yamakawa, Characterization of wave-front corrected 100 TW, 10 Hz laser pulses with peak intensities greater than $10^{20}$ W/cm2, Review of Scientific Instruments 77, 023102 (2006).
  * Rus  _et al._ [2015] B. Rus, P. Bakule, D. Kramer, J. Naylon, J. Thoma, J. T. Green, R. Antipenkov, M. Fibrich, J. Novák, F. Batysta, T. Mazanec, M. A. Drouin, K. Kasl, R. Baše, D. Peceli, L. Koubíková, P. Trojek, R. Boge, J. C. Lagron, Š. Vyhlídka, J. Weiss, J. Cupal, J. Hřebíček, P. Hříbek, M. Durák, J. Polan, M. Košelja, G. Korn, M. Horáček, J. Horáček, B. Himmel, T. Havlíček, A. Honsa, P. Korouš, M. Laub, C. Haefner, A. Bayramian, T. Spinka, C. Marshall, G. Johnson, S. Telford, J. Horner, B. Deri, T. Metzger, M. Schultze, P. Mason, K. Ertel, A. Lintern, J. Greenhalgh, C. Edwards, C. Hernandez-Gomez, J. Collier, T. Ditmire, E. Gaul, M. Martinez, C. Frederickson, D. Hammond, C. Malato, W. White, and J. Houžvička, ELI-Beamlines: development of next generation short-pulse laser systems, Research Using Extreme Light: Entering New Frontiers with Petawatt-Class Lasers II 9515, 95150F (2015).


# Wavelet–based tools to analyze, filter, and reconstruct transient
gravitational–wave signals.

Andrea Virtuoso<EMAIL_ADDRESS>Università di Trieste,
Dipartimento di Fisica, I-34127 Trieste, Italy;
INFN Sezione di Trieste, I-34127 Trieste, Italy  Edoardo Milotti  Università
di Trieste, Dipartimento di Fisica, I-34127 Trieste, Italy;
INFN Sezione di Trieste, I-34127 Trieste, Italy

###### Abstract

The analysis of gravitational–wave signals is one of the most challenging
application areas of signal processing, because of the extreme weakness of
these signals and of the great complexity of gravitational–wave detectors.
Wavelet transforms are specially helpful in detecting and analyzing
gravitational–wave transients and several analysis pipelines are based on
these transforms, both continuous and discrete. While discrete wavelet
transforms have distinct advantages in terms of computing efficiency,
continuous wavelet transforms (CWT) produce smooth and visually stunning
time–frequency maps where the wavelet energy is displayed in terms of time and
frequency. In addition to wavelets, short–time Fourier transforms (STFT) and
Stockwell transforms (ST) are also used, or the Q–transform, which is a Morlet
wavelet–like transform where the width of the Gaussian envelope is
parameterized by a parameter denoted by $Q$ [Chatterji et al., Class. Quant.
Gravity 21 (2004) S1809]. To date, the use of CWTs in gravitational–wave data
analysis has been limited by the higher computational load when compared with
discrete wavelets, and also by the lack of an inversion formula for wavelet
families that do not satisfy the admissibility condition. In this paper we
consider Morlet wavelets parameterized in the same way as the Q–transform
(hence the name wavelet Q–transform) which have all the advantages of the
Morlet wavelets and where the wavelet transform can be inverted with a
computationally efficient specialization of the non–standard inversion formula
of Lebedeva and Postnikov [Lebedeva and Postnikov, Royal Society Open Science,
1 (2014) 140124]. We also introduce a two–parameter extension (the wavelet
Qp–transform) which is well–adapted to chirping signals like those originating
from compact binary coalescences (CBC), and show that it is also invertible
just like the wavelet Q–transform. The inversion formulas of both transforms
allow for effective noise filtering and produce very clean reconstructions of
gravitational–wave signals. Our preliminary results indicate that the method
could be well suited to perform accurate tests of General Relativity by
comparing modeled and unmodeled reconstructions of CBC gravitational–wave
signals.

wavelets, Q–transform, inverse transform, denoising, Qp–transform,
gravitational waves

##  I Introduction

Since the first detection of gravitational waves (GW) in 2015 [5, 6, 4] the
LIGO–Virgo–KAGRA (LVK) collaboration has successfully analyzed and
reconstructed gravitational–wave transients (GWT) with different methods [1,
10, 12] which fall into two large classes, the modeled methods and the
unmodeled methods. Modeled methods use theoretical waveforms and
matched–filtering to detect and analyze GWTs [55, 53, 24, 17, 52] and they
have been hugely successful in the inference of the properties of merging
compact binary objects (see, e.g., [4], or [12] and references therein).
Unmodeled methods are based instead on the coherence of data in different
detectors, and make no assumptions on signal waveform [7, 2, 14, 30, 26, 37,
23, 48, 50]. Even though they are not as sensitive as modeled methods there
are two main reasons to continue developing and using them: first, their
non–reliance on prior models provides independent checks in the detection and
analysis of GWTs [6]. Second, unmodeled methods are effective also when
dealing with GW sources where we lack precise theoretical waveforms, such as
core–collapse supernovae (CCSN) [3, 49, 51].

The lack of well–defined reference waveforms requires the adoption of complex
mathematical tools to extract signals from background noise. The
representations of signals in the time–frequency (TF) domain are specially
useful to disentangle signals from noise and efficiently evaluate the
coherence of data in different detectors. Such representations are obtained
with a variety of methods, like short–time Fourier transforms (STFT),
Stockwell transforms (ST) [47], wavelet transforms (WT) [34], and several
variants that mix the different approaches.

While all transforms are superficially similar, they differ in important
details, such as the existence of an inverse transform. The invertibility of a
TF representation is crucial for further processing of the signals, such as
filtering out unwanted noise. Continuous wavelet transforms (CWT) do have an
exact reconstruction formula if and only if the mother wavelet $\psi(x)$
satisfies the admissibility condition

$C_{\psi}=\int_{-\infty}^{+\infty}\frac{|\tilde{\psi}(y)|^{2}}{|y|}dy<\infty$
(1)

where $\tilde{\psi}(y)$ is the Fourier transform of $\psi(x)$. This is
substantially equivalent to the condition that the mother wavelet has a
vanishing integral

$\int_{-\infty}^{+\infty}\psi(x)dx=0$ (2)

(see, e.g., [25], Paragraph 2.4). When the admissibility condition is
satisfied, the wavelet transform of the time–domain signal $x(t)$

$T_{x}(\nu,\tau)=\int_{-\infty}^{+\infty}x(t)\psi^{*}_{\nu,\tau}(t)dt$ (3)

is invertible with the standard formula

$x(t)=\frac{1}{C_{\psi}}\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}d\nu\;d\tau\;T_{x}(\nu,\tau)\psi_{\nu,\tau}(t)$
(4)

Interestingly, other inversion formulas exist (see, e.g., [25], Section 2.4)
which also require the admissibility condition.

A TF representation should also have an excellent resolution, both in time and
in frequency. However, time resolution $\sigma_{\tau}$ and frequency
resolution $\sigma_{\nu}$ can not be both arbitrarily small, they are limited
by the Gabor–Heisenberg principle [27]

$\sigma_{\tau}\sigma_{\nu}\geq\dfrac{1}{4\pi}$ (5)

The lower bound of this inequality is actually attained by the Morlet wavelets
(see eqs. (9) and (10)).

Finally, the signal energy may be spread over regions of arbitrary shape in
the TF plane, and this means that in order to achieve a sparse TF
representation — i.e., a representation that concentrates signal energy in as
few components as possible – we should be able to adapt the time and frequency
resolutions to follow the main features of the TF representation of a signal.
Sparse representations are specially useful in the analysis of GW signals,
because sparseness favors signal detectability and helps signal filtering and
reconstruction.

Since the product of time resolution and frequency resolution is fixed, it
should be possible to set a balance between resolutions to optimize the
analysis of an input signal. As an example, STFT lacks this tunability because
it has a fixed envelope window (and therefore fixed time and frequency
resolutions) for all frequencies: this implies that at different frequencies
the window encompasses a different number of oscillations, which is often
undesirable [42]. On the contrary, WTs feature envelope windows with variable
width such that the number of oscillations is frequency–independent, resulting
in variable time and frequency resolutions at different central frequencies.
This scale invariance is a precious feature which we retain in our analysis of
GWTs, just as in past works [38, 30, 23], and which helps achieving sparse TF
representations of signals.

The issue of setting an optimal balance of time and frequency resolution is
addressed by the Q–transform with a Gaussian window function [18, 21, 22, 43,
28], which is defined by the pair of equations

$\displaystyle X(\tau,\nu,Q)$
$\displaystyle=\int_{-\infty}^{+\infty}dt\,s(t)g(t;\tau,\nu,Q)e^{-2\pi i\nu
t}$ (6a) $\displaystyle g(t;\tau,\nu,Q)$
$\displaystyle=\left(\frac{8\pi\nu^{2}}{Q^{2}}\right)^{1/4}\,e^{-\left(\frac{2\pi\nu(t-\tau)}{Q}\right)^{2}}.$
(6b)

where $s(t)$ is the time–domain representation of the signal for which the
Q–transform is evaluated, $\tau$ is the central time, $\nu$ the central
frequency and $Q$ a parameter that sets the balance between time
($\sigma_{\tau}$) and frequency ($\sigma_{\nu}$) resolution. The Q–transform
is not a continuous wavelet transform: even if it inherits from wavelet
transforms the property that its window $g(t;\tau,\nu,Q)$ depends on the
frequency $\nu$ — a key property of wavelets — the oscillatory term $e^{-2\pi
i\nu t}$ is not centered at $t=\tau$, as for wavelets, but is the same for all
times.

Still, the Q–transform is very similar to the Morlet wavelet transform, and
here we consider a version of the Morlet wavelets with the addition of the $Q$
parameter as in the Q–transform for tunability of the time and frequency
resolutions. The wavelet version of the Q–transform that we consider in this
paper drops the requirement of a common oscillatory term in the Q–transform,
which is replaced by $e^{-2\pi i\nu(t-\tau)}$. With this substitution, we
obtain a wavelet version of the Q–transform, the _wavelet Q–transform_

$T(\tau,\nu,Q)=\int_{-\infty}^{+\infty}dt\,s(t)\psi^{*}(t;\tau,\nu)$ (7)

with the wavelets

$\psi^{*}(t;\tau,\nu,Q)=\left(\frac{8\pi\nu^{2}}{Q^{2}}\right)^{1/4}\,e^{-\left(\frac{2\pi\nu(t-\tau)}{Q}\right)^{2}-2\pi
i\nu(t-\tau)}.$ (8)

where ∗ denotes complex conjugate. The time and frequency resolutions depend
only on the Gaussian window, and are found to be the same for both Q–transform
and wavelet Q–transform

$\displaystyle\sigma_{\tau}^{\mathrm{(Q)}}$ $\displaystyle=\dfrac{Q}{4\pi\nu}$
(9) $\displaystyle\sigma_{\nu}^{\mathrm{(Q)}}$ $\displaystyle=\dfrac{\nu}{Q}.$
(10)

It is important to note that both resolutions depend on $\nu$. Moreover, while
low values of $Q$ ($\sim 2\pi$) favor a better time resolution with respect to
frequency resolution, high values of $Q$ ($\gg 2\pi$) give good frequency
resolution and poor time resolution. In summary, when computing the wavelet
Q–transform with high values of $Q$ the signal energy in the TF plane is more
spread in time compared to low values of $Q$. In contrast, with low values of
$Q$, the signal energy in the TF plane has a higher spread in frequency
compared to high values of $Q$. In all cases, the product of the resolutions
remains fixed at the lower bound of Gabor–Heisenberg uncertainty inequality
and endows this transform with the best possible TF resolution111In actual
implementations the tails of the Gaussian window are truncated for two
reasons: first, to improve computational efficiency, and second, because with
real signals which are time– and band–limited at some point this truncation
has to be done anyway. As a result, the truncated wavelets have a slightly
reduced resolution.

Unfortunately, the Morlet wavelets do not match the admissibility condition,
and the resulting WT is not invertible with the standard formula: since the
wavelet Q–transform is just a parameterized version of the Morlet wavelet
transform it is not invertible as well.222The BayesWave pipeline [37, 23],
which is based on Morlet wavelets, bypasses the non–invertibility issue by
picking one by one individual wavelets within a Bayesian framework which
utilizes a reversible–jump Markov chain Monte Carlo (RJMCMC), and then finally
superposing all of them to reconstruct the signal. This method produces good
reconstructions but the RJMCMC step is computationally very expensive. The
non–invertibility of Morlet wavelets can also be overcome by shifting them to
satisfy the admissibility condition [16]: in GW data analysis, that
possibility has been exploited in [44]. However, recent progress on a wavelet
inversion formula that does not require the admissibility condition [31, 40]
means that the wavelet Q–transform can provide both TF representations
comparable to those of the standard Q–transform and all the advantages of
invertibility. We make further progress in this direction and here we present
a more efficient inversion formula which also implements TF filtering. This
new inversion formula, which we discuss in Section II, is a real game–changer,
because it leads to computationally efficient denoising and waveform
reconstruction with the wavelet Q–transform.

Next, in Section III we describe a variant of the wavelet Q–transform that is
better adapted at describing GW chirps from CBCs and is still graced with an
inverse formula. With this parameterization, GW chirps are better modeled and
the Gaussian window can be made wider, as we discuss in Section III.3, giving
in many cases a better match of the wavelet with the chirping signal, as shown
in Section IV. Because of the new $p$ parameter, we call this new transform
the wavelet Qp–transform.

In Section IV we apply the wavelet Qp–transform and the associated denoising
method discussed in Section II to some well known GW events (GW150914,
GW170817 and GW190521 [5, 4, 6, 8, 11]), and compare the performance of the
wavelet Qp–transform to that of wavelet Q–transform. We also evaluate the
waveforms obtained with the denoising formula for both transforms, comparing
these with a modeled reconstruction based on matched filtering. Finally, we
conclude the paper with a short discussion in Section V.

##  II Non-standard reconstruction formula for the wavelet Q–transform

As explained in the Introduction, a major drawback of the Morlet wavelets, and
therefore of the wavelet Q–transform, is that they do not satisfy the
admissibility condition. Here we introduce a modified form of the alternative
wavelet reconstruction formula recently proposed by Lebedeva and Postnikov
[31, 40] which does not require the admissibility condition. A complete
derivation of the reconstruction formula for the wavelet Q–transform is given
in the Supplementary Text, and the result is

$s(\tau)=\dfrac{2}{\left(\dfrac{2}{\pi
Q^{2}}\right)^{1/4}\text{Re}\left[\text{erf}\left(\frac{Q}{2}\right)\right]}\,\text{Re}\left[\int_{0}^{+\infty}d\nu\,\dfrac{1}{i\sqrt{2\pi\nu}2\pi\nu}\dfrac{\partial}{\partial\tau}T(\tau,\nu,Q)\right].$
(11)

This new reconstruction formula is our first important result, which
represents a considerable advantage with respect to the traditional
Q–transform [21, 22, 43]. After discretization of the integral in eq. (11) as
usually done with the standard formula (4) [29], the new formula allows the
application of specific TF filters.

In general, we can perform filtering by multiplying the transform
$T(\tau,\nu)$ times a filter function. When filtering means satisfying a
predefined criterion C, e.g., passing an amplitude threshold, then the filter
is represented by an indicator function $\mathbf{1}_{C}(\tau,\nu)$ equal to 1
wherever the criterion is satisfied and 0 elsewhere. We find that an important
subclass of time–frequency filters, those that select only the regions bounded
by the $L$ time–dependent frequency intervals333We take into account the
possibility of having multiple unconnected intervals. This is important, e.g.,
when higher order modes are present in CBC chirps [54], or in the case of
overlapping signals [41].
$\left[\nu^{\mathrm{low}}_{l}(\tau),\nu^{\mathrm{high}}_{l}(\tau)\right]$ with
$l=1,\,...\,,L$, can be implemented in a computationally efficient way with
the following integral expression

$s_{D}(\tau)=\text{Re}\left[\int_{-\infty}^{+\infty}df\,\tilde{s}(f)e^{2\pi
if\tau}w(\tau,f,Q)\right],$ (12)

where $s_{D}(\tau)$ is the denoised signal and where we have introduced the TF
window

$w(\tau,f,Q)=\dfrac{1}{\text{erf}\left(\frac{Q}{2}\right)}\sum_{l}\left\\{\text{erf}\left[\dfrac{Q}{2}\left(\dfrac{f-\nu^{\mathrm{low}}_{l}(\tau)}{\nu^{\mathrm{low}}_{l}(\tau)}\right)\right]-\text{erf}\left[\dfrac{Q}{2}\left(\dfrac{f-\nu^{\mathrm{high}}_{l}(\tau)}{\nu^{\mathrm{high}}_{l}(\tau)}\right)\right]\right\\}.$
(13)

As an example, consider Figure 1 which shows the real part of a pair of
Q-wavelets, and Figure 2 which shows the corresponding window functions for a
single pair of values $\nu^{\mathrm{low}}$ and $\nu^{\mathrm{high}}$.

Figure 1: Real part of a pair of Q–wavelets. The wavelet on the left has been
generated with $Q=50$, $\tau=-1$ s, and $\nu=15$ Hz, the one on the right with
$Q=10$, $\tau=1$ s, and $\nu=15$ Hz. The standard deviation of the Gaussian
envelope is $\approx 0.34$ s for the wavelet on the left, and $\approx 0.07$ s
for the wavelet on the right. Figure 2: Window functions for the pair of
wavelets in Figure 1 and a single pair of values $\nu^{\mathrm{low}}=10$ Hz
and $\nu^{\mathrm{high}}=20$ Hz. The solid line corresponds to $Q=50$ and the
dashed line to $Q=10$. The wavelet with the longer duration produces a window
with sharper edges in the frequency domain.

The window functions are equivalent to narrow bandpass filters in the
frequency domain, with the larger $Q$’s producing a sharper filter functions.
Neither filter has a compact support, however it is clear that the decay of
the filter tails is very fast, and therefore both filters are very well
localized in the sense first explained by Slepian [45, 46]. Notice also that
both filter functions are real and therefore act as non–causal filters.

Just as the coherence among detectors, a careful selection of the frequency
intervals in the TF plane can also be used to reject background noise; in the
analyses that we report in this paper we retain only those TF intervals having
an energy above a fixed threshold, as described in Section IV.

##  III The wavelet Qp–transform

###  III.1 A wavelet transform with chirping–frequency wavelets

In this Section we introduce a variant of the wavelet Q–transform where we let
the frequency of the wavelet change to adapt to the shape of gravitational
wave chirps. A similar idea was developed previously, in the context of the
so–called “chirplets” [37, 35, 36, 19, 20, 39, 58], with the goal of following
the TF evolution of the CBC inspiral phase. This was achieved either by
introducing the first derivative of the wavelet frequency as a parameter

$\nu(t)=\nu(\tau)+\mathinner{\dfrac{\partial\nu}{\partial
t}\Biggr{\rvert}}_{t=\tau}(t-\tau).$ (14)

and by deploying a wide array of different methods to place chirplets in the
TF plane [37, 35, 36, 19, 20, 39, 58], or by modeling the TF evolution of the
signal with a more complex polynomial parameterization [57, 32]. The wavelet
Qp–transform is based on the same idea, but with an important difference: the
derivative ${\partial\nu}/{\partial t}|_{t=\tau}$ is itself a function of the
local wavelet frequency over the wavelet duration
$[\tau-\sigma_{\tau},\tau+\sigma_{\tau}]$:

$\displaystyle\mathinner{\dfrac{\partial\nu}{\partial
t}\Biggr{\rvert}}_{t=\tau}$
$\displaystyle=\dfrac{(1+p)\nu(\tau)-(1-p)\nu(\tau)}{(\tau+\sigma_{\tau})-(\tau-\sigma_{\tau})}=4\pi\nu^{2}(\tau)\dfrac{p}{Q}$
$\displaystyle=\dfrac{1}{2\pi}\left(\dfrac{2\pi\nu(\tau)}{Q}\right)^{2}2pQ,$
(15)

where we have used expression (9) for $\sigma_{\tau}$. The $p$ parameter
determines the fractional frequency change with respect to the central
frequency over the $[\tau-\sigma_{\tau},\tau+\sigma_{\tau}]$ time interval.
The resulting nonlinear frequency change helps these wavelets in adapting to
the CBC chirping gravitational–wave signals.

As a result, the (conjugate) Qp–wavelets are described by the following
expression:

$\psi^{*}(t;\tau,\nu,Q,p)=\left(\frac{8\pi\nu^{2}}{Q^{2}}\right)^{1/4}\,e^{-\left(2\pi\nu(t-\tau)\frac{\sqrt{1+2ipQ}}{Q}\right)^{2}-2\pi
i\nu(t-\tau)}.$ (16)

Figure 3 shows a pair of Qp–wavelets with different central frequencies and
positive $p$, which determines a chirp with frequency increasing in the
positive directions, just as the CBC gravitational–wave chirps.

Figure 3: Real part of a pair of Qp–wavelets. The wavelet on the left has been
generated with $Q=50$, $\tau=-1$ s, $\nu=15$ Hz, and $p=0.1$, the one on the
right with $Q=10$, $\tau=1$ s, $\nu=15$ Hz, and $p=0.1$. Because of the
positive value of $p$ the frequency increases in the positive direction of
time.

It is important to note that for Qp–wavelets the time resolution is the same
as in eq. (9), while the frequency resolution differs from that found in eq.
(10). Indeed, it can be shown that in the case of the Qp wavelets

$\displaystyle\sigma_{\tau}^{\mathrm{(Qp)}}$
$\displaystyle=\dfrac{Q}{4\pi\nu}$ (17)
$\displaystyle\sigma_{\nu}^{\mathrm{(Qp)}}$
$\displaystyle=\dfrac{\nu}{Q}\sqrt{1+(2pQ)^{2}}.$ (18)

so that the lower bound of the Gabor–Heisenberg uncertainty product is no
longer reached if $p\neq 0$. This reduced resolution comes from the frequency
spread due to the chirping behavior and is discussed further in Section III.3.
Using the Qp wavelets in eq. (16) we can write the wavelet Qp–transform of a
signal $s(t)$ as

$\displaystyle T(\tau,\nu,Q,p)$
$\displaystyle=\int_{-\infty}^{+\infty}dt\,s(t)\psi^{*}(t;\tau,\nu,Q,p)$ (19a)
$\displaystyle=\int_{-\infty}^{+\infty}dt\,s(t)\left(\frac{8\pi\nu^{2}}{Q^{2}}\right)^{1/4}\,e^{-\left(2\pi\nu(t-\tau)\frac{\sqrt{1+2ipQ}}{Q}\right)^{2}}e^{-2\pi
i\nu(t-\tau)}$ (19b)
$\displaystyle=\int_{-\infty}^{+\infty}df\,\tilde{s}(f)\tilde{\psi}^{*}(f;\tau,\nu,Q,p)$
(19c)
$\displaystyle=\int_{-\infty}^{+\infty}df\,\tilde{s}(f)\left(\frac{1}{2\pi\nu^{2}Q^{2}}\right)^{1/4}\frac{Q}{\sqrt{1+2iQp}}e^{-\left(\frac{Q}{2\sqrt{1+2iQp}}\frac{f-\nu}{\nu}\right)^{2}}e^{2\pi
if\tau}$ (19d)

where Fourier transforms are denoted by the overhead tilde. When computing the
wavelet Qp–transform, the $p$ parameter aligns the wavelets along the chirp:
positive values of $p$ are associated with a counterclockwise rotation, which
results in a better match of the transform with frequency–increasing signals;
negative values of $p$ are associated with a clockwise rotation, which results
in a better match of the transform with frequency–decreasing signals.
Therefore, the reduced frequency resolution is offset by the better match of
wavelet and signal. If $p=0$, i.e., if there is no frequency chirp, we recover
all the results obtained in the case of the wavelet Q–transform.

###  III.2 Non-standard inversion formula for wavelet Qp–transform

The inversion formula for the wavelet Q–transform, eq. (11), as well as its
denoising version, eqs. (12) and (13), can be extended to the wavelet
Qp–transform. The resulting formulas, which are proved in the Supplementary
Text, are:

$s(\tau)=\dfrac{2}{\left(\dfrac{2}{\pi
Q^{2}}\right)^{1/4}\text{Re}\left[\text{erf}\left(\frac{Q}{2\sqrt{1+2iQp}}\right)\right]}\,\text{Re}\left[\int_{0}^{+\infty}d\nu\,\dfrac{1}{i\sqrt{2\pi\nu}2\pi\nu}\dfrac{\partial}{\partial\tau}T(\tau,\nu,Q,p)\right].$
(20)
$s_{D}(\tau)=\text{Re}\left[\int_{-\infty}^{+\infty}df\,\tilde{s}(f)e^{2\pi
if\tau}w(\tau,f,Q,p)\right]$ (21)
$w(\tau,f,Q,p)=\dfrac{1}{\text{Re}\left[\text{erf}\left(\frac{Q}{2\sqrt{1+2iQp}}\right)\right]}\sum_{l}\left[\text{erf}\left(\dfrac{Q}{2\sqrt{1+2iQp}}\left(\dfrac{f-\nu^{\mathrm{low}}_{l}(\tau)}{\nu^{\mathrm{low}}_{l}(\tau)}\right)\right)-\text{erf}\left(\dfrac{Q}{2\sqrt{1+2iQp}}\left(\dfrac{f-\nu^{\mathrm{high}}_{l}(\tau)}{\nu^{\mathrm{high}}_{l}(\tau)}\right)\right)\right].$
(22)

Again, notice that setting $p=0$ we recover the previous reconstruction
formulas for the wavelet Q–transform. It is important to remark that they
critically depend on the parameterization of the frequency derivative, which
leads to a solvable Gaussian integral.

Figure 4 shows an example of the complex window $w(\tau,f,Q,p)$ in the case of
the wavelet Qp–transform, which is a complex function. We notice that the
window amplitude in the frequency domain (middle panel) is slightly irregular:
this is a manifestation of the Gibbs phenomenon for this particular kind of
wavelets. It is also important to note that the window phase is
frequency–dependent (bottom panel): this window acts again as a filter, but,
unlike the window for the wavelet Q–transform, it does not behave as a
non–causal filter as a consequence of the asymmetrical shape of the wavelet in
the time domain. Still, the phase is close to zero between
$\nu^{\mathrm{low}}$ and $\nu^{\mathrm{high}}$, and this means that for
frequencies in this range the window behavior is close to that of a non–causal
filter. Numerical experiments indicate that, as expected, the non–causal
behavior is more prominent for larger values of $p$ and for higher ratios
$(\nu^{\mathrm{high}}-\nu^{\mathrm{low}})/Q$. This means that by keeping these
values low, in specific applications we can safely use the phase information
provided by the wavelet Qp–transform to complement the amplitude information.

Figure 4: Window function determined by Qp–wavelets with different $Q$’s. Here
we take a single pair of values $\nu^{\mathrm{low}}=10$ Hz and
$\nu^{\mathrm{high}}=20$ Hz (marked by the red vertical lines), $\nu=15$ Hz,
and $p=0.05$. In each panel, the black solid line corresponds to $Q=100$, and
the black dashed line corresponds to $Q=5$. The top panel shows the real part
of the window function $w(\tau,f,Q,p)$, the middle panel shows its absolute
value, the bottom panel shows its phase (in radians). Outside the 10 Hz – 20
Hz interval, the phase of the window function oscillates wildly.

###  III.3 Practical implementation of the wavelet Q– / Qp–transform

The continuous wavelet transforms that we consider in this paper are applied
to data samples $s(t_{n})$ ($n=0,...,N-1$) taken with sampling rate $f_{s}$.
The samples are whitened using an interpolated version of the power spectral
density that defines the spectral sensitivity of the detector and which is
independent from the sample rate with which the power spectral density is
obtained. On the other hand, since the width of the frequency bins is a
function of the sampling rate each frequency bin in the transform of a sampled
signal picks up a total noise power which is a function of the sampling rate.
For this reason, a correct comparison with the noise background, and therefore
a correct whitening, requires a transform that does not depend on the sampling
rate. In the Supplementary Text we prove that the following scaled
(dimensionless) version of equation (19c) satisfies this requirement

$T_{\mathrm{nd}}(\tau,\nu,Q,p)=\dfrac{\sqrt{f_{s}}}{N}\sum_{m}\,\tilde{s}(f_{m})\tilde{\psi}^{*}(f_{m};\tau,\nu,Q,p)$
(23)

where $\tilde{s}(f_{m})=\sum_{n}s(t_{n})e^{-2\pi if_{m}t_{n}}$ and
$m=-N/2+1,...,+N/2$. As a result, the individual values of the wavelet
Qp–transform $|T(\tau,\nu,Q,p)|^{2}$ corresponding to the whitened noise
background, i.e., with Gaussian white noise $~{}N(0,1)$, follow a $\chi^{2}$
distribution with $2$ degrees of freedom, with mean value and standard
deviation equal to $2$ for all reasonable values of $Q$ and $p$ involved in GW
data analysis (see the Supplementary Text for details).

In addition to the obvious discretization of time $t$ and frequency $f$ due to
sampling, and to the numerical problems that this sampling carries with it
[29], we must also select a discrete tiling of the TF plane $(\tau,\nu)$ to
estimate the continuous wavelet transforms in a graphically satisfactory and
computationally efficient way. To this end, we address the following two
issues:

  * •

choice of a tiling represented by a discrete lattice of points
$(\tau_{i},\nu_{j})$ in the TF plane for the evaluation and representation of
the transform, assuming a time range
$[\tau_{\mathrm{min}},\tau_{\mathrm{max}}]$ and a frequency range
$[\nu_{\mathrm{min}},\nu_{\mathrm{max}}]$;

  * •

informed selection of $Q$ and $p$ to produce a signal representation as sparse
as possibile.

We can link the positions $(\tau_{i},\nu_{j})$ of the tiles to the TF
resolutions, eq. (17) and eq. (18), of the wavelets in the following way:
first, we consider a set of frequencies $\\{\nu_{j}\\}$ covering the
$[\nu_{\mathrm{min}},\nu_{\mathrm{max}}]$ range such that

$\nu_{j+1}=\nu_{j}+\alpha\sigma_{\nu_{j}}=\nu_{j}\left(1+\frac{\alpha}{Q}\sqrt{1+(2pQ)^{2}}\right)$
(24)

where the tiling parameter $\alpha$ has been introduced and where the notation
$\sigma_{\nu_{j}}$ stresses the actual dependence of the uncertainty on
$\nu_{j}$. Next, for each frequency $\nu_{j}$, the times $\tau_{i}(\nu_{j})$
are chosen accordingly

$\tau_{i+1}(\nu_{j})=\tau_{i}(\nu_{j})+\alpha\sigma_{\tau_{i}}(\nu_{j})=\tau_{i}(\nu_{j})+\alpha\dfrac{Q}{4\pi\nu_{j}},$
(25)

where, again, the notation stresses the dependence on frequency. The resulting
tiling is similar to that used for the multiresolution frames of discrete
wavelets [25]. Since the tile area is the constant
$\alpha^{2}\sqrt{1+(2pQ)^{2}}/4\pi$ and the time and frequency ranges
correspond to the TF area

$(\tau_{\mathrm{max}}-\tau_{\mathrm{min}})(\nu_{\mathrm{max}}-\nu_{\mathrm{min}})$

the total number of tiles is

$N_{\mathrm{tiles}}=\dfrac{(\tau_{\mathrm{max}}-\tau_{\mathrm{min}})(\nu_{\mathrm{max}}-\nu_{\mathrm{min}})}{\alpha^{2}}\dfrac{4\pi}{\sqrt{1+(2pQ)^{2}}}.$
(26)

Notice that a smooth representation of the continuous wavelet transform
requires a very redundant representation where $\alpha\leq 1$, and that the
computational load for the wavelet Q– and Qp–transforms scales as

$\dfrac{(\tau_{\mathrm{max}}-\tau_{\mathrm{min}})(\nu_{\mathrm{max}}-\nu_{\mathrm{min}})}{\alpha^{2}}.$

Moving on to the selection of optimal values of $Q$ and $p$, we note first
that we certainly need at least a few large–amplitude cycles in the wavelet
for a good determination of the instantaneous frequency of a signal. We can
compute the number of cycles $k$ in the range
$[\tau-n\sigma_{\tau},\tau+n\sigma_{\tau}]$ by integrating the phase
$2\pi\nu(t)$ over $t$ in the same range and dividing by $2\pi$, i.e.,

$k=\int_{\tau-n\sigma_{\tau}}^{\tau+n\sigma_{\tau}}\left[\nu(\tau)+\mathinner{\dfrac{\partial\nu}{\partial
t}\Biggr{\rvert}}_{t=\tau}(t-\tau)\right]\,dt.$ (27)

Using eq. (15), with little additional algebra we find

$Q=\dfrac{k4\pi}{2n+n^{2}p}.$ (28)

Setting $k>1$ and assuming $p\ll 1$, we find the bound

$Q>\dfrac{4\pi}{2+p}\sim 2\pi.$ (29)

Therefore, from now on we take $Q\gtrsim 2\pi$.

There are no similar bounds on $p$ although we expect a limited frequency
evolution within the same $[\tau-n\sigma_{\tau},\tau+n\sigma_{\tau}]$ range,
which leads to the rule–of–thumb $p\lesssim 1/Q$.

Since we aim at a sparse representation [34] after removing the noise
background, we filter out the tiles in the TF plane that are below a preset
energy threshold and select those values of $Q$ and $p$ that maximise the
energy density $\varepsilon$ of the surviving tiles

$\varepsilon=\dfrac{\int\int_{|T(\tau,\nu,Q,p)|^{2}>\mathrm{thr}}|T(\tau,\nu,Q,p)|^{2}d\tau
d\nu}{\int\int_{|T(\tau,\nu,Q,p)|^{2}>\mathrm{thr}}d\tau d\nu}.$ (30)

The procedure satisfies two goals: it maximizes the total energy of the
selected tiles and at the same time it minimizes their TF area, leading to a
representation of the signal as compact as possible.

In this context, it is interesting to compare the wavelet Q–transform and the
wavelet Qp–transform: since Qp wavelets better adapt to a chirp, they can have
a longer duration, i.e., a larger $Q$, and fewer above–threshold wavelets are
actually needed to represent a chirp. This gain in terms of sparsity of
wavelet representation is balanced by a loss in time resolution: a larger $Q$
corresponds to a larger TF uncertainty product.

##  IV Application to a few notable GW events

In this Section we apply the wavelet Q– and Qp–transforms to three well–known
events detected during the first three LIGO–Virgo observing runs: GW150914,
GW170817 and GW190521. The data are public and are taken from the
Gravitational–Wave Open Science Center website (GWOSC) [13, 15]. In each case,
data sampled at $4096\text{\,}\mathrm{H}\mathrm{z}$ are downloaded,
downsampled at $2048\text{\,}\mathrm{H}\mathrm{z}$, and finally whitened
before analysis 444For the decimation step we used the SciPy [56] function
scipy.signal.decimate, and for the whitening step the GWpy [33] method whiten
defined in the class TimeSeries..

###  IV.1 GW150914

The first GW event was detected by the LIGO–Virgo collaboration in 2015 [5]
and widely analyzed with both match–filtered methods and with unmodeled
methods [4, 6]. Here, we reanalyze the public data using both the wavelet Qp–
and the wavelet Q–transform on LIGO Livingston (L1) and LIGO Hanford (H1)
detectors’ data, evaluating the maximum of the energy (“energy peak”), the
energy density above a given energy threshold defined in eq. (30) and the
corresponding TF area. It is worth mentioning that even if obtained with the
wavelet Q–transform our TF map looks very similar to those made with the
conventional Q–transform because in these representations the only difference
is a phase term that disappears when evaluating the local signal energy
$|T_{\mathrm{nd}}(\tau,\nu,Q,p)|^{2}$.

We estimate the GW waveforms using the denoising formula for both wavelet
transforms, eqs. (21), (22), and eqs. (12), (13). In each case we evaluate the
residuals between data and the reconstructed waveforms: a good reconstruction
is expected to produce a Gaussian distribution of the residuals, due to the
whitening performed before the analysis. The resulting TF maps and the
reconstructed waveforms for the Hanford interferometer are shown in Figure 5.
The Supplementary Text includes an extended version of this Figure and a
summary table with evaluations of the performance of the transforms.

The wavelet Qp–transform produces a considerable increase in both peak energy
and energy density with respect to the wavelet Q–transform: the increase of
the energy density is partly due to the corresponding lowering of the TF–area,
meaning that the wavelet Qp–transform produces a more compact TF
representation of the signal, as expected. Note that these values depend on
the threshold chosen when calculating the energy density. Here, the energy
threshold has been set to $E_{\mathrm{thr}}=7$: indeed, for a Gaussian noise
background, we expect only $0.1\%$ of the transform values to be above that
threshold, as shown in the Supplementary Text.

The waveforms (bottom panel in Fig. 5) have been obtained applying the
denoising formula to the corresponding wavelet transform: in each case, pixels
have been selected if
$|T_{\mathrm{nd}}(\tau,\nu,Q,p)|^{2}>E_{\mathrm{thr}}=7$. The signal shapes
obtained with the wavelet Q– and Qp–transforms are comparable to those
obtained with LALInference, a pipeline based on matched–filtering and Bayesian
estimation [55], and reported in [5, 4]. We evaluate the similarity between
the waveforms reconstructed with the wavelet Q– and Qp–transforms and the LAL
reconstruction with the overlap:

$O\left[s_{\mathrm{Qp/Q}},s_{\mathrm{LAL}}\right]=\dfrac{(s_{\mathrm{Qp/Q}},s_{\mathrm{LAL}})}{\sqrt{(s_{\mathrm{Qp/Q}},s_{\mathrm{Qp/Q}})\cdot(s_{\mathrm{LAL}},s_{\mathrm{LAL}})}},$
(31)

where the subscripts denote the kind of reconstruction and
$(a,b)=\int_{t_{rs}}a(t)b(t)dt$ is the scalar product calculated over the time
duration of the reconstructed signal. Both transforms return an overlap very
close to $1$ (the ideal match), as reported in Table 1 in the Supplementary
Text: this is a remarkable result, considering that we used neither templates
as in modeled analyses [55] nor the likelihood formalism of unmodeled analysis
[30]. The effectiveness of the denoising procedure is highlighted by the
histograms of the residuals, which have a distribution close to Gaussian for
both transforms and is corroborated by the values of the standard deviations
reported in Table 1 in the Supplementary Text. Overall, there is no major
difference between the waveforms obtained from wavelet Qp– and wavelet
Q–transform: this is probably due to the fact that the signal energy is
sufficiently high to be mostly above threshold for both transforms. This is
not the case for the Virgo data of GW190521, which we discuss in the next
subsection.

Figure 5: Q– and Qp–reconstruction of the first event GW150914, using data
from H1. The top panel shows the wavelet Qp–transform and the middle panel
shows the wavelet Q–transform; the colorbars display the energy scale, i.e.,
the value of $|T_{\mathrm{nd}}(\tau,\nu,Q,p)|^{2}$. The bottom panel shows the
strain plot with the original data (green), the waveform reconstructed with
LALInference (blue) (from [5]), and the waveform from the denoising formula
for the wavelet Qp–transform (red), eqs. (21) and (22).

###  IV.2 GW190521

The analysis for GW190521 [11] is performed in the same way as done for
GW150914, with the addition of Virgo (V1) data. In this case, Figure 6 shows
the TF map and the reconstructed waveform obtained with the wavelet
Qp–transform and V1 data; an extended version of the Figure can be found in
the Supplementary Text, which also includes a quantitative summary of the
results in Table 2.

This event differs from the previous one in many ways, it has a lower
signal–to–noise ratio, its duration is shorter and the bandwidth is limited to
lower frequencies, and it has been detected by a three–interferometer network.
V1 is the least sensitive interferometer, and for this reason we lowered the
denoising threshold to $E_{\mathrm{thr}}=5$, mainly to avoid cutting too large
a portion of the signal. Still, this threshold is sufficiently high to
effectively remove Gaussian noise, since in that case only $0.7\%$ of the
transform values surpass this energy threshold. As a general rule, the choice
of the energy threshold can be made by evaluating the energy distribution of
the transform values, identifying the signal by distinguishing it from the
$\chi^{2}$ distribution of the Gaussian noise, and finally choosing the
threshold by setting a reasonable trade–off between the removal of Gaussian
noise and the preservation of the GW signal.

The complete results are shown in Figure 3 in the Supplementary Text, where
the thresholds for L1 and H1 have been kept at $E_{\mathrm{thr}}=7$ as for
GW150914. In general, we note that the wavelet Qp–transform always performs
better than the wavelet Q–transform, as reported in Table 2 in the
Supplementary Text. For Virgo the improvement of the wavelet Qp–transform with
respect to the wavelet Q–transform is more marked and it is important to note
that results are considerably improved by careful filtering, even though the
threshold has been lowered with respect to that used for the other detectors.
With the wavelet Qp–transform several pixels pick up a higher energy density
and cross the energy threshold. Overall, the reduction of the TF area due to
the compactness of the wavelet Qp–transform representation is overcompensated
by the larger number of above–threshold pixels, leading to a noticeable
increase of the TF area of the wavelet Qp–transform with respect to the
wavelet Q–transform. In other words, with the wavelet Q–transform many pixels
that are likely to originate from the signal remain below threshold, while
with the wavelet Qp–transform they move above threshold, a clear improvement
due to the introduction of the $p$ parameter, and overall, the energy density
also increases. We remark that the wavelets of both transforms are
energy–normalized, therefore any difference in energy is given only by the
better/worse match of the signal with the wavelets.

The waveforms obtained with the denoising formulas applied to the TF map
obtained with both wavelet transforms are quite good for L1 and H1 and
comparable with those obtained with LALInference and reported in [11], with a
high value of the overlap. The low energy values for both transforms affect
the Virgo reconstruction, and in that case the overlap with LAL reconstruction
is lower. However, the wavelet Qp–transform has a larger number of
above–threshold pixels and a much better waveform reconstruction, resulting in
a considerably higher overlap with the LAL reconstruction with respect to that
of Q–transform.

Figure 6: Q– and Qp–reconstruction of GW190521, using data from V1. The top
panel shows the wavelet Qp–transform and the middle panel shows the wavelet
Q–transform; the colorbars display the energy scale, i.e., the value of
$|T_{\mathrm{nd}}(\tau,\nu,Q,p)|^{2}$. The bottom panel shows the strain plot
with the original data (green), the waveform reconstructed with LALInference
(blue) (from [11]), and the waveform from the denoising formula for the
wavelet Qp–transform (red), eqs. (21) and (22).

###  IV.3 GW170817 and glitch denoising

A peculiarity of the GW170817 detection [8] is that about a second before the
merger a loud glitch was observed in the L1 data: in the analysis performed by
the LIGO–Virgo collaboration the glitch was successfully removed with noise
mitigation techniques [4, 53, 9]. Here, we focus on the removal of this glitch
using the Qp–wavelet denoising formula. We do so by performing first a wavelet
Qp–transform of the original data by using the combination of $Q$ and $p$
values giving the best representation of the GW signal: then, we apply the
denoising formula setting to zero all pixels above an energy threshold of
$25$. Again, this threshold has been chosen by evaluating the transform energy
distribution, identifying both the signal and the glitch distributions, and
finally taking a threshold that does not cut off the signal while removing the
glitch only. In this way, we obtain a denoised TF map and a filtered time
series, for which we compute again the wavelet Qp–transform with the same $Q$
and $p$ pair maximizing the energy density of the GW signal.

Our results are shown in the two panels of Figure 7: the wavelet Qp–transform
of the filtered series (bottom panel) shows that the glitch (well visible in
the upper panel) has been almost completely removed.

Figure 7: Wavelet Qp–transform of the L1 data for GW170817. The upper panel
shows the wavelet Qp–transform of the original whitened data which include the
glitch. The lower panel shows the filtered TF map, where the glitch has been
almost completely removed.

##  V Conclusions

In this paper we presented important improvements to the Q–transform. We
extended the Q–transform to the wavelet Q–transform — a parameterized version
of the Morlet wavelet transform — which does have an inversion formula, unlike
the Q–transform. This result extends the applicability of the formula
originally obtained by Lebedeva and Postnikov [31], leading to an original
denoising algorithm which is both effective and computationally efficient.
While the wavelet Q–transform, just as the Q–transform, has a minimal
uncertainty in the time-frequency plane, it does not have an optimal
performance in the case of chirping signals like those produced by CBCs. We
found a useful variant — which we dubbed the wavelet Qp–transform — which can
still be inverted and produces sparser representations of the
gravitational–wave signals from CBCs. We have fully characterized these
wavelets, both in the time and in the frequency domain, also studying their
statistical fluctuations and correlation properties in presence of a Gaussian
white noise background. We have produced a Python implementation of both
transforms and used it to analyze some important gravitational–wave signals
detected by the LIGO/Virgo Collaboration during the O1 and O2 observing runs.
These tests illustrate the higher efficiency of the wavelet Qp–transform for
chirping signals, and also the excellent performance of the transform as a
denoising tool.

The results presented here are encouraging, even though we have tested the
transforms on a very small set of GW events, and they suggest more systematic
studies involving a larger number of events or a set of simulations.

There remain several unanswered questions. The most pressing is how to combine
the data streams of several gravitational–wave detectors exploiting their
coherence.

Finally, we wish to remark that the inversion formula described in this paper,
in particular its denoising version of eqs. (12) and (13), which we developed
for the analysis of GW signals, can be applied to many other fields like
music, medicine, geophysics, engineering and in general in all those fields
which require the analysis of noisy, transient signals.

##  VI Code availability

A first implementation of the wavelet Q– and Qp–transforms in Python,
following the guidelines provided by the GWpy module qtransform.py [33], is
freely available at the link https://zenodo.org/doi/10.5281/zenodo.10649072,
including as an example the analysis of GW150914.

## Acknowledgements

We would like to thank the LIGO–Virgo groups of Trieste, Trento, Roma and
Zürich for useful discussions. This research has made use of data or software
obtained from the Gravitational Wave Open Science Center (http://gwosc.org), a
service of LIGO Laboratory, the LIGO Scientific Collaboration, the Virgo
Collaboration, and KAGRA. LIGO Laboratory and Advanced LIGO are funded by the
United States National Science Foundation (NSF) as well as the Science and
Technology Facilities Council (STFC) of the United Kingdom, the
Max–Planck–Society (MPS), and the State of Niedersachsen/Germany for support
of the construction of Advanced LIGO and construction and operation of the
GEO600 detector. Additional support for Advanced LIGO was provided by the
Australian Research Council. Virgo is funded, through the European
Gravitational Observatory (EGO), by the French Centre National de Recherche
Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN)
and the Dutch Nikhef, with contributions by institutions from Belgium,
Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain.
KAGRA is supported by Ministry of Education, Culture, Sports, Science and
Technology (MEXT), Japan Society for the Promotion of Science (JSPS) in Japan;
National Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in
Korea; Academia Sinica (AS) and National Science and Technology Council (NSTC)
in Taiwan. EM acknowledges support by ICSC – Centro Nazionale di Ricerca in
High Performance Computing, Big Data and Quantum Computing, funded by European
Union – NextGenerationEU.

## References

  * [1] B. Abbott, R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, R. Adhikari, V. Adya, C. Affeldt, et al.  GWTC-1: a gravitational-wave transient catalog of compact binary mergers observed by LIGO and Virgo during the first and second observing runs.  Physical Review X, 9(3):031040, 2019.
  * [2] B. Abbott, R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, R. X. Adhikari, V. Adya, C. Affeldt, et al.  All-sky search for short gravitational-wave bursts in the second advanced ligo and advanced virgo run.  Physical Review D, 100(2):024017, 2019.
  * [3] B. Abbott, R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, V. Adya, C. Affeldt, M. Agathos, et al.  Optically targeted search for gravitational waves emitted by core-collapse supernovae during the first and second observing runs of advanced ligo and advanced virgo.  Physical Review D, 101(8):084002, 2020.
  * [4] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, et al.  Gw150914: First results from the search for binary black hole coalescence with advanced ligo.  Physical Review D, 93(12):122003, 2016.
  * [5] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, et al.  Observation of gravitational waves from a binary black hole merger.  Physical review letters, 116(6):061102, 2016.
  * [6] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, et al.  Observing gravitational-wave transient GW150914 with minimal assumptions.  Physical Review D, 93(12):122004, 2016.
  * [7] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, et al.  All-sky search for short gravitational-wave bursts in the first advanced ligo run.  Physical Review D, 95(4):042003, 2017.
  * [8] B. P. Abbott, R. Abbott, T. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, V. B. Adya, et al.  GW170817: observation of gravitational waves from a binary neutron star inspiral.  Physical review letters, 119(16):161101, 2017.
  * [9] B. P. Abbott, R. Abbott, T. D. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, et al.  Effects of data quality vetoes on a search for compact binary coalescences in Advanced LIGO’s first observing run.  Classical and Quantum Gravity, 35(6):065010, 2018.
  * [10] R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, A. Adams, C. Adams, R. Adhikari, V. Adya, C. Affeldt, et al.  GWTC-2: compact binary coalescences observed by LIGO and Virgo during the first half of the third observing run.  Physical Review X, 11(2):021053, 2021.
  * [11] R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, R. Adhikari, V. Adya, C. Affeldt, M. Agathos, et al.  GW190521: a binary black hole merger with a total mass of 150 $M_{\odot}$.  Physical review letters, 125(10):101102, 2020.
  * [12] R. Abbott, T. Abbott, F. Acernese, K. Ackley, C. Adams, N. Adhikari, R. Adhikari, V. Adya, C. Affeldt, D. Agarwal, et al.  GWTC-3: compact binary coalescences observed by LIGO and Virgo during the second part of the third observing run.  Phys. Rev. X, 13:041039, Dec 2023.
  * [13] R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, R. X. Adhikari, V. B. Adya, C. Affeldt, M. Agathos, et al.  Open data from the first and second observing runs of Advanced LIGO and Advanced Virgo.  SoftwareX, 13:100658, 2021.
  * [14] R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, N. Adhikari, R. Adhikari, V. Adya, C. Affeldt, D. Agarwal, et al.  All-sky search for short gravitational-wave bursts in the third advanced ligo and advanced virgo run.  Physical Review D, 104(12):122004, 2021.
  * [15] R. Abbott, H. Abe, F. Acernese, K. Ackley, S. Adhicary, N. Adhikari, R. Adhikari, V. Adkins, V. Adya, C. Affeldt, et al.  Open data from the third observing run of LIGO, Virgo, KAGRA and GEO.  arXiv preprint arXiv:2302.03676, 2023.
  * [16] P. Addison, J. Watson, and T. Feng.  Low-oscillation complex wavelets.  Journal of Sound and Vibration, 254(4):733–762, 2002.
  * [17] F. Aubin, F. Brighenti, R. Chierici, D. Estevez, G. Greco, G. M. Guidi, V. Juste, F. Marion, B. Mours, E. Nitoglia, et al.  The MBTA pipeline for detecting compact binary coalescences in the third LIGO–Virgo observing run.  Classical and Quantum Gravity, 38(9):095004, 2021.
  * [18] J. C. Brown.  Calculation of a constant Q spectral transform.  The Journal of the Acoustical Society of America, 89(1):425–434, 1991.
  * [19] E. J. Candes, P. R. Charlton, and H. Helgason.  Detecting highly oscillatory signals by chirplet path pursuit.  Applied and Computational Harmonic Analysis, 24(1):14–40, 2008.
  * [20] E. Chassande-Mottin and A. Pai.  Best chirplet chain: Near-optimal detection of gravitational wave chirps.  Physical review D, 73(4):042003, 2006.
  * [21] S. Chatterji, L. Blackburn, G. Martin, and E. Katsavounidis.  Multiresolution techniques for the detection of gravitational-wave bursts.  Classical and Quantum Gravity, 21(20):S1809, 2004.
  * [22] S. K. Chatterji.  The search for gravitational wave bursts in data from the second LIGO science run.  PhD thesis, Massachusetts Institute of Technology, 2005.
  * [23] N. J. Cornish, T. B. Littenberg, B. Bécsy, K. Chatziioannou, J. A. Clark, S. Ghonge, and M. Millhouse.  BayesWave analysis pipeline in the era of gravitational wave observations.  Physical Review D, 103(4):044006, 2021.
  * [24] T. Dal Canton, A. H. Nitz, B. Gadre, G. S. C. Davies, V. Villa-Ortega, T. Dent, I. Harry, and L. Xiao.  Real-time Search for Compact Binary Mergers in Advanced LIGO and Virgo’s Third Observing Run Using PyCBC Live.  The Astrophysical Journal, 923(2):254, 2021.
  * [25] I. Daubechies.  Ten lectures on wavelets.  SIAM, 1992.
  * [26] M. Drago, S. Klimenko, C. Lazzaro, E. Milotti, G. Mitselmakher, V. Necula, B. O’Brian, G. A. Prodi, F. Salemi, M. Szczepanczyk, et al.  Coherent WaveBurst, a pipeline for unmodeled gravitational-wave data analysis.  SoftwareX, 14:100678, 2021.
  * [27] D. Gabor.  Theory of communication. Part 1: The analysis of information.  Journal of the Institution of Electrical Engineers-part III: radio and communication engineering, 93(26):429–441, 1946.
  * [28] C. Henshaw, M. Arogeti, A. Heranval, and L. Cadonati.  Visualization of frequency structures in gravitational wave signals.  arXiv preprint arXiv:2402.16533, 2024.
  * [29] D. Jordan, R. Miksad, and E. Powers.  Implementation of the continuous wavelet transform for digital time series analysis.  Review of scientific instruments, 68(3):1484–1494, 1997.
  * [30] S. Klimenko, G. Vedovato, M. Drago, F. Salemi, V. Tiwari, G. Prodi, C. Lazzaro, K. Ackley, S. Tiwari, C. Da Silva, et al.  Method for detection and reconstruction of gravitational wave transients with networks of advanced detectors.  Physical Review D, 93(4):042004, 2016.
  * [31] E. A. Lebedeva and E. B. Postnikov.  On alternative wavelet reconstruction formula: a case study of approximate wavelets.  Royal Society Open Science, 1(2):140124, 2014.
  * [32] L. Li, N. Han, Q. Jiang, and C. K. Chui.  A chirplet transform-based mode retrieval method for multicomponent signals with crossover instantaneous frequencies.  Digital Signal Processing, 120:103262, 2022.
  * [33] D. M. Macleod, J. S. Areeda, S. B. Coughlin, T. J. Massinger, and A. L. Urban.  GWpy: A Python package for gravitational-wave astrophysics.  SoftwareX, 13:100657, 2021.
  * [34] S. Mallat.  A Wavelet Tour of Signal Processing: The Sparse Way.  Academic Press, Inc., 2008.
  * [35] S. Mann and S. Haykin.  The chirplet transform: A generalization of Gabor’s logon transform.  In Vision interface, volume 91, pages 205–212. Citeseer Princeton, NJ, USA, 1991.
  * [36] S. Mann and S. Haykin.  The chirplet transform: Physical considerations.  IEEE Transactions on Signal Processing, 43(11):2745–2761, 1995.
  * [37] M. Millhouse, N. J. Cornish, and T. Littenberg.  Bayesian reconstruction of gravitational wave bursts using chirplets.  Physical Review D, 97(10):104057, 2018.
  * [38] V. Necula, S. Klimenko, and G. Mitselmakher.  Transient analysis with fast Wilson-Daubechies time-frequency transform.  In Journal of Physics: Conference Series, volume 363, page 012032. IOP Publishing, 2012.
  * [39] A. Pai, E. Chassande-Mottin, and O. Rabaste.  Best network chirplet chain: Near-optimal coherent detection of unmodeled gravitational wave chirps with a network of detectors.  Physical Review D, 77(6):062005, 2008.
  * [40] E. B. Postnikov, E. A. Lebedeva, and A. I. Lavrova.  Computational implementation of the inverse continuous wavelet transform without a requirement of the admissibility condition.  Applied Mathematics and Computation, 282:128–136, 2016.
  * [41] P. Relton, A. Virtuoso, S. Bini, V. Raymond, I. Harry, M. Drago, C. Lazzaro, A. Miani, and S. Tiwari.  Addressing the challenges of detecting time-overlapping compact binary coalescences.  Physical Review D, 106(10):104045, 2022.
  * [42] O. Rioul and M. Vetterli.  Wavelets and signal processing.  IEEE signal processing magazine, 8(4):14–38, 1991.
  * [43] F. Robinet, N. Arnaud, N. Leroy, A. Lundgren, D. Macleod, and J. McIver.  Omicron: a tool to characterize transient noise in gravitational-wave detectors.  SoftwareX, 12:100620, 2020.
  * [44] S. Roy.  Nonorthogonal wavelet transformation for reconstructing gravitational wave signals.  Physical Review Research, 4(3):033078, 2022.
  * [45] D. Slepian.  On bandwidth.  Proceedings of the IEEE, 64(3):292–300, 1976.
  * [46] D. Slepian.  Some comments on Fourier analysis, uncertainty and modeling.  SIAM review, 25(3):379–393, 1983.
  * [47] R. G. Stockwell, L. Mansinha, and R. Lowe.  Localization of the complex spectrum: the S transform.  IEEE transactions on signal processing, 44(4):998–1001, 1996.
  * [48] M. Szczepańczyk, S. Klimenko, B. O’Brien, I. Bartos, V. Gayathri, G. Mitselmakher, G. Prodi, G. Vedovato, C. Lazzaro, E. Milotti, et al.  Observing an intermediate-mass black hole GW190521 with minimal assumptions.  Physical Review D, 103(8):082002, 2021.
  * [49] M. J. Szczepańczyk, J. M. Antelis, M. Benjamin, M. Cavaglia, D. Gondek-Rosińska, T. Hansen, S. Klimenko, M. D. Morales, C. Moreno, S. Mukherjee, et al.  Detecting and reconstructing gravitational waves from the next galactic core-collapse supernova in the advanced detector era.  Physical Review D, 104(10):102002, 2021.
  * [50] M. J. Szczepańczyk, F. Salemi, S. Bini, T. Mishra, G. Vedovato, V. Gayathri, I. Bartos, S. Bhaumik, M. Drago, O. Halim, et al.  Search for gravitational-wave bursts in the third Advanced LIGO-Virgo run with coherent WaveBurst enhanced by machine learning.  Physical Review D, 107(6):062002, 2023.
  * [51] M. J. Szczepańczyk, Y. Zheng, J. M. Antelis, M. Benjamin, M.-A. Bizouard, A. Casallas-Lagos, P. Cerdá-Durán, D. Davis, D. Gondek-Rosińska, S. Klimenko, et al.  An Optically Targeted Search for Gravitational Waves emitted by Core-Collapse Supernovae during the Third Observing Run of Advanced LIGO and Advanced Virgo.  arXiv preprint arXiv:2305.16146, 2023.
  * [52] L. Tsukada, P. Joshi, S. Adhicary, R. George, A. Guimaraes, C. Hanna, R. Magee, A. Zimmerman, P. Baral, A. Baylor, et al.  Improved ranking statistics of the GstLAL inspiral search for compact binary coalescences.  arXiv preprint arXiv:2305.06286, 2023.
  * [53] S. A. Usman, A. H. Nitz, I. W. Harry, C. M. Biwer, D. A. Brown, M. Cabero, C. D. Capano, T. Dal Canton, T. Dent, S. Fairhurst, et al.  The PyCBC search for gravitational waves from compact binary coalescence.  Classical and Quantum Gravity, 33(21):215004, 2016.
  * [54] G. Vedovato, E. Milotti, G. Prodi, S. Bini, M. Drago, V. Gayathri, O. Halim, C. Lazzaro, D. Lopez, A. Miani, et al.  Minimally-modeled search of higher multipole gravitational-wave radiation in compact binary coalescences.  Classical and Quantum Gravity, 39(4):045001, 2022.
  * [55] J. Veitch, V. Raymond, B. Farr, W. Farr, P. Graff, S. Vitale, B. Aylott, K. Blackburn, N. Christensen, M. Coughlin, et al.  Parameter estimation for compact binaries with ground-based gravitational-wave observations using the LALInference software library.  Physical Review D, 91(4):042003, 2015.
  * [56] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [57] Y. Yang, W. Zhang, Z. Peng, and G. Meng.  Multicomponent signal analysis based on polynomial chirplet transform.  IEEE Transactions on Industrial Electronics, 60(9):3948–3956, 2012.
  * [58] G. Yu and Y. Zhou.  General linear chirplet transform.  Mechanical Systems and Signal Processing, 70:958–973, 2016.