A comprehensive observational study of the FRB 121102 persistent radio source

We have observed the persistent radio source (QRS121102) in the VLA K-band (18 to 26 GHz) and Ku-band (12 to 18 GHz) using the C configuration. The channel width was 2 MHz and the integration time was 3 s for the K-band observations and 2 s for the Ku-band observations. The observations include six epochs from 2017 May 29 to 2017 August 10 (Table 1). In each epoch, the observation started with a single scan of the primary calibrator 3C147 (Field 0) for flux scale and bandpass calibration, and then a few cycles of the phase calibrator 1 J0555+3948 (Field 1), phase calibrator 2 J0518+3306 (Field 2), and QRS121102 (Field 3) (Table 2).

The optical spectra used in this work were obtained using Keck/LRIS. We obtained two types of observations. One, with medium spectral resolution, was used to measure the spectral width of the H$\alpha$ emission line associated with the FRB 121102 host galaxy. The other, with low spectral resolution, used two slit orientations to obtain spectra of galaxies immediately adjacent to the FRB 121102 host.

The medium resolution spectrum of the FRB 121102 host was obtained on October 12, 2018 using the 1.0” longslit and the D560 dichroic. The red side used the grating with 1200 lines/mm blazed at 7500 Å  and targeted the H$\alpha$ emission.

Low resolution spectra of the host and nearby sources were obtained on 2017 January 26 using the 1.5” long slit and the 560 dichroic. The grating used on the red side has 400 lines / mm blazed at 8500 Å, covering wavelength from about 5462 Å  to 10318 Å, with a dispersion of 1.86 Å/pixel. We estimate the spectral resolution to be $\sim 4$ Å  (1$\sigma$) from the weighted average width of the four isolated skylines (5577.0 Å, 5898.0 Å, 6315.7 Å, 7257.4 Å). On the blue side the grism has 600 lines / mm blazed at 4000 Å, covering wavelengths from 3122  Å  to 5603  Å. Standard stars were observed for the flux response calibration, arc-lamp spectra were obtained for the wavelength calibration, and bias frames were taken for the bias-subtraction. The flat field was generated using dome flats.

In this section, we describe how we measure the radio flux density of QRS121102 in each epoch.

The visibility data were calibrated and imaged using CASA following the standard procedures (task names shown in parentheses).¹¹1The Common Astronomy Software Applications (CASA) is a software package developed by the NRAO. In the calibration, the VLA antenna positions were updated (gencal). The primary calibrator 3C147 was used to find the absolute scale of the gain amplitudes by referring to the standard flux density (setjy), and to correct for the instrumental delay and the complex antenna-response variation with frequency (bandpass). The complex gain solutions (both amplitude and phase) were obtained from the phase calibrator 1, J0555+3948, that is $\sim 8^{\rm o}$ apart from QRS121102 (CASA task gaincal), and the gain amplitudes were properly-scaled using the absolute scale just obtained using 3C147 (fluxscale). All calibration solutions were then applied (applycal) to QRS121102, as well as the fainter phase calibrator 2, J0518+3306 ($\sim 2^{\rm o}$ away from QRS121102). After calibration, potential radio frequency interference (RFI) and internally generated signals were removed by flagging out narrow-banded spikes in the spectra (flagdata). In epoch 3, one scan of QRS121102 was removed as it was several times brighter than the others (within 30 minutes apart) due to strong RFI, or perhaps a passing cloud.

The calibrated visibility data were binned to $20$ s and $500$ kHz to speed up synthesis imaging using the CASA task (split). The data were gridded and Fourier transformed and the synthesized beam was deconvolved (tclean). The FWHM of the synthesized beam was $\sim 1.5^{\prime\prime}$ (major axis) $\times 1.2^{\prime\prime}$ (minor axis) in the Ku-band and $\sim 1.3^{\prime\prime}$ (major axis) $\times 0.9^{\prime\prime}$ (minor axis) in the K-band in C configuration. Our images were created using a small cell size of $0.1^{\prime\prime}\times 0.1^{\prime\prime}$ and a Briggs robust weighting of 0.5 (Briggs, 1995).

We measured the flux density of each source from the images by fitting each with a 2-D Gaussian model (imfit; Fig.1). We obtained the flux density of those point sources (based on the criterion implemented by the CASA imfit task) from the peak, and of the others from the integrated flux density within each fitting region. Uncertainties were calculated by propagation of errors of the 2D fitting model. Table 1 lists the results and Fig. 2 shows the flux density light curves in both bands.

The angular size of QRS121102 has been reported to be under $\sim$ 0.2 mas at 5 GHz as measured by the Very-long-baseline interferometry (VLBI) (Marcote et al., 2017). It is expected to be unresolved in our observations, where the beam FWHM is $\sim$ 1.4” in the Ku-band and $\sim$ 0.9” in the K-band. In our images, the 2D Gaussian fitting results show that QRS121102 is a point source in all but epoch 3, where the image is marginally resolved as $1.13\pm 0.61$ arcsec along the major axis and $0.49\pm 0.22$ arcsec along the minor axis (FWHM, deconvolvd from beam). This is likely due to the remaining phase errors in the calibration solution, as the image of phase calibrator 2 J0518+3306 appears highly distorted. The errors were removed by phase-only self-calibration in field 2, but not in field 3, for reasons outlined below.

In field 3, our observations of QRS121102 were affected by a bright source, NVSS J053153+331014, that is $\approx$ 1.8’ (twice of the primary beam FWHM in the Ku-band) away from the pointing center (QRS121102) and over one order-of-magnitude brighter ($\approx$ 3 mJy in flux density before correcting for primary-beam attenuation). To reduce the associated errors, we performed self-calibration (phase only) for NVSS J053153+331014 and QRS121102 simultaneously. We also tried to remove the flux density contribution from NVSS J053153+331014 modeled from self-calibration, and then subtract the model visibility from the corrected visibility data. Neither attempt made a significant improvement on the image of QRS121102, since self-calibration failed to correct the beam model error far away from the pointing center.

In addition, there might be remaining calibration errors since the flux densities of the two phase calibrators, J0555+3948 and J0518+3306, also vary by $\sim 4\%$ and $\sim 10\%$ throughout the epochs in each band (Fig. 2). The variations are likely not intrinsic to the source for two reasons. (1) J0555+3948 has been reported to vary by $2.0\%$ on a time scale of 251 days at 33 GHz and $3.4\%$ on a time scale of 293 days at 16 GHz (Franzen et al., 2009). It is unlikely to show a greater variability on a timescale of days, as been observed in this work (e.g. epochs 3, 4, 5). (2) More importantly, the flux density light curve of QRS121102 shows a moderate positive correlation with that of J0555+3948, and a strong positive correlation with that of J0518+3306, with a correlation coefficient of 0.67 and 0.91, respectively. To reduce potential calibration errors, in each epoch we re-scaled the flux density measurements of QRS121102 by those of the two phase calibrators (Fig. 3). We adopt the flux density scaled by J0518+3306 (phase calibrator 2) thereafter, since it has a smaller angular separation from QRS121102.

The LRIS data were processed using the LRIS automated reduction pipeline (LPipe; Perley 2019) following standard procedures: subtract the bias, make flat fields and apply flat-field correction, remove cosmic ray pixels, model and remove sky lines, perform wavelength calibration by referring to the arc lamp spectra, and determine the flux response function by observing standard reference stars. The processed 2D spectrum targets the H$\alpha$ emission.

The width of the H$\alpha$ emission line is produced by multiple effects (Rybicki & Lightman, 1979)–

• •

  dynamical velocity dispersion due to gravity,

• •

  instrumental broadening,

• •

  random thermal motions,

• •

  natural broadening,

• •

  turbulent motions.

In this section, we test the hypothesis that QRS121102 is associated with an AGN by using the velocity dispersion to roughly constrain the mass of a potential supermassive BH. The velocity dispersion is estimated from the H$\alpha$ line width.

The H$\alpha$ emission line width of the host is determined by fitting the rest-frame 2D spectrum with a 2D Gaussian function whose rotation angle is fixed at zero, plus a constant offset. We convert the observed 2D spectrum of the host into the rest-frame wavelength using a previously reported redshift of 0.19273 (Tendulkar et al., 2017). The width (1$\sigma$) of the emission line is $0.9316\pm 0.0026$ Å when combining all four exposures, and $0.8947\pm 0.0015$ Å after removing the two exposures polluted by the nearby cosmic rays (Fig. 4). A bright pixel at the center of the H$\alpha$ emission is seen in the residual (bottom panel of Fig. 4), but no extended structure that might impact the emission width result is found.

We then determine the width of the instrumental broadening effects by collapsing the 2D spectrum into 1D and fitting each of the three isolated sky lines in the field with a 1D Gaussian function plus an offset. The instrumental spectral broadening width is given by the weighted mean of the three sky lines’ 1-$\sigma$ widths ($0.638725\pm 0.000039$ Å). The natural broadening width of the H$\alpha$ line is known to be 0.46 mÅ (e.g. Rybicki & Lightman 1979). The 1-$\sigma$ width of the line is $0.6265\pm 0.0021$ Å after quadratically removing the instrumental and natural broadening effects.

The rest frame line width of a Maxwellian velocity distribution (i.e. Gaussian along the line-of-sight) is given by (Rybicki & Lightman, 1979):

Here, $\Delta\nu$ Hz is the rest-frame line width (1$\sigma$) in frequency and $\nu_{0}=4.57\times 10^{14}$ Hz is the frequency of the H$\alpha$ emission in vacuum. $\sigma_{v}$ is the velocity dispersion due to gravity, $\sqrt{2kT/m_{\rm H}}\approx 12.8$ km/s is the most probable thermal velocity assuming a gas temperature of $10^{4}$ K (e.g. Draine 2011). $v_{\rm turb}$, the turbulent velocity, is weakly constrained to be $\lesssim 10^{3}$ km/s as inferred from the scattering measurements of FRB 121102 (Table 2 of Simard & Ravi 2021). We ignored its contribution and find an upper limit of $\sigma_{v}\lesssim 30\rm~{}km~{}s^{-1}$.

Assuming that the FRB 121102 host galaxy has a central BH, we estimate its mass using the empirical M-$\sigma$ relation reported in recent literature. The BH mass is $7.8^{+8.2}_{-5.2}\times 10^{4}M_{\odot}$ using the relation derived from a sample of 88 AGN covering a stellar velocity dispersion $\sigma_{*}$ of 30 – 268 km s^-1 (Greene & Ho 2006, error bars calculated from the intrinsic scatter found in the relation). A consistent BH mass of $8.8^{16.2}_{-5.8}\times 10^{4}M_{\odot}$ is found using the relation based on 93 low-mass active galaxies (Xiao et al. 2011, $\sigma_{*}$ from 31 to 138 km/s). Other reports show similar results, though an extrapolation of the M-$\sigma$ relation is required as the sampled objects cover higher $\sigma_{*}$: the mass is $4.2^{+7.5}_{-2.7}\times 10^{4}M_{\odot}$ (error bars from the intrinsic scatter in the relation) based on a sample of 49 BH mass dynamical measurements in spiral galaxies, S0 galaxies and elliptical budges ($\sigma=$ 67–385 km s^-1; Gültekin et al. 2009), and $\sim 10^{4}M_{\odot}$ from 72 similar objects ($\sigma=$ 75 – 347 km s^-1; Ferrarese & Ford 2005; McConnell & Ma 2013). Finally, we note that the measured velocity dispersion is lower than any of those measured from a sample of 35 tidal disruption events host galaxies ($\sigma>43$ km/s) reported by Wevers et al. 2019, 2017. This suggests that the BH mass in the FRB 121102 host is lower than the BH masses of the tidal disruption event galaxy sample.

We consider two potential errors in our BH mass estimation. First, the H$\alpha$ velocity dispersion measurement may not be suitable for the dynamical analysis. The BH mass is usually estimated from the stellar velocity dispersion measured from multiple absorption lines from an optical or IR spectrum. In our observation, these stellar absorption line widths were unavailable due to the limited SNR. Instead, we infer the velocity dispersion from the line width of a single compact H$\alpha$ emission region (radius $<0^{\prime\prime}.24$ at 1$\sigma$, Kokubo et al. 2017) that is offset from the stellar continuum centroid of the galaxy by $0^{\prime\prime}.29\pm 0^{\prime\prime}.05$ (Tendulkar et al., 2017). The H$\alpha$ line width reveals the dynamics of the partially ionized warm star forming gas formed in discrete clouds. If the gas pervades in the galaxy, it is expected to have larger velocity dispersion than the stars due to turbulent motions and provide an upper limit to the BH mass. However, the H$\alpha$ region is isolated to one part of the host and may not represent the global gas dynamics in the host. This could lead to systematic errors as reported in the dynamical analyses of galaxies with irregular gas and dust distributions (Ho et al., 2002).

However, the dynamical mass implied by the H$\alpha$ velocity dispersion is comparable to the stellar mass of the FRB 121102 host galaxy inferred using the optical/IR spectral energy distribution. We assume that the system is virialized for an order-of-magnitude estimation. In an ellipsoid, the kinetic energy is dominated by random motions. The virial theorem gives a total stellar mass of $M_{*}\sim\sigma_{v}^{2}~{}R_{\rm eff}/G\sim 10^{8}M_{\odot}$, adopting a half-light-radius of $R_{\rm eff}=0.68$ kpc (Bassa et al., 2017). In a pure rotational disk (e.g. disk of a spiral galaxy), a comparable value is expected. The inferred mass is consistent with the stellar mass reported by Bassa et al. (2017) from a spectral energy distribution fit ($(1.3\pm 0.4)\times 10^{8}M_{\odot}$). This suggests that the velocity dispersion may be useful for the dynamical BH mass estimation. We also note that QRS121102 is spatially associated with the H$\alpha$ emission region.

Second, the M-$\sigma$ relations could be less reliable at our velocity dispersion for two reasons. (1) Most reports derive the empirical M-$\sigma$ relation based on a sample of more massive galaxies with larger stellar velocity dispersions. It is unclear how well the relation extends to lower velocity dispersions. (2) Even in the few reports that cover velocity dispersions down to 30 km/s, the M-$\sigma$ relation is still less reliable due to the larger uncertainty in the empirical mass-luminosity (M-L) relation at the lower end. To derive the M-$\sigma$ relation, the dynamical mass of the BHs were estimated using the virial relation, where the virial radius is either measured directly from reverberation mapping (e.g. Peterson 1993) or indirectly from the empirical M-L relation that was derived from the reverberation mapped AGN. At the low-mass end, few AGN have been reverberation mapped, so an extrapolation of the M-L relation has to be made in the BH mass estimations, introducing extra uncertainty to the resulting M-$\sigma$ relation. The M-$\sigma$ relation might be flatter at the lower-mass end based on a few intermediate mass BHs in the sample (e.g. Greene & Ho 2006), although the flattening is inconclusive without reverberation mapping results of the lower-mass BHs.

We validate our BH mass estimation result by comparing to the empirical BH-galaxy mass correlation. Reines & Volonteri (2015) measure the correlation between galaxy stellar mass and BH mass based on a sample of 341 AGN host galaxies, including a sub-sample of dwarf galaxies (Reines et al., 2013). They find that $\rm log\left(M_{BH}/M_{\odot}\right)=7.45+1.05~{}\rm log\left(M_{*}/10^{11}M_{\odot}\right)$, with a scatter of 0.55 dex. A stellar mass of $10^{8}M_{\odot}$ gives an AGN BH mass of $10^{4\sim 5}M_{\odot}$ (1$\sigma$ uncertainty assuming log-normal mass distribution), consistent with our BH mass estimation.

Data were obtained with two longslit orientations (shown in Fig. 5). The LRIS data were bias subtracted, flattened, cosmic-ray removed, skyline subtracted and flux calibrated using the LPipe (Perley, 2019). The LPipe wavelength calibration failed on the red side, since the arc lamps spectrum missed too many of the expected reference lines. We calibrated the wavelength manually by fitting the sky spectra to six isolated, bright skylines. We then selected all nearby sources for which reliable spectra could be extracted– two sources near the host on slit 1 (S1, S2) and one on slit 2 (S3) (Fig. 5).

In this section, we evaluate whether or not the host and its nearby sources belong to the same galaxy group. The redshift of each spectrum was estimated as follows. The procedures were also tested on several SDSS spectra and yield results within 1% of the known redshift values.

1. 1.

   Remove strong sky line features and data near the boundaries of the wavelength coverage.

2. 2.

   Interpolate and smooth the spectrum using a 1D Gaussian kernel ($\sigma=3$ Å) to reduce random high-frequency noise. The kernel width was chosen such that it is narrower than the spectral line width at 8000 Å ($1\sigma\sim 5$ Å) for a typical galactic velocity dispersion of $\sim 200$ km/s (e.g. Honma et al. 2012). Somewhat different kernel width choices (1Å, 3Å, 5 Å) do not change the redshift results significantly.

3. 3.

   Fit for the continuum and subtract from the spectrum. The fit was done using the Astropy package specutils, which removes spikes using a median filter and fits the spike-removed spectrum with a list of models using the Levenberg-Marquardt algorithm.

4. 4.

   Extract a 60-Å-long segment of spectrum centered at each typical galaxy absorption (or emission) line and stack them. We considered the following absorption lines as they are near the visible wavelength at a redshift of $\sim$ 0 to 1: the Balmer series, Ca K&H (3934.777, 3969.588), the G-band (4305.61), Mg (5176.7), Na (5895.6), and Ca[II] (8500.36, 8544.44, 8664.52). We ignored emission lines since none of our sources show significant emission features. The segment width of 60 Å was chosen such that it were able to include at last the $\pm 5\sigma$ region of a spectral line broadened by a typical galactic vlocity dispersion of 200 km/s (e.g. Honma et al. 2012), but not too wide to contain multiple lines.

5. 5.

   Compare the amplitude of the stacked segment at a grid of trial redshift values. In this work, the best estimation occurs at the deepest valley since we only considered absorption lines.

For S1, we find $z\approx 0.5796$ and inspect absorptive features at the wavelengths of the Ca H&K, Mg, H$\gamma$ and possibly the G-band. For S2, we find $z\approx 0.4471$ and see absorptive features at the wavelengths of the Ca K, the G-band, Na, and potentially H$\alpha$. For S3, a few weak lines (the Ca H, Mg and Na) indicate that this source might be close to S2 in redshift, but we were unable to reach a conclusion due to the low SNR. We find no significant spectral feature at the expected wavelengths assuming that S3 were at the host redshift. None of the three sources is likely to be at the same redshift as the host of FRB 121102.

We have found that QRS121102 is unlikely to be an AGN based on the low inferred BH mass ($\lesssim 10^{4\sim 5}M_{\odot}$), high radio-to-X-ray luminosity ratio and the absence of companions from the same galaxy group. In this section, we discuss other possible sources for the compact persistent radio emission.

The size of the persistent radio emission could be explained by an isolated young neutron star with luminous synchrotron emission produced in a pulsar wind nebula (PWN), or plerion. The pulsar wind forms a terminal shock at a radius where the wind momentum flux and the confining pressure reach a balance, and forms a PWN further out. The shock radius is given by: $r_{w}=\sqrt{\dot{E}/(4\pi\eta cp)}$ (Slane, 2005), where $\dot{E}$ is the pulsar energy injection rate into the wind, $\eta$ is the fraction of area covered by the wind, $c$ is the speed of light, $p$ is the confining pressure outside the shock and is proportional to the electron number density $n_{e}$ for medium with the same components. The relation gives a radius of $\sim 0.1$ pc for a canonical isolated radio pulsar, and is confirmed to be $\approx 3\times 10^{17}$ cm from the X-ray images of the Crab nebula (Weisskopf et al., 2000). For the Crab pulsar, the spin down energy rate is $\dot{E}\sim 4.5\times 10^{38}$ ergs s^-1 (Staelin & Reifenstein, 1968), and the confining pressure outside the shock is proportional to the density of the medium, which can be approximated as a typical ISM ($n_{e}\sim 10^{-1}$ to $10^{-2}$ cm^-3; e.g. Draine 2011). In comparison, the energy ejection rate of FRB 121102 by the flares and the wind into the surrounding medium is estimated to be $\dot{E}\sim 10^{39}$ to $10^{40}$ ergs s^-1 (Fig. 5 of Margalit et al. 2020), and the medium ahead of the termination shock is likely denser than a typical ISM, as indicated from the high RM of FRB 121102 (Michilli et al., 2018).

One example of the PWN emission model that produce the observed radio luminosities and the size of QRS121102 is presented in Margalit & Metzger (2018), where the authors explain the persistent emission using a concordance FRB model. On a large scale ($\sim 10^{15}$ cm, Eqn. 4 of Metzger et al. 2019), the train of ion-electron shells merge into a steady wind and feed into a nebula via a terminal shock, which heats up electrons in the nebula and produces the persistent synchrotron radio source. Based on the lack of self-absorption feature down to 6.0 GHz in the spectrum of the persistent source (Chatterjee et al., 2017), Margalit & Metzger (2018) estimate that the emission region’s radius $R\gtrsim 0.46\times 10^{17}~{}\rm cm$ adopting the observed luminosity at 6 GHz. Moreover, Resmi et al. (2020) report the lack of self-absorption down to 400 MHz in their GMRT observations. Using their flux density measurement at 400 MHz and the scaling relation $R\propto L_{\nu,\rm obs}^{4/11}\nu_{\rm obs}^{-10/11}$ (Eqn. 21 in Margalit & Metzger 2018), we find that $R\gtrsim 5.4\times 10^{17}~{}\rm cm$. This emission size is consistent with the constraint implied by the lack of refractive scintillation in the Ku-band.

Among the more unique features of QRS121102 is its unusually flat radio spectrum at GHz frequencies (Chatterjee et al., 2017). As we have noted, several examples of compact radio sources of similar luminosities exist in dwarf galaxies. Indeed, two recent discoveries of transient radio sources not associated with BHs also reached similar radio luminosities. The first, FIRST J141918.9+394036, had a peak radio luminosity of $2\times 10^{29}\rm~{}erg~{}d^{-1}~{}Hz^{-1}$ at 1.4 GHz, and is hosted by a star-forming dwarf galaxy (Law et al., 2018). FIRST J141918.9+394036 is most likely the afterglow of an off-axis long-duration GRB (Mooley et al., 2021). The second, VT J121001+495647, had a peak radio luminosity of $1.5\times 10^{29}\rm~{}erg~{}s^{-1}~{}Hz^{-1}$ at 5 GHz, was associated with a star-forming region, and was ascribed to interaction with a dense circum-stellar medium ejected through binary interaction (Dong et al., 2021). In both cases, however, classical synchrotron spectral shapes were observed together with secular time-evolution, unlike in the case of QRS121102. We urge continued wideband monitoring of QRS121102, together with more detailed evaluation of potential empirical analogs.

Finally, we rule out a few other possible origins of QRS121102 based on our observations. First, the source cannot be a supernova remnant (SNR) due to its high luminosity. We have calculated the specific luminosity of QSR121102 at 3 GHz and 10 GHz as $L_{\nu}\sim 10^{22}$ W Hz^-1. Varenius et al. (2019) and Ulvestad (2009) recently report the radio luminosities (5 GHz and 8.4 GHz) of 102 SNRs in the merging galaxies Arp 220 and Arp 229. The brightest SNR in the sample is $L_{\nu}<10^{21}$ W Hz^-1, over one order of magnitude lower than that of QSR121102. Moreover, the luminosity is inconsistent with the star formation rate (SFR) of the host galaxy if the source were SNR(s). The brightest SNR and the SFR of a galaxy have been reported to be related (Chomiuk & Wilcots, 2009) as: $L_{1.4}^{max}=95^{+31}_{-23}\rm~{}SFR^{0.98\pm 0.12}$, where $L_{1.4}^{\rm max}$ is in a unit of $10^{24}$ erg s^-1 Hz^-1 and SFR is in $M_{\odot}$ year^-1. Adopting the SFR upper limit of 0.4 $M_{\odot}$ year^-1 based on the host H$\alpha$ emission flux (Tendulkar et al., 2017), the brightest SNR in the host would be $L_{1.4}^{\rm max}\approx 10^{18.6}$ W Hz^-1, over 3 orders of magnitude lower than that of QSR121102. The persistent source is too bright to be an SNR. Second, the source is too bright for a long-duration GRB (LGRB) radio afterglow. Adopting a typical LGRB peak radio luminosity of $L_{\nu,\rm 8.5GHz}\sim 2\times 10^{31}$ erg^-1 s^-1 Hz^-1 and a decay rate of $\propto t^{-2}$ from day 10 after the GRB (e.g. Berger et al. 2003), the radio luminosity would have reached the level of QRS121102 ($L_{\nu,\rm 10GHz}\approx 2\times 10^{29}$ erg^-1 s^-1 Hz^-1) within 3 weeks, while the radio luminosity of QRS121102 have been nearly constant below 10 GHz for years (e.g. Resmi et al. 2020; Tendulkar et al. 2017).