The CARMENES search for exoplanets around M dwarfs Behaviour of the Paschen lines during flares and quiescence ⋆

The hydrogen Paschen lines are known activity indicators, but studies of them in M dwarfs during quiescence are as rare as their reports in flare studies. This situation is mostly caused by a lack of observations, owing to their location in the near-infrared regime, which is covered by few high-resolution spectrographs. We study the Pa β line, using a sample of 360 M dwarfs observed by the CARMENES spectrograph. Descending the spectral sequence of inactive M stars in quiescence, we find the Pa β line to get shallower until about spectral type M3.5V, after which a slight re-deepening is observed. Looking at the whole sample, for stars with H α in absorption, we find a loose anti-correlation between the (median) pseudo-equivalent widths (pEWs) of H α and Pa β for stars of similar effective temperature. Looking instead at time series of individual stars, we often find correlation between pEW(H α ) and pEW(Pa β ) for stars with H α in emission and an anti-correlation for stars with H α in absorption. Regarding flaring activity, we report the automatic detection of 35 Paschen line flares in 20 stars. Additionally we found visually six faint Paschen line flares in these stars plus 16 faint Paschen line flares in another 12 stars. In strong flares, Paschen lines can be observed up to Pa14. Moreover, we find that Paschen line emission is almost always coupled to symmetric H α line broadening, which we ascribe to Stark broadening, indicating high pressure in the chromosphere. Finally we report a few Pa β line asymmetries for flares that also exhibit strong H α line asymmetries.


Introduction
Stellar activity is ubiquitous in M dwarfs. It is frequently studied by observing activity tracers, such as the X-ray flux (Pizzolato et al. 2003;Foster et al. 2022;Magaudda et al. 2022) or chromospheric line fluxes (Gomes da Silva et al. 2021). Many lines sensitive to the chromosphere and transition region are located in the ultraviolet (UV), such as the Lyα line (Youngblood et al. 2016), which, together with the bulk UV emission, is a crucial ingredient to assess the habitability of exoplanets (Youngblood et al. 2017). However, all important UV lines are not observ-able from the ground and partly not even with current satellites. Fortunately, also the observable visible range includes a number of activity-sensitive lines, such as the Ca ii H&K lines and the related near-infrared (NIR) triplet lines of Ca ii, which are not as sensitive to activity as the former though (Martínez-Arnáiz et al. 2011;Martin et al. 2017;Mittag et al. 2017;Pavlenko et al. 2019). Yet, the most widely used activity indicator in M dwarfs remains the Hα line (Gizis et al. 2002;Walkowicz & Hawley 2009;Lodieu et al. 2011;Schöfer et al. 2019).
Extensive studies of spectral activity tracers in the infrared have only recently become feasible with the advent of NIR high-resolution spectrographs, such as CARMENES (Quirrenbach et al. 2020). Although chromospheric indicator lines are often less prominent in the infrared, there are a few known ex-Article number, page 1 of 22 arXiv:2308.07685v1 [astro-ph.SR] 15 Aug 2023 A&A proofs: manuscript no. carmenes_pa amples such as the K i doublet, which was reported to be magnetically sensitive (Fuhrmeister et al. 2022;Terrien et al. 2022), and the He i infrared triplet (IRT), which is a long-known activity tracer in the Sun and other stars (Vaughan & Zirin 1968;Sanz-Forcada & Dupree 2008;Andretta et al. 2017). The Paschen (Pa) series of hydrogen, which is known to react to strong flares, is also located in the infrared regime. In the solar context, Paschen lines were used mainly to determine the electric field in prominences via the Stark effect (Foukal et al. 1987;Casini & Foukal 1996). For stars, there are a number of flare observations, in which members of the Pa series were detected in emission. For example, Liebert et al. (1999) observed the M9.5 dwarf 2MASSW J0149090+295613 during a mega-flare event, which allowed them to detect Pa 7 1 through 11 in emission, and observe their decay during another three observations within about 30 min. Schmidt et al. (2007) observed a large flare on the M7 star 2MASS J1028404−143843, which showed strong continuum enhancement even at 10 000 Å and exhibited Pa 8 to 11 in emission. Fuhrmeister et al. (2008) found Pa 7 to 11 in emission during a mega-flare on CN Leo (M5.5), which decayed fully within about 10 minutes. The Pa 7 to 9 lines were also reported in emission for a medium sized flare on Proxima Centrauri by Fuhrmeister et al. (2011). In a dedicated flare search covering about 50 hours of observation time of the active M dwarfs EV Lac, AD Leo, YZ CMi, and vB 8, Schmidt et al. (2012) detected 16 flares, out of which three showed infrared emission including Pa 5 (= Paβ) and Pa 6 (= Paγ). Kanodia et al. (2022) analysed high-resolution flare spectra of the M8 dwarf vB 10 covering Pa 6, 7, 11, and 12. While Pa 12 was only marginally detected, Pa 6 showed indications of a weak red asymmetry and persisted to be in emission in a second exposure about 20 minutes later. In a study of Hα line asymmetries using CARMENES data, Fuhrmeister et al. (2018) searched for Paβ line emission in 36 spectra exhibiting flare-induced Hα wing asymmetries and reported 9 weak detections of Paβ emission.
Examples of non-flare studies are more rarely found. In a publication by Klein et al. (2020), who studied the M1 dwarf AU Mic, it was found that the stellar rotation period could be recovered using the He i IRT and the Paβ lines. These authors also reported that the origin of the He i IRT lines seems to be more concentrated toward equatorial latitudes, while the Paβ line is primarily formed at polar regions on AU Mic. Turning to large stellar samples, Schöfer et al. (2019) used CARMENES data in a comprehensive activity study to compare the relation between different activity indicators in M dwarfs using a spectral subtraction technique. While they did not find any correlation between the pseudo equivalent-width (pEW) of Paβ and that of the Hα line, they reported deeper (excess) absorption in the Paβ line for the most active stars as measured by Hα and of spectral type earlier than M4.0.
In the following, we utilize the unique database of M dwarf spectra obtained with the CARMENES spectrograph to study the Paschen lines along the M dwarf sequence in more detail. First, we analyse their behaviour in the quiescent activity levels of the stars along the M sequence and, second, we investigate the Paschen lines for the detected flares. Specifically, we address the question of how often flares are detectable by Paschen emission in comparison to Hα emission and what physical parameters favour Paschen line emission. 1 We use the notation Pa N to indicate the Paschen line with the upper level N and only refer to Paβ, Paγ, and Paδ, which are equivalent to Pa 5, Pa 6 and Pa 7, using the Greek nomenclature, unless in direct comparison to higher level lines.

Observations
All spectra used in this study were taken with the CARMENES spectrograph, installed at the 3.5 m Calar Alto telescope (Quirrenbach et al. 2020). CARMENES covers the wavelength range from 5 200 to 9 600 Å in the visual channel (VIS) and from 9 600 to 17 100 Å in the near-infrared channel (NIR). The instrument provides a spectral resolution of ∼ 94 600 in VIS and ∼ 80 400 in NIR. While the CARMENES data are obtained mainly for planet search, they are also a resource for studies of stellar parameter determination and activity. A large part of the data (years 2016-2020) have become public (Ribas et al. 2023). The data are especially well suited also for other purposes, since the CARMENES sample is biased only marginally. Since Gaia data were not available at the time of building the CARMENES guranteed time observations M-dwarf sample (Alonso-Floriano et al. 2015;Reiners et al. 2018;Ribas et al. 2023), the CARMENES consortium selected the brightest stars (in J band) for each spectral subtype that were observable from Calar Alto (i. e. δ > -23 deg) and that did not have any known close companion at ρ < 5 arcsec. As a result, the only bias in our sample is Malmquist's, by which overluminous young stars in stellar kinematic groups are overrepresented in our target list. However, most of these are very active and have, therefore, a large RV jitter that impedes reaching the main scientific objective of CARMENES, which is the search for Earth-like planets in the habitable zone of M dwarfs. As a result, the consortium discontinued observations of a few "RV-loud" stars at the beginning of the program, after a minimum number observations (Tal-Or et al. 2018). Nevertheless, we also include 27 of the 31 "RV-loud" stars in the investigation in this work.
In our analysis, we considered a sample of 360 M dwarfs observed by CARMENES, resulting in more than 19 000 spectra taken before September 2022. We excluded known binaries (Baroch et al. 2018;Schweitzer et al. 2019;Baroch et al. 2021), which may hamper our analysis by orbit-induced line shifts. Moreover, the CARMENES consortium has invested a considerable effort in determining stellar parameters of the target stars, like spectral types (Alonso-Floriano et al. 2015), luminosities and colours (Cifuentes et al. 2020), photospheric parameters, namely T eff , log g, and [Fe/H] (Passegger et al. 2018Marfil et al. 2021;Passegger et al. 2022), rotation velocities , rotation periods (Díez Alonso et al. 2019;Shan et al. 2023), magnetic fields , and masses and radii . Furthermore, stellar activity has been studied using the CARMENES high resolution spectroscopic data. In particular, Hα was investigated (Fuhrmeister et al. 2018;Schöfer et al. 2019), along with other activity sensitive lines. Examples are the He i infrared triplet (Fuhrmeister et al. 2019 or the optical and infrared K i doublets (Fuhrmeister et al. 2022). Also other spectral indicators Schöfer et al. 2019Schöfer et al. , 2022 were studied. In this work we make extensive use of the results obtained in these publications and refer to them for further details.
The stellar spectra were reduced using the CARMENES reduction pipeline (Zechmeister et al. 2014;Caballero et al. 2016). Subsequently, we corrected them for barycentric and systemic radial velocity motions and carried out a correction for telluric absorption lines (Nagel et al. 2019) using the molecfit package 2 . No correction for airglow emission lines was attempted, although they can play a role near the Paschen lines and may be shifted into the integration ranges used for their analysis. Be- cause of the amount of available data and the difficulty in dealing with the lines automatically, we decided to identify and remove the affected spectra later on in the analysis.
The CARMENES instrument does not cover the Pa α line at 18756.4 Å, but all other members of the Paschen series are covered, including Pa 5 (Paβ) at 12 821.578 Å, Pa 6 (Pa γ) at 10 941.17 Å, and Pa 7 (Pa δ) at 10 052.6 Å. The Paschen series ends at about 8250 Å, which is located already in the VIS channel of CARMENES. O 2 airglow lines are found near the Paβ line at (vacuum) wavelengths of 12 819.46, 12 822.43, 12 824.78 Å, with the latter line being the strongest (Oliva et al. 2015).

pEW measurements of activity indicators
To assess the activity state of the stars in each spectrum, we employed pEW measurements. The spectra of M dwarfs do not show an identifiable continuum because of the abundance of molecular absorption lines. These pEW measurements were then used to search for flares in Hα and Paβ.
To give an overview of the appearance of the Paβ line, we show examples of it along the M dwarf sequence in Fig. 1 for stars with Hα in absorption. In this sequence, the strongest Paβ absorption is observed in the M0.0 star. Absorption subsequently weakens until a minimum is reached for the shown M3.5 star, after which the Paβ lines deepens again; a more detailed discussion is provided in Sect. 4.1.
To quantify the level of absorption or emission in the lines, we computed pEWs of the Hα line, the bluest and middle lines of the Ca ii IRT, as well as the Paβ, Paγ, and Paδ lines. We considered Hα to be in absorption if the pEW value was larger than −0.6 Å, while lower values marked Hα in emission. This threshold was already used by Fuhrmeister et al. (2019) and is in between other adopted values as −0.5 Å  or −0.75 Å (West et al. 2011). If Hα is in absorption or emission has traditionally been used to discriminate between inactive or active stars and we also use it here to split up our sample in this sense.
For the pEW computation, we list the central wavelength, full width of the line integration window, and the location of the two reference bands in Table 1; for a more detailed description of pEW measurements of chromospheric lines we refer to Fuhrmeister et al. (2023). While the reference bands are typically blue-and redward of the central wavelength, both are located blueward for the Paβ line, since it is located near the red edge of one of the spectral CARMENES orders and we did not want to use a reference band in a different order. While we used 1.5 Å wide line integration bands for the Paγ and Paδ lines, we opted for a narrower 0.6 Å wide integration band of Paβ to minimize contamination by airglow. This is similar to the even narrower 0.5 Å band used by Schöfer et al. (2019) for Paβ. The widths of 1.6 Å and 0.5 Å for the Hα and Ca ii IRT lines were also used by Fuhrmeister et al. (2020).
For all studied lines, emission during flares may be broader than the line integration band. Therefore, in extreme cases, the full variability range may not be represented with our choice of integration ranges. Moreover, rotation rates higher than about v sin i = 15 km s −1 will affect the pEW measurements by shifting flux out of the integration bands (this threshold is exceeded by 20 of our sample stars). Nevertheless, we found the chosen integration bands be suitable for identifying variability, which we are most interested in.

Search for Paschen and Hα line flares
To study the Paschen lines during flares, flares with a reaction of the Paschen lines need to be identified in the first place. The CARMENES observing schedule does rarely produce consecutive spectra of the same star, but observations of the same star are typically separated by some days. To facilitate a flare search, we computed for each star the median, µ, of the pEW measurements for each chromospheric line and the median average deviation about the median (MAD). We list these values together with some basic stellar parameters in Table 2 for each star. The MAD yields a robust estimator of the standard deviation. If MAD(pEW(Hα)) and σ(pEW(Hα)) denote the MAD and standard deviation of the time series of pEW measurements of Hα σ(pEW(Hα)) = 1.4826 × MAD(pEW(Hα)) . (1) The same nomenclature is used for the other lines.
In a first step, we searched for flares indicated by Hα and the Ca ii IRT lines. We accepted a spectrum as flaring (and call this an Hα flare) in case of combined Hα and Ca ii IRT excursions, viz., (i) pEW(Hα) < µ(pEW(Hα)) − 3σ(pEW(Hα)) and (ii) pEW(Ca IRT 1,2 ) < µ(pEW(Ca IRT 1,2 )) − 3σ(pEW(Ca IRT 1,2 )) , where the condition (ii) must apply to at least one of the considered Ca ii IRT lines. Using only the Hα line in the search worked well for inactive stars, which exhibit pronounced but seldom flares and show a rather stable Hα absorption line otherwise, but not for stars showing persistent, strong variability in Hα. Since for many of these stars, the Ca ii IRT lines showed less pronounced variations outside of flares, we coupled the search for flares showing up in Hα to Ca ii IRT. Flares also showing up in the Paβ line were identified by requiring that condition (i) and (ii) are met, (i. e. the star is flaring in Hα) and an equivalent 3σ condition for the pEW of the Paβ line is fulfilled. In these so identified Paβ flares we additionally searched for Paγ and Paδ flares, applying the 3σ condition one more time.
By this method, most cases of spurious flares induced by statistical noise are suppressed. Additionally, the coupling with an Hα criterion removes spurious flares caused by airglow contamination in Paβ. We nevertheless inspected all Paβ flare detections by eye. In the follwoing we call all Paβ flares fulfilling our flare criteria 'automatically detected' and flares that pass the visual inspection 'visually confirmed'.  The line is absent in photospheric spectra. Moreover, it can be seen in Fig. 1 that around spectral type M2.0 V another absorption feature blueward of the Paβ line at about 12821.4 Å starts to emerge, deepening for later spectral types. The feature can also be seen in the PHOENIX photospheric spectrum for the M5.0 V star shown in Fig. A.1, although it is not as deep as observed, which is not uncommon for a molecular feature.

Results and discussion
Concerning the origin of the Paschen lines, we distinguish between the lines observed in absorption or in emission. Cram & Mullan (1979) found that for the Hα line the line source function is controlled by photoionization (and recombination) for cases where it is in absorption. Therefore, the n = 3 level as groundstate of the Paβ line is populated by Hα absorption and by recombination. On the other hand, for flaring states, where we found the line in emission, it must be collisionally controlled, as Cram & Mullan (1979) argued for Hα as well.

The Paschen lines along the M dwarf spectral sequence
To generalise the impression from the example sequence in Fig. 1, we show in Fig. 2 the distribution of all median(pEW(Paβ)) per star as a function of the effective temperature T eff , adopted from Cifuentes et al. (2020, from spectral energy distribution fitting) and Marfil et al. (2021, from spectral synthesis). Fig. 2 shows that the Paβ line becomes shallower (i.e., yields lower pEWs) for lower T eff (or later spectral type) until a turning point is reached at about T eff < 3400 K, which corresponds to spectral types of about M4.0 V. In Fig. 3 we compare the median observed pEWs of the Hα and Paβ lines of all 360 sample stars. T eff is shown colourcoded in 100 K intervals, to emphasise the temperature dependence of the Paβ line. Looking at the stars with Hα in absorption (pEW(Hα > −0.6 Å; which we call inactive stars for the purpose of this study) in Fig. 3, an anti-correlation between µ(pEW(Paβ)) and µ(pEW(Hα)) can be noticed for each effective temperature interval. To quantify this impression, we calculated Pearson's correlation coefficients, r, for these samples, and obtained val- Table 2. The full table is provided at CDS. We show here the first ten rows as a guidance.
Article number, page 5 of 22 A&A proofs: manuscript no. carmenes_pa ues between −0.42 and −0.91 with p-values between 0.05 and 10 −5 , for 4000 K > T eff > 3200 K indicating fair to highly significant correlations. Only for the highest temperature interval with r = −0.45 and p = 0.13 and for the lowest temperature interval with r = 0.44 and p = 0.08 correlations are questionable (and would be positive for the lowest temperature stars). Thus, in general more absorption in Hα is on average associated with less absorption in Paβ in the stars with Hα in absorption for most temperatures. Since in our stellar sample typically a larger pEW(Hα) is connected to less activity, the Paβ line deepends for higher activity levels for the here considered inactive stars. This finding is in line with the analysis by Schöfer et al. (2019).
Stars with Hα in emission (which are called active stars traditionally) are only available in meaningful numbers in our sample for T eff < 3400 K and there is no comparable correlation for these; only for the 3200 < T eff < 3300 K interval we find a fair correlation with r = 0.52 and p = 0.01, the other two temperature intervals show no correlation with r between −0.03 and 0.15 and p > 0.40. Nevertheless, for stars with T eff >3600 K, the pEW(Paβ) increases further for stars with Hα in emission compared to the pEW(Paβ) of stars with Hα in absorption. For stars with 3200 < T eff < 3600 K and Hα in emission, pEW(Paβ) saturates at the highest values found for stars with Hα in absorption. For stars with T eff < 3200 K saturation effects play a role for some of the stars, while others show higher or lower values of pEW(Paβ).
Generally, for the coolest stars in our sample, the spread in pEW(Paβ) is largest. Moreover, for these stars the most inactive stars with the lowest pEW(Paβ) values may be not present causing the mean apparent re-deepening of the Paβ line together with the additional absorption feature at 12821.4 Å. Nevertheless, as can be seen in Figs. 2 and 3, some of the coolest stars resume the trend to lower pEW(Paβ). These low values are found in more active stars and we interpret this as fillin-in of the line suggesting that the Paβ line is very sensitive to the pressure in the chromosphere.
For 187 stars of our sample a rotational period is known. We therefore compare also the median(pEW(Paβ)) to the rotational period for these stars (see Fig. 4). As discussed above, for the stars with higher effective temperature, which are generally more inactive, a deepening of the Paβ line can be noticed towards shorter rotation periods (i. e. higher activity levels). Only for the coolest stars deepening and fill-in or saturation is observed as one proceeds to shorter periods.

Time series of individual stars
To study the relation between the pEWs measured in the Hα and Paβ lines in individual stars, we also computed Pearson's correlation coefficients for the time series of pEW(Hα) and pEW(Paβ) on a by-star basis and show the resulting values of r in Fig. 5, where significant results (p < 0.005) are highlighted by colour. Stars with Hα in emission tend to show a significant positive correlation between pEW(Hα) and pEW(Paβ), while stars with Hα in absorption more often show significant anti-correlation. The latter cases are mostly stars with T eff >3400 K, while the former are usually cooler stars, and one should keep in mind that our sample contains few stars with T eff >3400 K and Hα in emission. Also, many of the stars showing positive correlations are affected by flaring activity in the Paschen lines (see Sect. 4.2.1), which often dominates the variability in pEW(Paβ) and, thus, drives the correlation. Moreover, it is an interesting question, if detectable rotational modulation is imprinted on the Paβ time series. We therefore computed a generalized Lomb-Scargle periodogram (Zechmeister & Kürster 2009) for the Paβ time series and accepted periods between 1.5 and 150 days and a false-alarm probability smaller than 0.005 as significant rotational modulation. We found six stars fulfilling these criteria, which we list in Table 3. For three of these stars, there is no known rotation period. For one star a conflicting rotation period is known; for another star (J16343+571 / CM Dra; spectral type M4.5 V, eclipsing binary) we found a period of 2.4 days, while a period of about half this value of 1.27 days is known for this star for the mutual orbital period (Doyle et al. 2000). For the last star (J03133+047 / CD Cet; spectral type M5.0 V) we found a period of 127.2 days, while a period of 126.2 days was determined by Newton et al. (2016). A more detailed study found a period of 170 +19 −38 days using photometry and about 134 days using spectroscopy (Bauer et al. 2020). These findings show again that the Paβ line is sensitive to activity, but not as sensitive as other tracers. This may -especially for period search -be partly caused by the problems of obtaining pEWs free of the influence of telluric and airglow lines or the artefacts from their removal. Anyway, rotational modulation in M dwarfs is traced best with photometric variations (Irwin et al. 2011;West et al. 2015;Suárez Mascareño et al. 2016;Díez Alonso et al. 2019).

Overview
Applying the automatic flare search described in Sect. 3.2, we found 357 Hα flares in 153 stars, and 46 Paβ flares in 30 stars. We summarize these and all the numbers in this Section in Table 4. We examined all Paβ flares by eye and removed 11 for which we found airglow to remain a problem. The other stars all show the Paβ line in emission. This excludes confusion with high amplitude rotational modulation, since Paβ emission certainly involves high pressure, see also the discussion in Sect. 4.2.3. For these stars we automatically detect 15 stars with 24 Paγ flares and 15 stars with 24 Paδ flares. These are not the same, though. While for most flares with Paγ emission, Paδ emission is also detected, there are six flares, where no Paδ was detected. There are another six flares, where Paδ was detected despite no Paγ emission. For these latter six flares the Paδ detection is correct and Paγ was not detected due to noise in some spectra of the three affected stars, which enlarges the MAD incorrectly and hinders the automatic Paγ detection. Therefore we manually correct the number of Paγ flares to 30 flares in 18 stars. As a typical exam-  ple of the outcome of the automatic search, we show in Fig. 6, the M3.5 star J07319+362N / BL Lyn. Both, the Paβ and the Paγ lines can be seen in broad emission exceeding 5 Å. Other spectral features in the region are still imprinted on the broad emission lines. The (about) Gaussian shape of these lines is revealed by spectral subtraction of the quiescent spectrum (see Sect. 4.2.4).
An additional visual inspection of the stars with automatic Paβ flare detections yielded another six small Paβ flares in four of the stars. Therefore, especially lower amplitude Paβ flares may be missed by the automatic detection, and we therefore screened our whole sample also by eye. Thus, we identified 12 additional stars with 16 small Paβ flares. We show an example spectrum of a visually found Paβ flare in Fig. B.1. Generally, these visually found flares are comparable in strength to the smallest flares found by the automatic detection, and the reasons why they were missed are manifold. One star, for instance, shows a large red asymmetry, which shifts the Paβ emission out Fig. 5. Pearson correlation coefficient r between pEW(Hα) and pEW(Paβ) shown as a function of µ(pEW(Hα)). T eff is colour coded as shown in the legend for stars with a significant Pearson correlation (p < 0.005). Stars with automatically or visually found flares (see Sect. 4.2.1) are marked with crosses and plusses, respectively. of the integration range. In other cases, the low number of spectra, large ranges of variability, and airglow or a combination thereof confound the search. This altogether leaves us with 32 stars showing 57 Paβ flares in comparison to 153 stars showing 357 Hα flares. Since our flare classification relies on relative variation based on the MAD of the time series, we show in Fig. 7 also the absolute deviation of the flare related pEW(Paβ) measurements compared to the median of pEW(Paβ). We caution that these values are systematically underestimated, since our pEW integration range is not broad enoug to cover the whole line during the flare. Nevertheless, flares with ∆pEW(Paβ)> 0.03 Å are typically detected automatically, while for smaller flares a non-detection by the search algorithm gets more probable. Furthermore, we note that some of the detected Paschen flares have  (Fuhrmeister et al. 2018, but no detailed discussion of the specific properties of the Paschen lines was provided there.

Statistical properties of the Paβ flares
The accumulated exposure time of all spectra considered here is 213.02 days. Comparison with the total exposure time for spectra with automatically detected Hα flares leads to a "flare  (18) 24 (15) 357 (153) a The number of stars in which these flares are detected is given in parenthesis. b In the stars with automatically Paβ flare detection Fig. 7. Deviation of the flare related pEW(Paβ) from the median of the pEW(Paβ) of the respective time series for the automatically detected flares (black dots) and the visually found flares (red diamonds). We caution that ∆pEW is systematically understimated, since the integration width is not broad enough to cover the flaring line.
duty cycle" (Hilton et al. 2010) of 2.26 %. In contrast, the flare duty cycle for automatically detected Paβ flares, with all cases of airglow contamination excluded by visual inspection, is only 0.19 %, about an order of magnitude smaller. In Fig. 8, we show the flare duty cycle as a function of spectral subtype for automatically detected Hα and Paβ flares, which increases toward later type stars for both lines. For Hα flares, such an increase was already described by Hilton et al. (2010), who found, however, lower duty cycles of 0.02 % for early M stars and 3 % for late M dwarfs in time resolved spectra of the Sloan digital sky survey. We ascribe these different numbers to the different sensitivities and flare detection methods.
Additionally, we show in Fig. 8   We also compare the flare duty cycle of all active stars, which we found to be 4.7 % for Hα and 1.0 % for Paβ, while for the inactive stars it is 2.2 % for Hα and 0.03 % for Paβ. These numbers are comparable to the values for mid-type M dwarfs and earlytype M dwarfs, since there are very few early-type M dwarfs among the active stars, while the mid-type M dwarfs have many active stars among them. Moreover, we compute the flare duty cycle only considering the stars flaring in Hα. Then the flare duty cycle becomes 4.0 % for Hα and 0.3 % for Paβ.

Stark broadening in Hα for Paβ flares
Of the spectra exhibiting Paβ emission, the vast majority shows relatively symmetric Hα line broadening. In their analysis of line asymmetries, Fuhrmeister et al. (2018) reported that red asymmetries occur frequently, blue asymmetries are more rarely observed, and symmetric line broadening is the most rarely observed variant, which is only about half as frequent as red asymmetries. Therefore, we consider a chance finding unlikely and conclude that Paβ emission is likely coupled to the occurrence of symmetric broadening. Like other authors such as Kowalski et al. (2017) and Wu et al. (2022), we consider Stark (pressure) broadening the most plausible explanation for the rather symmetric line profiles, which may alternatively be attributed to turbulent broadening or an observational time integration effect, caused by the blurring of a blue and a red asymmetry during the exposure.
Stark broadening is a consequence of high pressure in the chromosphere and, therefore, is expected to be associated with material showing larger collision rates, which lead to a larger population of higher hydrogen excitation levels. Consequently, we attribute the Paβ line emission during flares to high pressures in the chromosphere and lower transition region. Notably, flares with comparable or even higher amplitudes in Hα but no line broadening lead to Paschen line emission as exemplified by the case shown in Fig. 9, where the higher amplitude flare marked in blue does not show Hα broadening and also no enhancement of the flux in the Paschen lines. If anything, marginal excess absorption may be present in Paβ during this flare. We caution, nevertheless, that Stark broadening of the Hα line does not necessarily lead to Paschen line emission. A case in point was presented by Paulson et al. (2006), who reported on Stark broadening of the Balmer lines during a flare on the M4.0 M dwarf Barnard's star, but did not detect Paschen line emission, although they covered Paδ and higher.

Hα and Paβ line profiles during flares
Since the pEW(Paβ) is measured using a narrow integration band (see Table 1), it captures only a fraction of the flux if the line is broad. Likewise, the Hα line also exceeds the integration width used to obtain its pEW during many of the observed flares. Therefore, the pEW values do not fully characterise the strengths of broad lines.
To obtain a better understanding of the line profiles during flares, we first obtained excess spectra by subtracting the median spectrum from the flaring spectrum and, subsequently, fitted the resulting lines using Gaussians. Specifically, we used a narrow and a broad Gaussian component for the Hα line as we did in Fuhrmeister et al. (2018) and a single Gaussian component for the Paβ line. We list the best-fit parameters of our model for each flare in Table C.1 for the automatically detected flares and in Table C.2 for all flares found by visual inspection.
The two-Gaussian model reproduces the spectral line shape of Hα fairly well. In particular, the two Gaussians can account for mutual shifts of the broad and narrow component, which is a measure of the line asymmetries. For the Paβ line, we used a single Gaussian model, which is appropriate in most cases. There are eight Paβ flares, where the fit of the excess flux in Paβ was not adequate. Six of these are visually found flares, where a combination of low amplitudes and high width prevents a good fit. For some of the visually found flares also the broad component of Hα cannot be fit for similar reasons. Additionally, there are a few examples where a single Gaussian does not seem to be a suitable model for the Paβ line profile; these are marked in the Tables C.1 and C.2.
With these fits we proceeded to investigate the correlation behaviour between the Hα and Paβ line properties. We found that the strength of the narrow Hα component is not correlated with that of the Paβ line (Pearson's r = 0.23, p = 0.11). Neither is the total strength of the narrow and broad Hα components correlated with that of Paβ (r = 0.38, p = 0.007). However, the strength of the broad Hα component is correlated with that of the Paβ line (Pearson's r = 0.54 and p = 7.7 · 10 −5 ). Likewise, there are correlations between the width of the broad Hα component and that of the Paβ line (r = 0.57 and p = 2.2 · 10 −5 ) and the shift of the broad component of Hα and the line shift of Paβ (r = 0.51 and p = 0.0002). This clearly shows that the broad component of Hα and the Paβ emission are intimately related during flares and that the emission originates most probably from the same material.
From the Gaussian fits also the luminosity of Paβ L Paβ can be computed using PHOENIX photospheric models (Husser et al. 2013) for T eff and log g and the radius of the respective star . Using luminosities by Cifuentes et al. (2020) this can be converted into log L Paβ /L bol . We list these values in Tables C.1 and C.2. Values of log L Paβ /L bol range from -4.03 to -8.16 but mostly concentrating between -5.5 and -6.5. We caution that these values strongly rely on theoretical assumptions and may therefore have large systematic errors.

Asymmetries in Paβ flares
Hα often exhibits asymmetric line profiles during flares, which can manifest either as blue or red wing emission with various velocity shifts and amplitudes (for asymmetries during flares on M dwarfs see Fuhrmeister et al. (2018), and references therein; for the Sun see for example Berlicki (2007)). Asymmetric Hα line spectra during flares are usually attributed to mass mo-tions. Specifically, blue asymmetries are thought to be caused by chromospheric evaporation during flare onsets (Li et al. 2022) or prominence eruption with possible coronal mass ejections (Honda et al. 2018;Notsu et al. 2021). The origin of the red asymmetries is less certain and may, indeed, vary depending on whether the asymmetry is observed during the impulsive or the decay phase of the flare. While red asymmetries may be caused by chromospheric condensations, mainly expected to happen in the impulsive phase, in the decay phase they may be caused by coronal rain (Wu et al. 2022) or are associated with post flare loops (Namizaki et al. 2023).
Looking at the Paβ line fits, we selected all spectra showing a shift of more than 15 km s −1 (= 0.65 Å) and found two blue asymmetries and one red asymmetry among the automatically detected flares and additionally three red asymmetries among the visually found flares. Although these are small numbers, there seem to be more red than blue asymmetries, in agreement with what Fuhrmeister et al. (2018) found for Hα asymmetries. This seems to indicate again, that the shifted Paβ and Hα emission originate in the same regions. The low number of spectra with shifted Paβ emission (compared to Hα) seem to indicate, that the pressure is usually not high enough in these regions to produce a measureable amount of Paβ emission.
Out of the six detected Paβ asymmetries we discuss here the two blue asymmetries and the largest red asymmetry in more detail. While the latter belongs to J01352−072 / Barta 161 12, the two blue asymmetries occured on J01033+623 / for Paβ. For J01352−072 / Barta 161 12, we find 147.9 km s −1 for Hα and 46.7 km s −1 for Paβ. Given that the fit of very broad lines typically produces larger uncertainties on the central wavelength, we consider the velocity shifts for the first two stars to be in agreement. In the case of the red asymmetry, where the difference is about 100 km s −1 , the broad Hα profile may actually be composed of more than one component, which is not accounted for in the modelling and could, thus, explain the difference. We show all three examples in Figs. 10 and 11.
There are more example of asymmetric Hα line shapes in the spectra of the stars, which exhibit Paβ flares at some point during the time series. However, in these instances usually no Paβ emission is detectable at all, neither at the nominal wavelength nor at a shift. In these cases, the densities in the moving material are likely too low to produce Paβ emission.

Higher Paschen series lines
We visually screened the spectra with detected flares also for lines higher up in the Paschen series. For seven stars, we could find Pa 10 (at 9017.8 Å), Pa 13 (at 8667.40 Å), and Pa 14 (at 8600.75 Å) unambigiously. Pa 13, Pa 15 (at 8547.73 Å), and Pa 16 (at 8504.83 Å) are blended with the wings of the Ca ii IRT at 8500. 35,8544.44,and 8664.52 Å. Pa 17 (at 8469.59 Å) is blended with two Ti i absorption lines at 8469.474 and 8470.797 Å. These Paschen lines are therefore hard to detect, especially next to the broad and highly variable Ca ii IRT lines, which during flares usually show strong emission. We show the  Table 5.
These higher Paschen lines have the potential to trace the gas conditions in the chromosphere. A detailed chromospheric modelling with a stellar atmosphere code would yield the best results. This was done for example using PHOENIX (Hauschildt et al. 1999) by Hintz et al. (2020) for the He i infrared triplet. Such a modelling is beyond the scope of this paper and therefore we stick here to a much easier and simpler analysis using the highest observed line, which was also used by Paulson et al. (2006) for a flare on Barnard's star using hydrogen Balmer lines. The method of a pressure estimate using the highest resolved line is described in Kurochka & Maslennikova (1970). For the flares, where we detected only Pa 10 as highest line, we argue that we are sensitivity limited. Therefore, our highest resolved Paschen line is Pa 13 or Pa 14. Since for these high pressures broaden-ing by the Doppler effect plays a minor role, the Stark effect dominates the broadening, which leads to a merging of the lines. Using Equ. 9 from Kurochka & Maslennikova (1970) leads to log n e ⩽ 14.0. This value compares well with the one found by Paulson et al. (2006) using the same method. It agrees also with the electron pressure found for the higher chromosphere by detailed flare modelling for a flare on CN Leo (Fuhrmeister et al. 2010).

Consecutive Paschen line flares
There are three Paschen line flares with two consecutive spectra: one flare on YZ Cet and two flares on AU Mic. In the case of YZ Cet, the two spectra were taken about one hour apart. For the flares on AU Mic, the temporal offsets are only 7 min and 14 min. The observed evolution is diverse. For the YZ Cet flare, the combined Hα flux of the narrow and broad component decayed by only ten percent within an hour, the Paβ line flux approxi-A&A proofs: manuscript no. carmenes_pa   Paβ emission dropped by about 95 percent. Although the sample is small and the sampling sparse, we conclude tentatively that the Paβ emission during flares likely decays as fast as or faster than the Hα emission, but not significantly more slowly.

Outstanding examples of Paβ line flares
Among the identified Paβ flares, there are a number of exceptional examples. We present here the stars with the largest Paβ amplitude in our fit (see Table C.1). There are six flares with an amplitude larger than 1.0 Å belonging to four stars. All four stars have rotational periods of less than 15 days and have generally high activity levels.
The star exhibiting the flare with the overall largest amplitude is the M5.0 V star J11474+667 / 1RXS J114728.8+664405, which shows two Paβ flares, with also the second flare having a considerable amplitude. We show the two flaring spectra in Fig. 13.
The star with the second highest amplitude is the young M0.5 V star J20451−313 / AU Mic, which has three out of four automatically detected flares among the large amplitude flares. This is also the star of earliest subtype exhibiting a Paβ flare and all three large Paβ flares of this star have the three broadest Paβ lines found. We show the flaring spectra in Fig. 14. Also J22468+443 / EV Lac is an exceptional star. It has one flare among the large Paβ flares and another three automatically detected Paβ flares. Considering the three additional manually identified Paβ flares, it is the star with the largest number of Paβ flares found (followed by J20451−313 / AU Mic). We show five (out of seven) flaring spectra of various strength in Fig. B.6. The weakest of these flares shows Paβ emission, but no detectable Paγ emission while all other exhibit notable Paγ emission.
The last star with a large amplitude Paβ flare is the M4.0 V star J13536+776 / RX J1353.6+7737. The spectra are shown in Fig. B.7. J13536+776 / RX J1353.6+7737 displays a second Paβ flare, which is quite a small one and cannot be fitted properly.

Summary and conclusions
In our study we analysed the Paschen lines, which are purely of chromospheric origin, in a sample of 360 M dwarfs, which provide together more than 19 000 CARMENES spectra. We specifically used the pEW(Hα) and pEW(Paβ) to characterise the behaviour of the Paβ line in non-flaring state along the M dwarf spectral sequence. We found, that on average the Paβ line becomes more shallow for later spectral types until about spectral type M3.5; for even later spectral types the line re-deepens. Comparing the pEWs of Hα and Paβ showed, that for inactive stars with Hα in absorption in a certain T eff range, the median(pEW(Hα)) per star is anti-correlated to the median(pEW(Paβ)). Only our hottest (T eff > 4000 K) and our coolest (T eff < 3200 K) temperature interval showed no correlation. For the active stars with Hα in emission, there is in contrast only one T eff interval, where a fair correlation between median(pEW(Hα)) and median(pEW(Paβ)) could be found. Nevertheless, for time series measurements of individual stars with Hα in emission we often found correlations between pEW(Hα) and pEW(Paβ). On the other hand, for time series measurements of pEW(Hα) and pEW(Paβ) we found an anti-correlation for many stars with Hα in absorption. For both cases -looking at the median values of the stars for comparing the stellar sample and also for looking at time series of individual stars -we caution, that there are no stars with Hα in emission for T eff > 3400K.
Regarding the flaring activity of the sample stars, we found 357 Hα flares in 153 stars in comparison to 30 (57) Paβ flares in 18 (32) stars with the number in brackets including flares found only by visual inspection. Out of the automatically found Paβ flares, 86% and 69 % also show Pa γ and Pa δ in emission. Even higher Pa lines could be found unambigously up to Pa 14 for three flares (9%). The detection of even higher Pa lines is hampered by their blending with the Ca ii IRT or other stronger absorption lines. Since our pEW integration width is chosen mainly to identify flares for further characterization we applied Gaussian fitting to the Paβ line. We demonstrate the quality of the Gaussian fit of the flare excess flux density by showing some of these fits in Figs. 6,9,10,11,and 13. Both, Hα and Paβ flares are more often found in later spectral types (75% of Hα flares and 90% of Paβ flares are in stars with spectral type of M3.0 V or later). The 'flare duty cycle' (as a measure for the time fraction the star spends flaring) also increases for later spectral types, as was found for Hα already by Hilton et al. (2010). Moreover, the stars with Paβ flares nearly all show Hα in emission; only two show Hα in a transition state from absorption to emission and one star shows weak Hα absorption outside the flares. Therefore, not surprisingly, the stars with the most exceptional flares in amplitude, number and width (which we show in Figs 13, 14, B.6 and B.7) are well known very active stars as J22468+443 / EV Lac or J20451−313 / AU Mic. Additionally, we found some examples of asymmetries in the Paβ lines during flares, clearly more often associated with red asymmetries, than with blue ones.
Even more interestingly, Paβ emission during flares seems to be coupled to high densities, because almost all cases of Paβ flaring occur, when Hα exhibits symmetric broadening, which is indicative of Stark broadening (Kowalski et al. 2017) and therefore high densities. Higher amplitude flares without Stark broadening typically do not lead to Paβ emission, but Stark broadening in Hα also need not to necessarily lead to Paβ emission. As an indication of the strong coupling between the broad (Stark) component of the Hα line and the Paβ emission, we found a correlation between amplitude, width and shift of these. This sensitivity to chromospheric densities of the Pa lines deserves further investigation. For this purpose, dense time series of spectra covering a Paβ flare are needed. We identified here a number of promising candidates for such a project. These stars seem to show Paβ flares more often than the majority of M dwarfs. Such flare observations would allow to investigate, during which flare stage the Paβ emission starts and if there is a time lag to the reaction of Hα like seen for other chromospheric lines in flare studies. Together with dedicated chromospheric flare modelling this would lead to a better understanding of the density variation during a flare.
Acknowledgements. This publication was based on observations collected under the CARMENES Legacy+ project. CARMENES is an instrument at the Centro Astronómico Hispano en Andalucía (CAHA) at Calar Alto (Almería, Spain), operated jointly by the Junta de Andalucía and the Instituto de Astrofísica de Andalucía (CSIC). The authors wish to express their sincere thanks to all members of the Calar Alto staff for their expert support of the instrument and telescope operation. The authors also thank the referee for the careful reading and the suggestions for improvement. CARMENES was funded by the Max-Planck-Gesellschaft (MPG), the Consejo Superior de Investigaciones Científicas (CSIC), the Ministerio de Economía y Competitividad (MINECO) and the European Regional Development Fund (ERDF) through projects FICTS-2011-02, ICTS-2017-07-CAHA-4, and CAHA16-CE-3978, and the members of the CARMENES Consortium (Max-Planck-Institut für Astronomie, Instituto de Astrofísica de Andalucía, Landessternwarte Königstuhl, Institut de Ciències de l'Espai, Institut für Astrophysik Göttingen, Universidad Complutense de Madrid, Thüringer Landessternwarte Tautenburg, Instituto de Astrofísica de Canarias, Hamburger Sternwarte, Centro de Astrobiología and Centro Astronómico Hispano-Alemán), with additional contributions by the MINECO, the Deutsche Forschungsgemeinschaft through the Major Research Instrumentation Programme and Research Unit FOR2544 "Blue Planets around Red Stars", the Klaus Tschira Stiftung, the states of Baden-Württemberg and Niedersachsen, and by the Junta de Andalucía. We acknowledge financial support from the Agencia Estatal de Investigación (AEI/10.13039/501100011033) of the Ministerio de Ciencia e Innovación and the ERDF "A way of making Europe" through projects PID2021-125627OB-C31, PID2019-109522GB-C5[1:4], and the Centre of Excellence "Severo Ochoa" and "María de Maeztu" awards to the Instituto de Astrofísica de Canarias (CEX2019-000920-S), Instituto de Astrofísica de Andalucía (SEV-2017-0709) and Institut de Ciències de l'Espai (CEX2020-001058-M). This work was also funded by the Generalitat de Catalunya/CERCA programme and the Agencia Estatal de Investigación del Ministerio de Ciencia e Innovación (AEI-MCINN) under grant PID2019-109522GB-C53 and by the Deutsche Forschungsgemeinschaft under grant DFG SCHN 1382/2-1. This work made use of PyAstronomy , which can be downloaded at https://github.com/sczesla/PyAstronomy. A&A proofs: manuscript no. carmenes_pa

Appendix A: Comparison to PHOENIX models
Here we compare the spectral wavelength range around the Paβ line of the M0.5 V star J02222+478 / BD+47 612 and the M5.0 V star J18165+048 / G 140-051 (both also shown in Fig. 1) to PHOENIX purely photospheric models from the library of Husser et al. (2013). PHOENIX is a stellar atmosphere code (Hauschildt et al. 1999), which is widely used to compute photospheric stellar models and their synthesised stellar spectra and has been applied to CARMENES spectra to establish stellar parameters (Passegger et al. 2018;Schweitzer et al. 2019;Cifuentes et al. 2020;Marfil et al. 2021). We use a model with T eff =3900 K for the M0.5 star, a model with T eff =3200 K for the M5.0 star and log g=5.0 in both cases. Marfil et al. (2021) listed T eff = 3894±11 K and 3240±36 K and log g = 4.99±0.09 and 4.97±0.13, for the two stars, respectively. In Fig. A.1 the generally good resemblance for both spectra can be seen. Although there are some weaker lines in the PHOENIX spectra which are not seen in the observed spectra and vice versa, the stronger atomic lines match quite well. However, the Paβ line is not present in the PHOENIX spectra, making it evident, that it is a purely chromospheric line.  There are only few spectra available for the star and the variation is high, which prevents the program to automatically detect the relatively small flare.  [Å]