Bertille Martinez

Julien Ramade

Clément Livache

Nicolas Goubet

Audrey Chu

Charlie Gréboval

Junling Qu

William L Watkins

Loïc Becerra

Erwan Dandeu

Jean Louis Fave

Christophe Méthivier

Emmanuelle Lacaze

Emmanuel Lhuillier

HgTe Nanocrystal Inks for Extended Short Wave Infrared Detection

Keywords: HgTe, nanocrystals, short wave infrared, photodiode, focal plane array

Short wave infrared (IR) detection is currently driven by InGaAs technology which has limited the perspective of cost effectiveness and consequently slowed the development of IR sensors. Since organic electronics are ineffective in this wavelength range, an alternative to conductive polymers is the use of colloidal quantum dots (CQDs) which exhibit strongly tunable IR absorption. In this paper, we focus on the extended short-wave IR (2.5 µm cut-off) to expand the capabilities of InGaAs, while using HgTe nanocrystals as the active material. Previous devices based on this material suffered from a low responsivity due to weak absorption (few percents). We present the integration of HgTe nanocrystals in ink form to build thick (up to 600 nm), strongly absorbing nanoparticle films. This ink is integrated into a diode, allowing one to boost the responsivity by two orders of magnitude and the detectivity by one order of magnitude compared to previous HgTe devices at the same wavelength. Detectivity reaches 3 × 10 9 Jones, while the time response is found to be 370 ns for room temperature and 0 V bias operation. Finally, the material is integrated into a focal plane array which is used to determine laser beam profile.

INTRODUCTION

The interest for nanocrystals [1][START_REF] Konstantatos | Colloidal Quantum Dot Optoelectronics and Photovoltaics[END_REF][START_REF] Klimov | Nanocrystal Quantum Dots[END_REF] has been mostly driven by their luminescent properties in applications such as displays, [START_REF] Bourzac | [END_REF] single photon sources [5] or fluorescent imaging, [6] though their absorption features appear equally promising. In particular, nanocrystal absorption is a robust property (i.e. absorption barely photobleaches and is not strongly affected by a change in the dielectric environment or surface chemistry) which can be broadly tuned from UV to THz ranges. [7][8][9][10][11] This robustness combined with the possibility to design narrow energy transition absorption which, unlike organic electronics, enables nanocrystals to be used for infrared sensing.

Short wave infrared (SWIR: 1-3 µm) is a part of the electromagnetic spectrum invisible to the naked eye finding applications in night vision, LIDAR imaging or active imaging. Infrared sensing in the SWIR is currently led by InGaAs semiconductors. The high performances achieved by this material are currently strongly inhibited by its cost. Indeed, road maps for this technology only plan marginal cost reduction. InGaAs also suffers from a lack of spectral tunability, especially when wavelengths longer than 1.7 µm are required. Significant efforts have been dedicated in the design of alternative infrared sensors based on colloidal nanocrystals. PbS has demonstrated significant performances in the near IR and up to telecom wavelengths. [12,13] Nevertheless, above 1.5 µm, quite large nanoparticles have to be grown exhibiting lower colloidal stability, [14] which is highly detrimental for the design of pinhole-free vertical geometry devices. An alternative is to use narrower band gap semiconductors and semimetals. [15,16] Under colloidal form, HgTe is certainly the most mature material combining infrared absorption and photoconduction. [17] On the other hand, efforts dedicated to the integration of HgTe nanocrystals into photodiodes [18,19] are much fewer than those reported for PbS. As a result, most devices based on HgTe nanocrystals with long wavelength absorption are based on planar geometries. [15,[20][21][22][23][24] Recently, mid-wave infrared photodiodes [25] , with high photoresponse [26] (up to 0.8 A/W), have been proposed by Guyot-Sionnest's group. Part of this success is based on the use of plasmon resonance to enhance the absorption of thin nanocrystal films with thicknesses below the absorption depth. [20,21,26,27] Yet, their 5 µm cut-off wavelength means that these devices still need to be cooled down to cryogenic temperatures to achieve a reasonable signal/noise ratio. In the SWIR and extended SWIR (up to 2.5 µm), targeted performances need to be achieved without cooling or with a Peltier system at most. Thus, new designs must be proposed. In this wavelength range, Jagtap et al [28,29] have proposed the use of a basic diode structure inspired by solar cells : a solid-state ligand exchanged HgTe layer is coupled with a TiO2 electron transport layer which is used to induce the rectifying behavior. Later, it was proposed to add a unipolar barrier which had the purpose to selectively decrease the electron dark current. Room temperature detectivity up to 3 × 10 8 Jones has been demonstrated. [23] In this device, the performance remains clearly limited by the low responsivity (100 µA/W at 0 V) rather than by the signal/noise ratio. Indeed, the TiO2 layer, with its 3.2 eV band gap used for the diode, presents a poor band alignment with HgTe nanocrystals. As a result, not only is the dark hole current filtered but so is a significant part of the photocurrent. Furthermore, only thin films (<200 nm) of HgTe have been used as absorbing layers leading to a low absorption. The absorption coefficient [30] α of HgTe nanocrystal with a cut-off at 2.5 µm has been estimated to be 3±1 × 10 3 cm -1 , thus we can estimate the absorption of a 200 nm thin film to be 1-10 -𝛼𝐿 ≈ 12%. Indeed, solid state ligand exchange, which is generally used for HgTe nanocrystals devices, makes it challenging to fabricate thick, smooth and strongly electronically coupled films. In the case of PbS, the recent trend is to switch from the deposition of nanocrystals in an apolar solvent to nanocrystal inks [31][32][33] where ligand exchange is undertaken through a phase transfer. [34] This effort is still quite marginally used for HgTe. Here, we explore the potential of a HgTe nanocrystal-based ink for the design of extended SWIR sensors. The capability of this ink is then tested for the design of a photodiode and a multipixel device. We particularly demonstrate the capability of this focal plane array to determine the infrared laser beam shape.

DISCUSSION

At first, HgTe nanocrystals were synthesized with an excitonic feature at 4000 cm -1 (2.5 µm or 0.5 eV) using the Keuleyan's procedure. [35] The nanoparticles, later named HgTe 4k, undertake a tetrapodic shape, as shown in the transmission electron microscopy (TEM) image in Figure 1a. Secondly, an ink was prepared from this solution. Since mercury has a strong affinity with thiols, a short thiol (ethanedithiol or mercapto alcohol) was used as a transfer agent to disperse the nanocrystals in the polar phase solvent dimethyl formamide (DMF) from the initial toluene solvent. The addition of HgCl2 was found to further enhance electronic passivation and allowed good conductivity. [25,[36][37][38] We observed that the nanoparticle spectrum was barely affected by this ligand exchange procedure, see Figure 1b. This contrasts with the solid state ligand exchange (SSLE) based on ethanedithiol (EDT) which led to a redshift and broadening of the band edge transition. [39] Spin coating of this solution at low speeds (<1000 rpm) yielded thick (up to 600 nm) and smooth films of HgTe, with a roughness of the order of the nanocrystal size, see Figure 1c. A vacuum drying step was then required to remove all the left over DMF. S4b andS4c).

Figure 1 : a) TEM image of HgTe 4k used in this study. b) Infrared spectra of HgTe 4k as-synthetized (DDT ligands, toluene, in red) and ligand-exchanged (ink, DMF, in black). c) Profilometry of a HgTe ink film, insert: macroscopic image of HgTe ink film. The film the quality of the film is highlighted by its reflectivity. d) Transfer curves measured for ink device (in black) and EDT on film ligand exchange device (in red). The two curves were normalized by their value at + 2 V. The current leakage was at least one order of magnitude below the channel current (see Figure

The transport properties of the layer, probed in an electrolyte gated field effect transistor configuration, [40] shows that the ambipolar character of the film is mostly unaffected, see Figure 1d and Figure S4. However the mobility of the film is significantly improved, see Figure S5. The determination of the mobility was done using conventional back gating through a SiO2 layer. We determined the electron mobility to be 3 × 10 -2 cm 2 V -1 s -1 , while the mobility of HgTe nanocrystals capped with EDT was previously determined [39] to be 10 -4 cm 2 V -1 s -1 . This is indeed promising to achieve longer diffusion lengths and collect photocarriers over a thicker layer.

In addition, the ink exhibits an enhanced stability compared to SSLE prepared films, see Figure 2a. Indeed, the resistance of a film stored in ambient conditions over a week only dropped by a factor two, while a factor 1000 drop was observed for a SSLE prepared film. With the SSLE process, coating the device with water and oxygen repealing layers was mandatory to achieve stable performances. [28] Understanding the increase in dark current, gives insight on how the optical feature diminishes over time. The film was annealed at 120 °C, to mimic clean room processing. This temperature is similar to what it is expected to be exposed during the lithography process or atomic layer deposition process. With conventional SSLE with EDT, we observed a severe redshift of the exciton, see Figure 2c, while the energy of the band edge remains almost unaffected for the ink, see Figure 2b. Indeed, upon aging and annealing, HgTe nanocrystals prepared by SSLE sinter which reduces the confinement and the optical band gap. As the dark current is exponentially dependent on the band gap, the reduced band gap leads to an increase of the dark current.

Films prepared by SSLE EDT and inks were also analyzed before and after aging using photoemission spectroscopy. As the ratio Te/Hg is constant in the case of ink based films, it increases by 63 % for SSLE EDT films. On these SSLE films, we see a new contribution of Te 3d core levels at 576.1 eV appearing after 72 hours aging, see Figure 2e. This contribution can be attributed to the formation of TeO2. [41] Mercury core levels were also investigated, but neither aging nor preparation methods (SSLE or ink) had a significant influence, see Figure 2d.

In the case of SSLE films, at low binding energies, the valence band shifts towards the Fermi level, see Figure S3b andc. This shift leads to a p-type doping consistent with the increase of Te/Hg ratio. On films prepared by ink, the energy difference between the valence band and Fermi level stays constant after 72 hours of aging and no oxidation of tellurium was observed. The p-doping of the nanocrystals also participates in the increase of dark current in the case of SSLE preparation. This thick absorbing layer was then integrated into a photodiode, see Figure 3a-c. The diode has the following structure: glass/FTO/HgTe 4k/HgTe 6k/Au. In the extended SWIR, the absorption of conventional transparent conductive oxide (TCO) starts to be non-negligible. To avoid an excessive absorption, a thin (≈ 80 nm) fluorinated tin oxide (FTO) is used, its transmission at the band edge of the nanocrystals being above 60%, see figure S6. A thick layer (≈ 600 nm) of HgTe 4k is used as an absorbing material, see Figure 3b andc. The diode rectifying behavior is introduced by a unipolar barrier made of HgTe nanocrystals with a wider band gap as proposed by Jagtap et al. [29] The role of the unipolar barrier is to let photocurrent flow while preventing the dark current injection for the other carrier [39][40][41].

To reach this goal, the barrier's and absorbing material's electronic levels needed to be aligned for the carriers to flow (in this case, holes), while the electronic levels need to be offset for the carriers to be filtered (i.e. electon). Such barrier can be obtained with smaller HgTe nanocrystals exhibiting a band edge energy of ≈720 meV (6000 cm -1 -later called HgTe 6k). Ultraviolet photoemission (Figure S1d) reveals a near resonance between their valence band and that of the absorbing layer made of HgTe 4k, while the conduction band is offset by more than 200 meV. It is also worth noticing that the HgTe 6k CQDs display only a p-type conduction (Figure S1) which further favors electron filtering. A top gold contact is used as a hole extractor.

The I-V curve reveals a rectifying behavior even at room temperature (see Figure 3d), with a weak dark current (<10 -6 A.cm -2 ) around 0 V. Under black body illumination (65 µW, corresponding to integrated power from 1.9 µm to the HgTe 4k cut-off), the responsivity of the device has been determined to be 20 mA.W -1 , see Figure 3e and S9. This is almost two orders of magnitude higher than previously reported devices at the same wavelength based on HgTe nanocrystals, and only diodes with optical resonators [19] have achieved higher responsivities. This is the result of (i) a higher absorption thanks to a thicker HgTe film. A L = 600 nm thick film is expected to absorb 1-10 -𝛼𝐿 ≈ 33% of light for each time it passes through the film. Experimentally, we determine that a film of 600±20 nm has an optical density of 0.23 at 4000 cm -1 m (2.5 µm), which corresponds to a light absorption of 40 %, see Figure S2. (ii) Better photocharge extraction due to the removal of the TiO2 layer. For the same diode structure, we observe a drop by four orders of magnitude of the responsivity in the presence of the TiO2 electron transport layer, see Figure S7. It is also worth noting that the frequency response of the device is flat in the 10 Hz to 1 kHz range, suggesting a fast response time (see Figure 3e). [45] The time response has actually been determined to be τ ≈ 370 ns under nanosecond pulse excitation (Figure 2g), this corresponds to a cut-off frequency of 𝑓 𝑐 = 1 2𝜋𝜏 = 430 kHz.

From this long lifetime, we can estimate the carrier diffusion length LD [46] to be given by LD= √ µ𝑘 𝑏 𝑇𝜏 𝑒 with µ the mobility, kb the Boltzman constant, T the temperature and e the proton charge. Using the value of 10 -4 cm 2 V -1 s -1 , as measured for EDT capped HgTe nanocrystals, we obtain a 10 nm diffusion length, barely larger than the nanocrystal size. Fabricating thick film is not useful in this condition since photocarriers will not be collected. On the other hand, the ink, with its 300 times higher mobility, leads to diffusion length of ≈170 nm, which can lead to an efficient charge collection over the 500 nm thick diode film.

Again, for nanocrystal based devices, the noise is found to be 1/f limited, [47][48][49] (see the inset of Figure 3f and Figure S10). The detectivity is calculated as follows:

𝐷 * = 𝑅√𝐴 𝑖 𝑛
with R the responsivity in A.W -1 , A the optical area of the device and in the current spectral density in A.Hz -1/2 . The measured detectivity reaches 3 × 10 9 Jones at room temperature, under 0 V and for detection at 1 kHz (see Figure 2f).

b) SEM cross-section of the diode. c) Band alignment of the different components of the diode. d) IV curves in the dark and under illumination by a 1.55 µm laser at different powers, insert: evolution of open circuit voltage (VOC) as a function of increasing power. e) Responsivity under black body illumination at 980 °C, corresponding to 65 µW incident power, at different frequencies. A Ge filter is used to remove photons with a wavelength below 1.9 µm, as well as to ensure that the responsivity only comes from HgTe 4k excitation. f) Detectivity as a function of frequency, insert: noise signal as a function of detection frequency (black line) and f -1/2 fit (red dash line). g) Time response of the diode measured with 355 nm pulsed laser (2 ns pulse). The time response value is obtained through single exponential fitting.

The operation of the diode and the unipolar barrier can be better appreciated thanks to photocurrent spectra (see Figure 4a). Since the unipolar barrier is made from an active material (HgTe 6k), photocurrent is likely to occur from different parts of the diode under broad band excitation. In the reverse bias regime, the HgTe 6k barrier prevents electron injection from the gold side (i.e., dark current) but does not prevent electronic current coming from HgTe 4k excitation (i.e., photocurrent) nor hole current (see the left inset of Figure 4b). This tends to maximize the photocurrent resulted from the HgTe 4k layer. Under low forward bias, a significant part of the electronic photocurrent coming from the HgTe 4k layer gets blocked by the barrier. A decrease of the relative weight of HgTe 4k in the photocurrent was observed (see the middle inset of Figure 4b). When a large bias is applied, the Stark effect increases and thus the transparency of the barrier and the photocurrent from the HgTe 4k layer can again flow in the system (see the right inset of Figure 4b). In this sense, the unipolar barrier plays exactly expected role as long as the device is operated under low voltage. This happens when the energy drops along the barrier e.F.L (where e is the elementary charge, F is the applied electric field and L is the length of the barrier, around 30 nm) is far below the band offset between HgTe 4k and HgTe 6k (≈ 200 meV) The diode capability for low flux detection has been tested by detecting a black body far away from the source (≈4 m, see Figure 5a andc). It is worth mentioning that extended SWIR diodes are not optimized for thermal imaging since only the photons in the abrupt rise of the black body spectrum between the Ge filter cutoff (1.9 µm to suppress visible photon) and the diode cut off wavelength (2.5 µm) get absorbed. At this distance, warm black bodies (>700 °C) can be clearly discriminated from the dark current (see Figure 5b). Another important characteristic, often overlooked for nanocrystal-based devices, is the angular response determination. The latter well follows the expected cosine dependency (see Figure 5d). As it may have been anticipated, though not confirmed, no angular selection rules are introduced by the HgTe nanocrystals or by any waveguide effect within the diode stacking. The last part of the paper is dedicated to the integration of this strongly absorbing film of HgTe within a multipixel geometry. HgTe, as all Hg containing materials, is relatively soft and the nanocrystals' growth temperature is relatively low (80 °C) making them easy to sinter when exposed to high temperatures. This leads to a redshift of the optical feature and a dramatic increase of the dark current. [28] The latter process can be partly prevented by using HgTe/CdTe core-shell systems. [50] Regardless, it remains that device geometries and fabrication steps need to account for the material's sensibility. Here, we demonstrate the feasibility of a 10 × 10 pixel matrix and its use for infrared laser beam shape determination. To limit the exposition of the HgTe nanocrystals to any annealing step, the electrodes are built a priori and the nanocrystals are deposited during the last step of fabrication. The matrix fabrication is described in Figure 6, S11 and S12. Typically, each pixel has a 27 µm size, with a 33 µm size, corresponding to a filling factor of 67% over the whole matrix. One should note that this pixel size is similar to what is currently used for InGaAs technology. Each pixel is made up of two metallic contacts separated by a SiO2 insulating barrier. The matrix' homogeneity can be seen from the dark current histogram, see Figure 7a. Out of 100 pixels, only 3 were shorted (even before the ink deposition). The full width at half maximum (FWMH) of the distribution is 45% of the average value. This may seem important, however, the fluctuations are mostly the result of the electrodes dispersion quality rather than film quality since we can spatially correlate the less resistive pixels to be in the vicinity of the shorted electrodes. We have then tested the capability of this focal plane array (FPA) for infrared beam shape determination. We have used a 1.55 µm laser diode source, with a beam linewidth determined to be around 60 µm, see Figure 7b. Once the laser is shone on the FPA, we recovered a similar linewidth 75 µm, see Figure 7c. Hence, the FPA is sensitive enough to detect higher order modes present in the laser beam, see Figure S13.

CONCLUSION

We have explored the possibility to use HgTe nanocrystal-based ink to design extended short wave infrared sensors with enhanced absorption. Compared to solid state ligand exchange film, the ink preparation preserves the optical spectrum and leads to higher electrical stability for in air operation. The material has been integrated into a diode, where the rectifying behavior is the result of a unipolar barrier. With this geometry, we have improved the responsivity of the film by two orders of magnitude and the detectivity by one order of magnitude compared to previous devices at the same wavelength and based on the same material. The detectivity now reaches 3 × 10 9 Jones at room temperature, while the time response can be as short as 370 ns. The material has also been tested as an active material of a focal plane array (FPA), with a pixel pitch down to 30 µm, with the fabrication designed to account for HgTe softness. We demonstrate the use of this FPA for infrared laser beam shape determination.

Figure 2 :

 2 Figure2: a) Evolution of resistance under 1 V bias for two photoconductive devices: ink device (in black) and EDT on film ligand exchange (in red). b) Absorption spectrum of a HgTe nanocrystal film obtained from an ink before and after annealing at 120 °C for 7 min. c) Absorption spectrum of a HgTe nanocrystal film capped with EDT before and after annealing at 120 °C for 7 min. d) Mercury 4f core levels of films prepared by EDT ligand exchange and by ink, just after their fabrication (t = 0), and after 72 hours spent in air (t = 72 h). e) Tellurium 3d core levels of films prepared by EDT ligand exchange and by ink, just after their fabrication (t = 0 ), and after 72 hours spent in air (t = 72 h). On EDT capped HgTe film, after 72 hours in air, we see a second contribution, in red on the signal, highlighted by the black dotted lines.

Figure 3 :

 3 Figure 3: a) Diagram of the diode structure. b) SEM cross-section of the diode. c) Band alignment of the different components of the diode. d) IV curves in the dark and under illumination by a 1.55 µm laser at different powers, insert: evolution of open circuit voltage (VOC) as a function of increasing power. e) Responsivity under black body illumination at 980 °C, corresponding to 65 µW incident power, at different frequencies. A Ge filter is used to remove photons with a wavelength below 1.9 µm, as well as to ensure that the responsivity only comes from HgTe 4k excitation. f) Detectivity as a function of frequency, insert: noise signal as a function of detection frequency (black line) and f -1/2 fit (red dash line). g) Time response of the diode measured with 355 nm pulsed laser (2 ns pulse). The time response value is obtained through single exponential fitting.

Figure 4 :

 4 Figure 4 : a) Photocurrent spectra of the diode under different biases. b) Measured ratio of photocurrent at 4000 cm -1 by photocurrent at 6000 cm -1. as a function of applied bias, inserts: band diagrams at different biases. The black line is the conduction band, the red line is the valence band. The black and red arrows correspond to the electron and hole trajectories respectively.

Figure 5 :

 5 Figure 5: a) Set up schematic for long distance detection. The measured signal is the 100 Hz current contribution modulated at 1 Hz. b) Signal current measured at different black body (BB) temperatures. c) Photograph of the experimental set up. d) Angular detection with a 1.55 µm, 4 mW laser. The black dots correspond to the experimental ratio between the responsivity at a given angle and at 0° (i.e. normal incidence), the red circle corresponds to a cosine evolution.

Figure 6 a

 6 Figure 6 a) Schematic of the matrix at different zoom levels. b) SEM image of the matrix (magnified 200x). c) SEM image of the matrix with faux colors corresponding to different layers: yellow is the internal gold contact (columns), green is the insulating silica and red is the external gold contact (lines) (magnified 3000x). Blue relates to the glass substrate.

Figure 7 :

 7 Figure 7: Matrix response. a) Dark current under 0.1 V for all the pixels on the matrix except the three shorted. b) Determination of 1.55 µm focused laser spot size using a razor blade. The experimental points are the black dots, the Gaussian fit is the red line. c) Intensity of the 97 working pixels with laser spot focused on the matrix. The Gaussian shape has a full width at half maximum (FWMH) of 75 µm.

ACKNOWLEDEGMENTS

EL thanks the support ERC starting grant blackQD (grant n°756225). We acknowledge the use of clean-room facilities from the "Centrale de Proximité Paris-Centre". This work has been supported by the Region Ile-de-France in the framework of DIM Nano-K (grant dopQD). This work was supported by French state funds managed by the ANR within the Investissements d'Avenir program under reference ANR-11-IDEX-0004-02, and more specifically within the framework of the Cluster of Excellence MATISSE and also by the grant IPER-nano2. NG and JR thank Nexdot for post doc funding.

Supporting information

Additional information about nanocrystal synthesis, material characterization, diode fabrication and characterization as well as focal plane array fabrication and characterization are given in supporting information.

Supporting Information is available from the Wiley Online Library.


Gaspard Farge
email: farge@ipgp.fr

Nikolaï M Shapiro

William B Frank

Moment-duration scaling of Low-Frequency Earthquakes in Guerrero, Mexico

Low-frequency earthquakes (LFEs) are detected within tremor, as small, repetitive, impulsive low-frequency (1-8 Hz) signals. While the mechanism causing this depletion of the high-frequency content of their signal is still debated, this feature may indicate that the source processes at the origin of LFEs are different from those for regular earthquakes. Key constraints on the LFE-generating physical mechanisms can be obtained by establishing scaling laws between their seismic moment and source durations. Here we apply a simple spectral analysis method to the S-waveforms of LFEs from Guerrero, Mexico to measure their seismic moments and corner frequencies, a proxy to source duration. We find characteristic values of M0 ∼ 3 × 10 12 N.m (Mw ∼ 2.3) and fc ∼ 3.0 Hz with the corner frequency very weakly dependent on the seismic moment. This moment-duration scaling observed for Mexican LFE is similar to one previously reported in Cascadia and is very different from the established one for regular earthquakes. This suggests that they could be generated by sources of nearly constant size with strongly varying intensities. LFEs do not exhibit the self-similarity characteristic of regular earthquakes, suggesting that the physical mechanisms at their origin could be intrinsically different.

Plain language summary

Low-frequency earthquakes are unusual, small earthquakes that are detected on the deep end of plate boundary faults. They occur during periods of slow slip, when the plates start sliding against each other much slower than during earthquakes, but long enough to release the energy that might otherwise have caused large earthquakes. Thus, they carry precious information about the fault behavior while slow slip is active.

In this study, we analyze a selection of low-frequency earthquakes from a slow slip zone of the Mexican subduction, beneath the state of Guerrero. We find that their magnitude is on average higher than in other regions where they have been detected (Mw1.5 --3). We also find that whatever the magnitude of these events, they always last about 0.3 s.

On the contrary, regular earthquakes are known to last longer when they grow larger, and lowfrequency earthquakes in Nankai (Japan) appear to behave similarly. We discuss a physical mechanism explaining such a discrepancy with regular earthquakes and such strong regional variations. The implication of crustal fluids circulating at extremely high pressure in between plates might be one of the key components of this mechanism, as it often is for slow slip processes.

Introduction

Low-frequency earthquakes (LFEs) are observed in association with volcanic and tectonic processes as impulsive, coherent wave arrivals with discernible S and sometimes P phases and a low-frequency content (1-8 Hz) relative to their small magnitudes. Tectonic LFEs have been documented in various subduction zones [e.g 7, [START_REF] Frank | Low-frequency earthquakes in the mexican sweet spot[END_REF][START_REF] Shelly | Low-frequency earthquakes in shikoku, japan, and their relationship to episodic tremor and slip[END_REF] and on strike-slip plate boundaries [e. g 12, 62]. They are detected in the transitional zones just beneath the regular seismogenic portions of the faults and are often associated with zones of relatively low shear-wave velocity, high V P /V S ratio, low quality factor, all indicators of fluid rich media [START_REF] Audet | Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing[END_REF][START_REF] Bostock | Low frequency earthquakes below southern vancouver island[END_REF][START_REF] Shelly | Low-frequency earthquakes in shikoku, japan, and their relationship to episodic tremor and slip[END_REF]. LFE activity is closely correlated in time and space with slow slip activity, and is thus used to track the spatial extent of slow slip phenomena [START_REF] Ghosh | Rapid, continuous streaking of tremor in cascadia[END_REF][START_REF] Houston | Rapid tremor reversals in cascadia generated by a weakened plate interface[END_REF][START_REF] Poiata | Imaging different components of a tectonic tremor sequence in southwestern japan using an automatic statistical detection and location method[END_REF][START_REF] Shelly | Complex evolution of transient slip derived from precise tremor locations in western shikoku, japan[END_REF] but also to detect short-duration and low-amplitude slow slip events that fall below the detection threshold for techniques relying on geodetic data [START_REF] Frank | Slow slip hidden in the noise: The intermittence of tectonic release[END_REF][START_REF] William B Frank | Revealing the cluster of slow transients behind a large slow slip event[END_REF][START_REF] Lengliné | Imbricated slip rate processes during slow slip transients imaged by low-frequency earthquakes[END_REF].

In analogy with volcanic tremor, it was first suggested that tectonic tremors are generated by fluids circulating in the vicinity of the fault zone [START_REF] Kao | A wide depth distribution of seismic tremors along the northern cascadia margin[END_REF][START_REF] Katsumata | Low-frequency continuous tremor around the moho discontinuity away from volcanoes in the southwest japan[END_REF][START_REF] Obara | Nonvolcanic deep tremor associated with subduction in southwest japan[END_REF]. Subsequent research based on the characteristics of tectonic tremors and LFEs signals -such as dominant S-wave content and shear double-couple mechanisms [START_REF] Bostock | Low frequency earthquakes below southern vancouver island[END_REF][START_REF] Frank | Low-frequency earthquakes in the mexican sweet spot[END_REF][START_REF] Ide | Mechanism of deep low frequency earthquakes: Further evidence that deep non-volcanic tremor is generated by shear slip on the plate interface[END_REF][START_REF] Shelly | Non-volcanic tremor and low-frequency earthquake swarms[END_REF] and association with slow slip events [START_REF] William B Frank | Uncovering the geodetic signature of silent slip through repeating earthquakes[END_REF][START_REF] Kostoglodov | The 2006 slow slip event and nonvolcanic tremor in the mexican subduction zone[END_REF][START_REF] Rogers | Episodic tremor and slip on the cascadia subduction zone: The chatter of silent slip[END_REF] -built empirical support for the hypothesis that they are generated by shear slip on the fault.

LFEs are a member of the "slow earthquake" class of phenomena, regrouping events which seem to be associated with slow slip processes on faults [START_REF] Beroza | Slow earthquakes and nonvolcanic tremor[END_REF]. Those events exhibit a wide range of durations (from slightly less than 1 s for LFEs, up to nearly a year for slow slip events) and magnitudes (from around M w 1 up to M w 7.5). When accounted as a one family of fault slip processes, their seismic moment M 0 appears to be directly proportional to their duration T [START_REF] Ide | A scaling law for slow earthquakes[END_REF], whereas regular earthquakes align along a M 0 ∝ T 3 scaling law [START_REF] Allmann | Global variations of stress drop for moderate to large earthquakes[END_REF][START_REF] Houston | Influence of depth, focal mechanism, and tectonic setting on the shape and duration of earthquake source time functions[END_REF][START_REF] Prieto | Earthquake source scaling and self-similarity estimation from stacking p and s spectra[END_REF] due to the self-similarity of the rupture process [START_REF] Kanamori | Theoretical basis of some empirical relations in seismology[END_REF]. However, considering slowearthquakes as one entity relies on bridging two different scales of observations: geodetic observations of large-scale, long-term (sub-daily at minima) slow slip events and seismic observations of short term events as LFEs, and longer events (30-60 s) called very-low-frequency earthquakes (VLFEs). There seems to be evidence that the shortest slow slip transients measured to this day align along a self-similar momentduration scaling M 0 ∝ T 3 [START_REF] William | Daily measurement of slow slip from low-frequency earthquakes is consistent with ordinary earthquake scaling[END_REF], and that only at the largest scales it would transition to a M 0 ∝ T scaling law [START_REF] Gomberg | Reconsidering earthquake scaling[END_REF]. This observation isolates LFEs and VLFEs from the bulk of geodetic observations of slow slip. It brings into question whether they are simply smaller scale manifestations of slow slip, or they have a more complex coupling relationship with it. In particular, LFEs co-occur with slow slip events, and investigating their properties provides key insight to understanding fault state, geometry of activation and overall rupture process during slow slip activity [5, 13, 21, 30, e.g.].

In this study, we investigate the moment-duration relationship for LFEs in Guerrero, Mexico. As for regular earthquake, evaluating how moment release scales with source duration is key to understand the dynamics of the process generating LFEs. Similar studies have already been carried out in regions witnessing tremor and slow slip. In Cascadia, LFEs source duration has been found to be very weakly dependent on their seismic moment, scaling along M 0 ∝ T 10 [8]. On the other hand, in Nankai, LFEs exhibit a self-similar behavior, scaling along M 0 ∝ T 3 [START_REF] Supino | Self-similarity of low-frequency earthquakes[END_REF]. To our knowledge this is the first published work to measure source parameters for LFEs in Guerrero, Mexico. We start with a large LFE catalog compiled by [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in guerrero, mexico[END_REF] and select those whose waveforms have a sufficiently high quality to measure their seismic moments and corner frequencies. The latter are estimated with two different methods to check the robustness of the obtained results. Similar to the study of [START_REF] Bostock | Magnitudes and momentduration scaling of low-frequency earthquakes beneath southern vancouver island[END_REF] for LFEs in Cascadia, we find that for the LFEs in Guerrero the corner frequency is very weakly dependent on the seismic moment.

Low-frequency earthquakes in Guerrero, Mexico

Slow-slip events (SSE), tectonic tremor, and low-frequency earthquakes have been observed in the state of Guerrero, Mexico, in the subduction zone of the Cocos plate under the North American plate [START_REF] Frank | Low-frequency earthquakes in the mexican sweet spot[END_REF][START_REF] Kostoglodov | The 2006 slow slip event and nonvolcanic tremor in the mexican subduction zone[END_REF][START_REF] Payero | Nonvolcanic tremor observed in the mexican subduction zone[END_REF]. They have been detected in the Guerrero seismic gap -a corridor extending from Acapulco to Mexico City -that has not witnessed major earthquakes (M w > 7) for more than 100 years [START_REF] Kostoglodov | Cien años de sismicidad en méxico, instituto de geofísica, universidad nacional autónoma de méxico[END_REF]. This longer recurrence time of major earthquakes compared to adjacent regions can be related to the large SSEs that accommodate a significant fraction of the plates convergence in the Guerrero seismic gap [START_REF] Radiguet | Slow slip events and strain accumulation in the guerrero gap, mexico[END_REF].

We use the data collected during the Meso-American Seismic Experiment (MASE) [START_REF] Husker | Wilson: The wirelessly linked seismological network and its application in the middle american subduction experiment[END_REF][START_REF] Perez-Campos | Horizontal subduction and truncation of the cocos plate beneath central mexico[END_REF], between 1 January 2005 and 15 April 2007. Following the work of [START_REF] Frank | Automatic detection of low-frequency earthquakes (lfes) based on a beamformed network response[END_REF], we selected a subset of 10 stations out of the 100 total based on their high signal-to-noise ratio (SNR), the continuity in time of their records, and their proximity to the LFE activity (see Figure 1). [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in guerrero, mexico[END_REF] used this dataset and the method described in [START_REF] Frank | Automatic detection of low-frequency earthquakes (lfes) based on a beamformed network response[END_REF] to compile a large catalog of 1.8 million LFEs grouped in 1120 families of repeating events. We [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in guerrero, mexico[END_REF] is in gray. Subduction interface depth contours are displayed in black [START_REF] Kim | Geometry and seismic properties of the subducting cocos plate in central mexico[END_REF]. The focal mechanism corresponds to an average mechanism computed from cataloged LFE detections [START_REF] Frank | Low-frequency earthquakes in the mexican sweet spot[END_REF]. (c.) shows the time-dip distribution of the events. use this catalog as a starting point in our study. It should be noted that the events in the catalog are detected by a matched-filter search using template events, filtered between 1 and 2 Hz. This band-pass is narrower than the range most studies of LFEs using matched-filter search have used up to now, usually 1-8 Hz [8, 12, e.g] Nonetheless, we show in appendix A.1 that performing the matched-filter search in a narrow band does not bias the event detection and measurements of events' seismic moments and corner frequencies. Within the range of measured seismic moments, events with a corner frequency higher than 1 Hz do not have a lower probability of being detected because of the band-limited nature of detection.

Selection of LFEs with high-quality waveforms

The catalog of [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in guerrero, mexico[END_REF] has been created on the network-based template matching detection method [27, for instance]. The advantage of this method is that it detects many events with very weak signals hidden in the noise. At the same time, the signal-to-noise ratio for most LFEs in the final catalog is too low to be suited for a determination of their source parameters. Therefore, we selected LFEs with high-quality waveforms that could be used to robustly measure their seismic moments and corner frequencies. Some single-component waveforms of catalog detections are contaminated with very high amplitude spikes, instrumental noise or waves from strong earthquakes. In these cases, useful signals are masked and corresponding waveforms cannot be used for determination of the LFE source parameters. We remove those seismograms from the analysis. Determination of seismic moment requires knowing the earthquake source location. However, LFE family locations determined by [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in guerrero, mexico[END_REF] contain a strong uncertainty on position in the direction perpendicular to the station array. This uncertainty is inherent to the linear observation geometry of the MASE network. Consequently, a few families of LFEs have templates with diffuse seismic phases, indicating that their measured position is loosely constrained. We therefore remove those families of events from our analysis.

We establish a set of criteria to automatically select LFEs and waveforms suitable for measuring the source parameters. We select 5 s long signals starting 2 s before the cataloged detection time of LFE S-waves.

This allows us to capture the most impulsive part of the S-wave signal, and allows for slight location and detection time imprecisions. This is done so as to trade-off between keeping a maximum of information and not allowing for too much diffuse seismic energy in the form of other seismic phases and events. By isolating the most impulsive part of the wave train, we best constrain the characteristic frequencies and seismic moment of the event.

In a next step, we reject signals with low amplitude and impulsivity to exclude false detections and contaminations with surface waves. First we define two frequency ranges in which we determine both an impulsivity I and an amplitude A: they will allow the criterion to be more flexible as for which events are taken into consideration. The impulsivity is calculated as the ratio of the seismic energy in the waveform time window W to the seismic energy contained in a three times larger background time window B around the detection, containing it:

I = W s(t) 2 dt B s(t) 2 dt (1)
where s(t) is the velocity signal. I varies between 0, for least impulsive, and 1, for most impulsive. It is used as a proxy of the detection SNR, but it also addresses the difficulty of defining what is noise within the tremor, where the LFE rate is very high. The amplitude is calculated as the maximum of the envelope of the filtered signal within the 5 s time window. For this purpose, the waveforms are converted into complex analytical signals s A (t):

s A (t) = E(t) . e i2πf i (t) A = max W (E(t)) (2)
E(t) is the signal envelope and f i (t) is its instantaneous frequency. The latter will be used later on to determine the corner frequency of the selected LFEs. Amplitudes and impulsivities are computed for all events on each channel of each station in the 1-2 Hz frequency band. This frequency range is characteristic of Mexican tremor signal in velocity, and has been shown to yield the best SNR for LFEs in this region [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in guerrero, mexico[END_REF][START_REF] Payero | Nonvolcanic tremor observed in the mexican subduction zone[END_REF].

So as to robustly ensure that the detection is a correctly located LFE we also base the selection criterion on the correlation coefficient used in the making of the LFE catalog. It is a measure of the correlation of an event waveforms on all 15 detection channels to the stacked-waveforms template of its family [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in guerrero, mexico[END_REF]. We consider the highest absolute values of the correlation coefficient to represent a wellconstrained location and an impulsive event. As a sum of 15 normalized correlation coefficients, it is normalized by the number of channels used in the detection of the event and thus varies between 0 and 1.

The selected waveforms have to fulfill the following criterion, where I is the waveform impulsivity and A its amplitude:

• The LFE should be impulsive, but waveforms with very high impulsivity are instrumental noise or earthquake wave arrivals: 0.5 < I < 0.92

• Higher-amplitude LFEs are generally detected in distant earthquakes waveforms. Therefore, we empirically limit the amplitude of a selected detection to 2 × 10 -6 m.s -1 . A lower threshold for selection T low (st, ch) is fixed for each station and component (approximately 5 × 10 -9 m.s -1 ). Both bounding values are manually determined using the impulsivity-amplitude distribution on each channel of each station (see Figure 2). For our source analysis, we select LFEs that contain at least 10 different channels satisfying the described criteria and have a network correlation coefficient superior to 0.4. This value has been chosen empirically, as the median of the correlation coefficients in the catalog.

T low (st, ch) < A < 2 × 10 -6 m.s -1 (4)
A benefit of the network template matching method is that it allows to detect earthquakes with very close origin times, and thus with overlapping signals at some stations. Overlapping signals are not suitable for the source parameter analysis. Therefore, if several LFEs are detected within a 10 s time window, only the LFE with the highest correlation coefficient, and thus most reliable waveforms, is kept. Figure 2 sums up the selection process.

The selection yields 3498 LFEs representing 822 families out of a total of 1120. Each event is recorded on 10 selected channels over the network and shown with blue circles in Figure 1. The most impulsive waveforms of two selected LFEs are displayed in Figure 3. The selected LFEs exhibit a similar activity behavior in space and time as [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in guerrero, mexico[END_REF] observed for the whole catalog. First, they occur mainly in bursts within the sweet spot region. Second, during the 2006 SSE, bursts are no longer limited to the sweet spot, but extend updip in the transient zone. The activity in those two regions resembles the description [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in guerrero, mexico[END_REF], the sweet spot has a nearly constant activity over the 2-year span, whereas the transient zone is active mainly during the 2006 SSE and more sparse in the inter-SSE period.

Characterizing seismic moment and corner frequency 4.1 Seismic moment

We estimate seismic moment M 0 from the displacement spectra of S-waves u(f ) that converges to a constant value at low frequencies [START_REF] Brune | Tectonic stress and the spectra of seismic shear waves from earthquakes[END_REF]. This constant is proportional to M 0 [1]. For every selected Vel. (m/s) The blue section of the waveform shows the time window used to define the detected signal, the red patch shows which channel's waveforms are used as an example for source parameters characterization in Figure 4. Impulsivity and amplitude are indicated for each channel.

waveform, we can determine the low-frequency spectral asymptotic value Ω 0 i,j as:

Ω 0 i,j = lim f →0 u i,j (f ) = C i,j . M 0 i,j , (5)
where i corresponds to the LFE number and j is an index given to each of the corresponding selected waveforms. The C i,j factor accounts for the source radiation pattern and the seismic wave propagation.

To measure these low-frequency asymptotic values, the 5 s long S-wave velocity waveforms selected in the previous section are first detrended and high-pass filtered above 1 Hz, to remove any influence of spectral leakage of microseismic noise, and then integrated into displacement. The power spectral density (PSD) of the displacement waveforms is computed with the multi-taper method [START_REF] Thomson | Spectrum estimation and harmonic analysis[END_REF]. It gives more robust spectral estimations than a simple discrete Fourier transform, especially for short signals.

Taking the square root of this spectral estimation yields the displacement spectral amplitude. The lowfrequency value of the displacement spectral amplitude Ω 0 i,j is estimated as the geometric mean of the 1-2 Hz portion of the spectral amplitude of displacement, filtered between 1 and 8 Hz (Equation 5, and Figure 4).

Then, seismic moment estimates from every waveform M 0 i,j can be retrieved by computing the propagation effects factor C i,j . For this purpose, we compute synthetic waveforms for each sourcestation couple, using LFE family locations from [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in guerrero, mexico[END_REF] as hypocenters, an average regional seismic velocity model [START_REF] Iglesias | S wave velocity structure below central mexico using high-resolution surface wave tomography[END_REF] and an average LFE focal mechanism, determined by [START_REF] Frank | Low-frequency earthquakes in the mexican sweet spot[END_REF]. We used the axitra software [START_REF] Coutant | axitra software[END_REF] based on the discrete wavenumber method [START_REF] Bouchon | A review of the discrete wavenumber method[END_REF].

After the synthetic seismograms are computed, we select a 5 s long time window around the S-waves, ) Ω 0 is determined as the red plateau of the displacement spectrum. (b.) s Ω 0 is determined as the maximum of the synthetic displacement spectrum, taking into account attenuation. (c.) The corner frequency f c is measured as the mean instantaneous frequency (MIF) over the waveform time window (see Figure 3 for the LFE waveforms), (d.) and as the maximum of the velocity spectrum (MVS) of the waveform. In gray, the un-filtered velocity spectrum, allowing to assess the effects of the 1-8 Hz filter on the spectrum.

process them in the same way as the LFEs' seismogram and compute their displacement spectrum without attenuation s 0 u(f ), where f is the frequency. The seismic wave attenuation effects are modeled in the spectral domain by introducing a(f ). A corrected synthetic spectral displacement s u(f ) becomes:

s u(f ) = a(f ) . s 0 u(f ) a i,j (f ) = exp( -2πf X i,j 2V S Q(f ) ), where Q(f ) = 273 f 0.66 (6)
where the frequency-dependent quality factor Q(f ) for the Guerrero region has determined by [START_REF] Ordaz | Source spectra and spectral attenuation of seismic waves from mexican earthquakes, and evidence of amplification in the hill zone of mexico city[END_REF], the average shear wave velocity V S is chosen to be 3.5 km.s -1 and X i,j is the source-receiver distance. We determine the low-frequency asymptote of the synthetic spectrum s Ω 0 i,j as the maximum of the synthetic spectrum, filtered between 1 and 8 Hz. Knowing the synthetic source moment s M 0 and following equation ( 5), C i,j can be eliminated:

Ω 0 i,j s Ω 0 i,j = M 0 i,j s M 0 , (7)
and a seismic moment estimated for each waveform of a selected LFE becomes:

M 0 i,j = s M 0 Ω 0 i,j s Ω 0 i,j . (8)

Corner frequency

We use two different techniques to estimate the corner frequency from all selected waveforms for every LFE, f c i,j . For both methods, the waveforms are detrended, corrected for attenuation (see equation 6) and band-pass filtered between 1 and 8 Hz to remove any influence of micro-seismic signal and higher frequency noise from the frequency band that correspond to LFEs frequencies. Further discussion on the choice of this frequency band can be found in the appendix, section A.2.

The first approach consists in automatically picking f c i,j as the frequency corresponding to the maximum in the velocity amplitude spectrum of every waveform. This method is labeled MVS (maximum of velocity spectrum) in the illustrations. The bandpass filter creates a slight amplitude reduction up to 1.5 Hz, which can be seen when comparing the filtered and unfiltered spectra in Figure 4. This could shape the spectrum into a peak slightly above 1 Hz, if the actual peak frequency between 1 and 8 Hz is close to 1 Hz. However, we find that the 1-8 Hz filter band is the best compromise between allowing microseismic noise to bias the corner frequency measurement and slightly constraining the shape of the spectrum. The effect on measured corner frequencies of the low-cutoff of the filter is discussed more extensively in section A.2 and can be assessed in figure S2 of the supporting information file.

In the second approach, we estimate the corner frequency as the dominant instantaneous frequency of the S-waves. The waveforms are converted into complex analytical signals (see equation 2) and f c i,j is estimated as the mean of the smoothed instantaneous frequency over the 5 s time window containing the S-waves. This method is labeled MIF (mean of instantaneous frequency) in the illustrations. For this method too, the effects of the filter band can be assessed in figure S2 of the supporting information file.

For each LFE, we obtain a final estimate of their moment M 0 i and two estimates of their corner frequency f c i , all computed as medians of those measured from individual waveforms selected for this LFE. For each event, the uncertainty on its source parameters is estimated as the standard deviation of the values of the given source parameter measured on each record of the event.

Results

Low-frequency earthquakes source parameters: seismic moment, corner frequency

Seismic moments and corner frequencies estimated for the selected LFEs are shown in Figure 5 along with the schematic imprint of the M 0 -f c distribution for classical earthquake, which follows f c ∝ M 0 -1 3 , for stress drops ranging from 0.1 to 100 MPa [START_REF] Allmann | Global variations of stress drop for moderate to large earthquakes[END_REF], and values of M 0 and f c determined for LFEs in previous studies. Seismic moments and magnitudes estimated for Mexican LFEs range between M 0 = 4.9 × 10 11 N.m (M w ∼ 1.7) and M 0 = 2.3×10 14 N.m (M w ∼ 3.5). These values are higher than previous estimations in other regions [START_REF] Bostock | Magnitudes and momentduration scaling of low-frequency earthquakes beneath southern vancouver island[END_REF][START_REF] Ide | A scaling law for slow earthquakes[END_REF]. The median uncertainty on moment magnitude of δM w ≈ 0.23. It is mainly due to errors in modeling the velocity structure of the crust, assuming a generic radiation pattern for all events and large uncertainties on the hypocentral locations [20, supplementary figures 3 to 17]. The distribution of seismic moments above a threshold magnitude of M wthr ∼ 2.32 follows a power-law, with an exponent β = 1.62 (corresponding b-value b = 2.43, see Appendix, Figure 7).

The measured corner frequencies mean and extrema for each method are displayed in Table 1. They are approximately centered around f c ∼ 3.0 Hz and fall within the accepted range of characteristic frequencies for LFEs of 1-6 Hz [START_REF] Bostock | Magnitudes and momentduration scaling of low-frequency earthquakes beneath southern vancouver island[END_REF][START_REF] Ide | A scaling law for slow earthquakes[END_REF][START_REF] Shelly | Non-volcanic tremor and low-frequency earthquake swarms[END_REF][START_REF] Thomas | Constraints on the source parameters of low-frequency earthquakes on the san andreas fault[END_REF]. The median uncertainty on corner frequency is δf c ≈ 1.6 for the maximum of velocity spectrum method, and δf c ≈ 0.8 for the mean instantaneous frequency method.

Moment-duration scaling

The M 0 -f c distribution in Figures 5 and6 shows that the corner frequencies seem virtually independent of seismic moment. In order to estimate the scaling law between the source parameters, M 0 -f c points are b.

fit, slope=-0.12 other studies [66], 3. [32], 4. [START_REF] Bostock | Magnitudes and momentduration scaling of low-frequency earthquakes beneath southern vancouver island[END_REF]) The self-similarity scaling law for classic earthquakes is represented by the purple patch. The higher bounding line corresponds to a stress drop of 100 MPa and the lower to a stress drop of 0.1 MPa [START_REF] Allmann | Global variations of stress drop for moderate to large earthquakes[END_REF][START_REF] Douglas | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF][START_REF] Madariaga | Dynamics of an expanding circular fault[END_REF].

f c ∝ M -1/3 0 current study
binned into 35 bins of ascending moment magnitude, containing 100 points each. Moment magnitude of the bin is considered to be the mean M w of the events in the bin. Corner frequency of the bin f c bin is the weighted mean of the measured corner frequency of the binned events. The weight associated to the corner frequency of event k is chosen to be the inverse of the variance of corner frequencies measured on each record for event k, w k = 1/σ k 2 . Figure 6 shows the bin distribution and the associated estimation of variability of corner frequency within bin, calculated as a weighted standard deviation:

σ bin = k (f c bin -f c k ) 2 × w k k w k (9)
with f c k the measured corner frequency for event of index k in the considered bin, w k = 1/σ k 2 its associated weight (see above), and f c bin as the bin's corner frequency.

Using the bins M 0 and f c and their associated weight w bin = 1/σ bin 2 , a weighted least-square fit to the following power law is performed, in log-log space:

f c ∝ M 0 -α ⇐⇒ log(f c ) = A -αlog(M 0 ). ( 10
)
The resulting values of the slopes confirm that corner frequencies are weakly dependent on seismic moment : for the mean instantaneous frequency method (MIF) used to obtain f c : α = 0.052, for the 10 12 10 13 10 14

Seismic moment M 0 (N.m) maximum of velocity spectrum method (MVS): α = 0.115, (see Figures 5 and6). This result is similar to the moment-duration scaling observed for LFEs in the Cascadia subduction zone by [START_REF] Bostock | Magnitudes and momentduration scaling of low-frequency earthquakes beneath southern vancouver island[END_REF], where α ∼ 1 10 . We proceed to estimate the probability density function (PDF) of the exponent of the scaling using a bootstrapping method. For each bin, we keep its moment but pick a new corner frequency from a normal distribution parametrized by the bin's corner frequency µ = f c bin and by the bin's weighted standard deviation σ = σ bin (see equation 9). An unweighted, least-square fit to a power law (equation 10) is then performed on the bin's seismic moment and re-sampled corner frequency. After 200,000 iterations of the previous steps, we obtain an estimate of the probability distribution function for the scaling parameter α. The mean and standard variation of the obtained distribution give a new estimation of the slope and its associated uncertainty. For the MIF method, α = 0.054±0.034, for the MVS method, α = 0.12±0.11. Assuming that α is normally distributed, the probability that it is higher than 0.25 (thus closer to 1 3 ) for our f c -M 0 distribution is around 0.12 for the MVS method (α = 1 3 is ∼ 1.2σ away from the mean α), and lower than 10 -8 for the MIF method (α = 1 3 is ∼ 5.7σ away from the mean α). The scaling exponent α we measure is thus significantly lower than the earthquake power law exponent α = 1 3 , for the MIF method at least. The MVS method gives more ambiguous results with regard to the value of the exponent α, due to the wider distribution on f c . This work is published in Journal of Geophysical Research: Solid Earth, DOI:10.1029/2019JB019099.

Stress-drop variations

The stress-drop for an event of seismic moment M 0 modeled by a circular dislocation of radius r is expressed as follows [START_REF] Madariaga | Dynamics of an expanding circular fault[END_REF]:

∆σ = 7 16 M 0 r 3 . ( 11
)
Assuming a rupture velocity proportional to the shear wave velocity, v r ∝ V S , it follows that the source radius is inversely proportional to the corner frequency [START_REF] Madariaga | Dynamics of an expanding circular fault[END_REF]:

r ∝ V S f c . ( 12
)
From equations 11 and 12, and keeping the assumptions that our events originate from regions of constant V S and share the average corner frequency of our distribution f c ∼ 3 Hz, we can express the stress-drop ratio for two events of seismic moment M 01 and M 02 :

∆σ ∝ M 0 =⇒ ∆σ 1 ∆σ 2 = M 01 M 02 (13)
Using the extrema of the seismic moment distribution, and within the limit of our assumptions, we observe a stress-drop variation of a factor 300 across the range of measured seismic moments. For illustration purposes, we calculate what values would the stress-drop of our events take under the assumption that they have the same rupture speed as classical earthquakes, v r = 0.9V S . It yields r = 0.21 V S fc [START_REF] Madariaga | Dynamics of an expanding circular fault[END_REF], and using equation 11, the selected events have stress-drop ranging from 3 kPa to 1 MPa, with most events around 17 kPa (M 0 ∼ 2 × 10 12 N.m). In the following section, we discuss a possible mechanism for the observed range of stress-drops.

6 Discussion and interpretation of the observed scaling law

Scaling laws of earthquakes, slow earthquakes and LFEs

The results of our analysis show that the selected LFEs follow a moment-duration scaling f c ∝ M 0 -α with an exponent significantly lower than the value for regular earthquakes α = 1/3 [2, 29, 54, among others], around α ∼ 0.1 and α ∼ 0.05. The scaling found for selected LFEs in Guerrero concurs with the one found for LFEs in Cascadia by [START_REF] Bostock | Magnitudes and momentduration scaling of low-frequency earthquakes beneath southern vancouver island[END_REF], who found a power law exponent of α = 1 10 , quite similar to the scaling described in this paper. It is noteworthy that the detection method used by [START_REF] Bostock | Magnitudes and momentduration scaling of low-frequency earthquakes beneath southern vancouver island[END_REF] relies on a matched-filter search, similarly to the catalog used in our study. However, we argue that performing a matched-filter search in the 1-2 Hz band sets a lower bound, but no upper bound for the characteristic frequency, for a given seismic moment. Therefore, the values of corner-frequencies should not be the product of the detection method, but reflect the real moment-corner frequency distribution of events, within the measured magnitude range. This argument is developed in appendix A.1. Interestingly enough, [START_REF] Supino | Self-similarity of low-frequency earthquakes[END_REF] found the moment-duration scaling for LFEs in Nankai tremor areas to be similar to regular earthquakes, α = 1 3 . This discrepancy will not be thoroughly discussed here, but brings out intriguing questions, about the possibly different nature of the processes recorded in Guerrero, Cascadia and Nankai, and on the other hand, the potential detection and method biases that can lead to measuring different scalings.

Using a stochastic model describing a slow earthquake rupture (the Brownian slow earthquake model), [START_REF] Ide | A brownian walk model for slow earthquakes[END_REF][START_REF] Ide | Seismic moment, seismic energy, and source duration of slow earthquakes: Application of brownian slow earthquake model to three major subduction zones[END_REF] showed that large-scale slow earthquakes as slow slip events' duration and seismic moment should follow a T ∝ M 0 scaling, whereas on the scale of LFEs, slow earthquake events should follow a T 2 ∝ M 0 scaling, corresponding a value of α = 1/2. In [START_REF] Ide | A brownian walk model for slow earthquakes[END_REF], it is noted that, due to the fact that detection is only possible above noise levels, the apparent duration of LFEs should be shorter than expected, bringing them closer to a moment-duration scaling with α = 1. The LFEs we analyzed in Guerrero and the LFEs of Cascadia [START_REF] Bostock | Magnitudes and momentduration scaling of low-frequency earthquakes beneath southern vancouver island[END_REF] cannot be described by either scaling laws.

There is evidence that the shortest slow slip transients measured to this day are self-similar, exhibiting a T 3 ∝ M 0 scaling law [START_REF] William | Daily measurement of slow slip from low-frequency earthquakes is consistent with ordinary earthquake scaling[END_REF]. The exact proportionality measured for larger scales slow earthquakes [START_REF] Ide | A scaling law for slow earthquakes[END_REF] could be due to the fact that the largest ruptures are bounded in their growth by the rheological properties of the fault zone, whereas smaller events growth is un-bounded, and thus self-similar [START_REF] Gomberg | Reconsidering earthquake scaling[END_REF]. This geometric argument is valid both for slow and regular earthquakes, and further isolates LFEs momentduration characteristic values and scaling from the ones of slow and regular slip. For LFEs in Cascadia and Guerrero, a satisfactory physical description that could reproduce the observation of the constant duration and relatively low-frequency content still needs to be devised.

Potential physical mechanisms accounting for magnitude-independent duration of LFEs

A source of constant size, but variable slip

The scaling with α = 1 3 for regular earthquake is explained by the self-similarity of the seismic rupture. It develops with a constant ratio between the characteristic length of the rupture and fault slip, and produces a magnitude-independent stress drop [START_REF] Allmann | Global variations of stress drop for moderate to large earthquakes[END_REF][START_REF] Kanamori | Theoretical basis of some empirical relations in seismology[END_REF][START_REF] Prieto | Earthquake source scaling and self-similarity estimation from stacking p and s spectra[END_REF]. Events with a wide range of magnitudes but with nearly constant durations could thus be generated by variable amounts of slip on asperities of nearly constant sizes, producing a very limited range of recorded signal duration (f c ) but variable seismic moments [START_REF] Bostock | Magnitudes and momentduration scaling of low-frequency earthquakes beneath southern vancouver island[END_REF][START_REF] Lengliné | Fluid-induced earthquakes with variable stress drop[END_REF][START_REF] Lin | Evidence for non-self-similarity of microearthquakes recorded at a taiwan borehole seismometer array[END_REF]. Independently from the current study, [START_REF] William | Daily measurement of slow slip from low-frequency earthquakes is consistent with ordinary earthquake scaling[END_REF] found that, in Guerrero, the seismic moment rate, calculated from LFEs median amplitude during rapid slow slip transients scales with the moment rate calculated using GPS displacements during the same transient events. This reflects that the higher the seismic moment of LFEs, the higher the slip would be around the LFEs source. Combining those two observations leads to thinking of LFE sources in Guerrero as patches of fault of nearly constant area, that are forced to rupture by the surrounding slow, aseismic slip, with variable amplitudes of slip directly related to the surrounding fault slip. A similar model has been designed to explain the observation of a size-limitation for LFEs in Cascadia, describing the LFE family patch as a collection of subpatches randomly activated, totaling to an apparent source radius of about 300 m [START_REF] Sr Chestler | Evidence for a scale-limited low-frequency earthquake source process[END_REF].

Potential effects of fluid circulation and high fluid pressures

Early work on tectonic tremor and LFEs made the hypothesis that they are generated by fluids circulating in the vicinity of the fault zone [START_REF] Kao | A wide depth distribution of seismic tremors along the northern cascadia margin[END_REF][START_REF] Katsumata | Low-frequency continuous tremor around the moho discontinuity away from volcanoes in the southwest japan[END_REF][START_REF] Obara | Nonvolcanic deep tremor associated with subduction in southwest japan[END_REF]. More recent studies have strengthened the link between nonstationnary, metamorphic fluid flow and migrations of LFEs and slow slip activity, based on numerical models of pore-pressure diffusion [START_REF] Víctor | Rapid tremor migration and porepressure waves in subduction zones[END_REF] and on analysis of the large-scale, spatio-temporal behavior of LFE activity [START_REF] Frank | Alongfault pore-pressure evolution during a slow-slip event in guerrero, mexico[END_REF].

A study by [START_REF] Lengliné | Fluid-induced earthquakes with variable stress drop[END_REF] reported that microseismic events occurring during water circulation tests display the same independence of corner frequency on seismic moment that we observe for our LFEs. As for the LFEs in this study, they find that stress drop varies on more than two orders of magnitude across the range of moments of their events. Likewise, [START_REF] Lin | Evidence for non-self-similarity of microearthquakes recorded at a taiwan borehole seismometer array[END_REF] reports that microseismic events detected in Taiwan exhibit the same characteristics. Both studies suggest that these observations can be interpreted as a systematic control of the events' stress-drop and duration by dynamic fluid injection. Indeed, fluid injection can reduce the effective normal stress on asperities of the fault where the events occur. As fluid is injected or when aseismic slip is active, rapid, strong, heterogeneous increase of pore pressure in the fault can dramatically decrease the effective normal stress and thus the strength of the fault, triggering seismic slip on localized asperities. Under these conditions, [START_REF] Lengliné | Fluid-induced earthquakes with variable stress drop[END_REF] and [START_REF] Lin | Evidence for non-self-similarity of microearthquakes recorded at a taiwan borehole seismometer array[END_REF] argue that the heterogeneity of fluid injection is likely to produce the wide range of stress-drops they observe. As noted by the previous papers, reduced effective normal stress have also been suggested as a cause of the low values of stress drop of VLFEs in Japan [START_REF] Ito | Very low frequency earthquakes within accretionary prisms are very low stress-drop earthquakes[END_REF]. In the light of the similarity of our observations and the ones presented in the discussed studies, dynamic triggering of LFEs by fluid injection could be a relevant mechanism to explain the LFEs' stress drop values.

In order to account for the observation of events' constant duration, [START_REF] Lin | Evidence for non-self-similarity of microearthquakes recorded at a taiwan borehole seismometer array[END_REF] proposes that slip could only be maintained while the pore pressure is high enough to allow for seismic rupture, that is until the injection stops and the high pore pressure transient diffuses. The duration of the recorded seismic event would thus be controlled by the fluid transport properties in the fault zone. Those conditions of rapid, localized fluid injection are plausible within large active fault zones [START_REF] Shapiro | Low-frequency earthquakes and pore pressure transients in subduction zones[END_REF]. In subduction zones, fluids released by the dehydration of hydrous minerals in the downgoing slab circulate within a permeable channel formed along or within the fault zone [START_REF] Angiboust | Effect of fluid circulation on subduction interface tectonic processes: Insights from thermo-mechanical numerical modelling[END_REF][START_REF] Van Keken | Subduction factory: 4. depth-dependent flux of h2o from subducting slabs worldwide[END_REF]. If the permeability or width of the channel has strong lateral heterogeneities, [START_REF] Shapiro | Low-frequency earthquakes and pore pressure transients in subduction zones[END_REF] suggests that strong, localized pore pressure transients can be formed. The strong, transient gradients of pore pressure can act as a source of seismic waves. The duration of such a source would thus be the characteristic time taken by a pore pressure transient to diffuse. [START_REF] Shapiro | Low-frequency earthquakes and pore pressure transients in subduction zones[END_REF] shows that as with such a mechanism, this duration is only governed by the transport properties of the fault zone and independent of the source magnitude. Although more precise observations of the source radiation pattern of LFEs and a substantial modeling effort are needed to confirm this hypothesis, it is consistent with our observations. This work is published in Journal of Geophysical Research: Solid Earth, DOI:10.1029/2019JB019099.

Another lead to account for the apparent magnitude-independent, low-frequency content of LFEs is to interpret it as the consequence of specific, near-source propagation effects [START_REF] Bostock | On corner frequencies, attenuation, and low-frequency earthquakes[END_REF]. In this study the authors assume the existence of a narrow region with a very strong shear wave attenuation around LFEs source and suggest that the observed depletion in high frequencies can be related to this strong near-source attenuation of seismic waves. Based on accounts of high compressional to shear wave velocity ratio V P /V S in seismic cross-section of subduction zones [START_REF] Bostock | Low frequency earthquakes below southern vancouver island[END_REF][START_REF] Shelly | Low-frequency earthquakes in shikoku, japan, and their relationship to episodic tremor and slip[END_REF], tidal triggering of LFE activity [START_REF] Nicholas | Inferring fault rheology from low-frequency earthquakes on the san andreas[END_REF][START_REF] Royer | Tidal modulation and triggering of low-frequency earthquakes in northern cascadia[END_REF][START_REF] Rubinstein | Tidal modulation of nonvolcanic tremor[END_REF], e.g.] and dynamic triggering by teleseismic waves [START_REF] Rubinstein | Seismic wave triggering of nonvolcanic tremor, episodic tremor and slip, and earthquakes on vancouver island[END_REF], the source region of LFEs is thought to be under nearly lithostatic fluid pressures. Laboratory studies indicate that such high fluid content can also be the cause of strong shear wave attenuation [START_REF] Tompkins | Ultrasonic p-and s-wave attenuation in oceanic basalt[END_REF].

Conclusions

We analyzed a subset of carefully selected LFEs from the cataloged events of Guerrero, Mexico [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in guerrero, mexico[END_REF]. The subset of LFEs samples a relatively broad range of seismic moments, between M w 1.7 and M w 3, with corner-frequencies clustering around 3 Hz, between 1.5 and 6 Hz. Overall, we find corner frequences to be very weakly dependent on seismic moments for LFEs in Guerrero. When trying to describe this relationship with the power law f c ∝ M 0 -α , we retrieve α = 0.12 ± 0.11 and α = 0.054 ± 0.034 for the two sets of corner frequency measurements. In both cases those scalings are significantly different from a self-similar rupture, for which α = 1/3.

Recent studies show that the clear departure from self-similarity measured for LFEs in Guerrero and Cascadia is not a necessary feature of LFEs in all tremor regions across the world. Exploring this regional specificity might bring significant insight on LFEs source physics and lead us to understand more about how deep fault environment affects slow slip processes.

A Effects of bandpass filtering on events corner frequencies A.1 Detection bandpass filter

The LFEs of the catalog used in this study were detected using a matched-filter search approach, with template waveforms band-pass filtered between 1-2 Hz [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in guerrero, mexico[END_REF].

A potential concern is that a matched-filter search using a narrow band-pass would only identify events with a characteristic frequency between 1 and 2 Hz. Assuming it was the case, it would bias our analysis towards the observed magnitude-independent corner frequencies.

However, the use of matched-filter search does not constrain the characteristic frequencies of the event between 1-2 Hz, it simply limits detections to signals with energy within this frequency band. Several studies using matched-filter search to detect classic earthquakes are able to detect events with a range of magnitudes, that, through the earthquake moment-duration scaling law, corresponds to a range of characteristic frequency wider than the chosen band-pass [71, figure 7.(b)],[22, figure 5]. In our case most events indeed have corner frequencies above 2 Hz and up to 6.5 Hz (see Table 1, Figure 5).

A compact form of the Brune spectrum for ground motion velocity can be used to compare the detectability of two events: v(f ):

v(f ) = 2πC.M 0 .f 1 + ( f fc ) 2 , ( 14
)
where C is a constant relating to the velocity structure of the medium and the source focal mechanism, M 0 the seismic moment of the event and f c its corner frequency. For any two events with f c1 < f c2 but identical moment, for any given frequency of detection f , the corresponding velocity amplitudes will follow v 1 (f ) < v 2 (f ). In other words, for a given magnitude, events with a higher characteristic frequency than events we detected do not have a lower probability to be detected because of the band-pass filter. Therefore, within the range of moments that we measured, the detection band-pass filter should not bias corner frequencies by constraining them at low values, close to 2 Hz. Using the same demonstration, events of characteristic frequency below 1 Hz and seismic moment within the range of moments we measure have a low probability to be detected by the matched-filter search. The matched-filter search band-pass thus imposes a lower bound on detections' corner frequencies in the range of moments we measure, but no higher bound. For events above 8 Hz however, the presence of high-frequency noise would bias the measure of corner frequencies, motivating our choice to limit our analysis to the 1-8 Hz frequency band. As a conclusion, in the seismic moment range of our events (8 × 10 11 -1 × 10 13 N.m), the measure of corner frequency is not biased by the detection process, for values of corner frequency between 1 and 8 Hz. A frequency range of around 1-10 Hz and range of moments 1.5 orders of magnitude wide were enough for [START_REF] Supino | Self-similarity of low-frequency earthquakes[END_REF] to measure a steeper dependence of corner-frequency on seismic moment, and should thus be enough in this study to estimate a moment-duration scaling without significant artifacts.

A.2 Analysis bandpass filter

In Figure S1, we see that the tremor seismic energy seems really bounded between 0.8 Hz and 8 Hz. Above this band, high-frequency noise, possibly anthropogenic, is active, below this band, the seismic time-series are dominated by low-frequency noise leaking from microseismic modes. The choice of the 1-8 Hz band-pass filter used to measure corner-frequencies is motivated by this observation.

It is possible that corner frequencies of the event we detected are lower than 1 Hz, and that the 1-8 Hz band-pass filter we use biases our measure of the corner frequencies to constrain them in a narrow band above 1 Hz. We tested the influence of the lower cutoff frequency of the band-pass filter on the measured corner frequency. Results of this test are shown in Figure S2.

The distribution of corner frequencies is not drastically modified down to a lower cutoff of 0.2 Hz, where corner frequencies suddenly regroup around 0.2 Hz. This sudden shift is due to the filter band not blocking pollution of micro-seismic noise anymore.

In addition to it, corner frequencies are slightly dragged to lower values, as the lower cutoff decreases from 1 Hz to 0.5 Hz. This should be expected: more and more pollution from micro-seismic noise is allowed in records of events on stations that witness higher amplitudes of micro-seismic noise, due to proximity to the ocean or site effects. On those records, the characteristic frequency will be picked lower, due to the influence of micro-seismic noise. Figure 4 (d.) is a good illustration of this: the peak of the un-filtered spectrum is located at 1 Hz due to micro-seismic noise influence.

The 1-8 Hz is thus chosen as a compromise. The lower cutoff is taken far enough from micro-seismic frequencies so that the event signal is not polluted, even at stations where micro-seismic noise has the highest amplitude. Ultimately, it does constrain the corner frequency between 1 and 8 Hz, but this seems reasonably wide.

B Moment-frequency distribution of LFEs

In Figure 7, we represent the complementary cumulative distribution function (abbreviated as CCDF, also called survival function) of seismic moments, P (M 0 > M 0 ). It is computed with M 0 regularly spaced in the magnitude range. The higher-moments tail of the distribution can be described as a power-law distribution, with the following probability density function (PDF) and CCDF:

p(M 0 ) = β M 0thr M 0 M 0thr -β-1 P (M 0 > M 0 ) = M 0 M 0thr -β (15)
where M 0thr is the lower bound for which the power-law distribution is valid. The power law exponent β relates to the b-value of the Gutenberg-Richter law following: b = 3 2 β. For moments above a threshold M 0thr , we perform an un-weighted, least-square fit of the logarithm of the moments CCDF to a powerlaw. We find that β = 1.62, and a corresponding b-value of b = 2.43. For comparison, b-values of classical seismicity are usually around 1. The threshold magnitude M wthr = 2.32 (corresponding to M 0thr ) is selected so as to minimize the Kolmogorov-Smirnov statistic for the moments' distribution and the fitted power law.

The deviation from a power law at lower moments is often attributed to observational limitations: below a completeness magnitude, events are too weak to all be detected, and it manifests as a deviation from the otherwise assumed power law distribution of moments. Here, the completeness magnitude would thus be M wthr = 2.32. The departure from a power-law at lower moments can also be the manifestation of a characteristic scale of magnitude of the events. Several studies report observations of LFEs size distribution being best modeled by exponential laws and thus of a characteristic scale controlling the LFE size distribution [START_REF] Chamberlain | Low-frequency earthquakes reveal punctuated slow slip on the deep extent of the alpine fault, new zealand[END_REF][START_REF] Sr Chestler | Evidence for a scale-limited low-frequency earthquake source process[END_REF]. In our case, a normal distribution can be fitted to the distribution of moments below M wthr , parametrized by a mean magnitude M w = µ = 2.

Figure 1 :

 1 Figure1: Spatial and temporal distribution of the selected LFEs -Selected LFEs are plotted as white and blue dots, the bulk catalog of[START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in guerrero, mexico[END_REF] is in gray. Subduction interface depth contours are displayed in black[START_REF] Kim | Geometry and seismic properties of the subducting cocos plate in central mexico[END_REF]. The focal mechanism corresponds to an average mechanism computed from cataloged LFE detections[START_REF] Frank | Low-frequency earthquakes in the mexican sweet spot[END_REF]. (c.) shows the time-dip distribution of the events.

5 <

 5 I 1-2Hz < 0.92 T low < A < 2.10 -6 m.s -1

Figure 2 :

 2 Figure 2: Selection of LFEs -(a.) Schematic recapitulation of the selection procedure. (b.) Amplitudeimpulsivity distribution of all detections recorded on the East component of the station TONA. Each dot is an LFE, its color corresponds to the number of channels on which it is recorded that pass the criteria. Selected LFEs are white dots: they are recorded on 10 or more channels passing the criterion and fulfill the recurrence and correlation requirements. The amplitude and first impulsivity criteria (equations 4 and 3) are shown with dashed lines. LFEs in the criterion domain pass the criterion on TONA.HHE, but 9 other records passing the criterion on other channels of the network are needed for an LFE to be selected.

Figure 3 :

 3 Figure 3: Selected LFEs' waveforms -The 4 most impulsive waveforms of LFE 379 and 687 are displayed here, 1-2 Hz band-pass filtered ((a.) and (c.)), and between 1 and 5 Hz ((b.) and (d.)).The blue section of the waveform shows the time window used to define the detected signal, the red patch shows which channel's waveforms are used as an example for source parameters characterization in Figure4. Impulsivity and amplitude are indicated for each channel.

Figure 4 :

 4 Figure 4: Determination of the seismic moment M 0 and corner frequency f c , for an example lowfrequency earthquake -Example LFE: nb. 379, channel HHE of station TONA. (a.) Ω 0 is determined as the red plateau of the displacement spectrum. (b.) s Ω 0 is determined as the maximum of the synthetic displacement spectrum, taking into account attenuation. (c.) The corner frequency f c is measured as the mean instantaneous frequency (MIF) over the waveform time window (see Figure3for the LFE waveforms), (d.) and as the maximum of the velocity spectrum (MVS) of the waveform. In gray, the un-filtered velocity spectrum, allowing to assess the effects of the 1-8 Hz filter on the spectrum.

Figure 5 :

 5 Figure5: Distribution of the measured source parameters -The distribution is displayed and colored according to the point density for both methods of determining f c : (a.) with mean instantaneous frequency (MIF), (b.) as the frequency of the maximum of the velocity spectrum (MVS). Weighted least-square fits to the binned moments (see Figure6) are shown as dotted lines, in each respective plot. Generic values for M 0 and f c measured in other regions are displayed along the distribution (1.[START_REF] Thomas | Constraints on the source parameters of low-frequency earthquakes on the san andreas fault[END_REF], 2.[66], 3. [32], 4.[START_REF] Bostock | Magnitudes and momentduration scaling of low-frequency earthquakes beneath southern vancouver island[END_REF]) The self-similarity scaling law for classic earthquakes is represented by the purple patch. The higher bounding line corresponds to a stress drop of 100 MPa and the lower to a stress drop of 0.1 MPa[START_REF] Allmann | Global variations of stress drop for moderate to large earthquakes[END_REF][START_REF] Douglas | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF][START_REF] Madariaga | Dynamics of an expanding circular fault[END_REF].

Figure 6 :

 6 Figure6: Binned distribution and uncertainties estimation -The measured M 0 -f c is displayed and colored according to point density, for both the mean instantaneous frequency (MIF) and maximum of velocity spectrum (MVS). Seismic moment and corner frequency are binned, and the weighted standard deviation is shown with the error bars (see text for details). Weighted least-square fits to the binned data are displayed, the fit resulting in the purple scaling has α = 1 3 as a constraint.

Figure 7 :

 7 Figure 7: Moment-frequency distribution of the selected LFEs -The seismic moment complementary cumulative distribution function (CCDF) of our LFEs is modeled by a Gaussian below M wthr 2.32, with mean magnitude M w 2.4, and by a power law above M wthr , with an exponent corresponding to a b-value of 2.4.

  23 and a standard deviation This work is published in Journal of Geophysical Research: Solid Earth, DOI:10.1029/2019JB019099. σ = 0.19. The mean magnitude M w = 2.23 defines a characteristic scale of the lower-bound of the moment distribution. It is either due to a scale limitation of LFEs in Guerrero, or a manifestation of the limitation of events' detection.

Table 1 :

 1 Measured LFE source parameters statistics analysis values Method \Frequencies (Hz)

	Min. Mean Max.
	Instantaneous frequency 2.0	3.18	5.52
	Maximum of velocity sp. 1.07	2.81	6.47

This work is published in Journal of Geophysical Research: Solid Earth, DOI:10.1029/2019JB019099.

Acknowledgements Data and materials availability: The source parameters data we generate in this study is available on the Slow Earthquake Database (http://www-solid.eps.s.u-tokyo.ac.jp/ sloweq/). The complete low-frequency earthquake catalog compiled by [20] can be found on the same archive. The MASE seismic dataset (https://doi.org/10.7909/C3RN35SP) used in this study is openly available from IRIS data center under the network code TO. Software use: One figure was made with the Generic Mapping Tools (GMT) [72]. We used the Python library ObsPy [44] for seismic data handling and the mtspec library [55] for multitaper spectral analysis. Research funding and support: This study was supported by the Russian Ministry of Education and Science (grant N 14.W03.31.0033) and by the European Research Council under the European Union Horizon 2020 research and innovation program (grant agreement no. 787399 -SEISMAZE). Numerical computations were performed on the S-CAPAD platform, at Institut ie Physique du Globe de Paris (IPGP), France. The authors confirm that there are no real or perceived conflicts of interest associated with this work.