The rare extended radio-loud narrow-line Seyfert 1 galaxy SDSS J1030+5516 at high resolution

Recently, Rakshit et al. (Astrophys. J. 869:173, 2018) reported the discovery of SDSS J103024.95+ 551622.7, a radio-loud narrow-line Seyfert 1 galaxy having a ∼100\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}${\sim}100$\end{document} kpc scale double-lobed radio structure. Here we analyse archival radio interferometric imaging data taken with the Very Large Array (VLA) at 5 GHz, and with the Very Long Baseline Array (VLBA) at 4.3 and 7.6 GHz. Two hotspots and a compact core are detected with the VLA at arcsec scale, while a single milliarcsec-scale compact radio core is seen with the highest resolution VLBA observations. The Fermi Large Area Telescope did not detect γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\gamma $\end{document}-ray emission at the position of this source. In the mid-infrared, the Wide-field Infrared Survey Explorer satellite light curve, covering more than 7 years and including the most recent data points, hints on flux density variability at 3.4 μm. Our findings support the notion that this source is a young version of Fanaroff–Riley type II radio galaxies.

by their narrow permitted optical lines (the full width at half maximum, FWHM, of Hβ line is below 2000 km s −1 , Goodrich 1989), a flux ratio of [O III]λ5007 to Hβ smaller than 3 (Osterbrock and Pogge 1985), and the strong emission of the Fe II multiplets. However, Cracco et al. (2016) showed that having strong iron lines may not be a distinctive property of NLS1 sources in accordance with the study of quasar emission lines of Boroson and Green (1992), which showed the anticorrelation between the strength of Fe II and O III lines.
The narrow permitted lines of NLS1 sources are explained with their relatively lower-mass central black holes, 10 6 -10 8 M (Mathur 2000), which consequently means high accretion rates close to the Eddington limit (Boroson and Green 1992). Based upon these, Mathur (2000) proposed that NLS1 sources can be young AGN residing in rejuvenated galaxies. Alternatively, the narrow Hβ lines of NLS1 sources can be due to orientation effect, if their disklike broad-line regions are seen pole-on (Decarli et al. 2008).
Similarly to AGN in general, a small fraction, ∼7% of NLS1 sources are radio-loud (Komossa et al. 2006;Zhou et al. 2006), where radio-loudness is determined by the ratio of the 6 cm radio to the 4400 Å optical flux density following Kellermann et al. (1989). Singh and Chand (2018) studied the radio properties of a large sample of opticallyselected NLS1 sources, and found that the radio-detected ones have small sizes, <30 kpc. The most radio-loud NLS1 sources (and the ones with the highest radio luminosity investigated by Singh and Chand 2018) show blazar-like properties: flat radio spectrum, compact radio cores, high brightness temperatures, significant variability, flat X-ray spectra, and altogether blazar-like spectral energy distribution (SED; e.g. Yuan et al. 2008). Several of them were also detected in γ -rays with the Fermi satellite (for a full list, see Romano et al. 2018), and in a few of them superluminally moving radio jet components were imaged with very long baseline interferometry (VLBI) technique (Lister et al. 2016). Therefore these sources, similarly to blazars, are thought to possess relativistic radio jets inclined at small angle to the line of sight.
Few of the radio-loud NLS1 sources have kpc-scale (from a few tens of kpc to ∼100 kpc) radio structures. Doi et al. (2012) found that the detection rate of extended radio emission in NLS1 sources is lower than in broad-line AGN. This is confirmed more recently by Berton et al. (2018). They observed 74 NLS1 sources and found that the majority of flat-spectrum radio-loud NLS1 sources have compact morphology on kpc scale. In most of the extended radioloud NLS1 sources, the radio emission is two-sided (Doi et al. 2012;Richards and Lister 2015;Congiu et al. 2017;Gabányi et al. 2018b).
Recently, Rakshit et al. (2018) reported the discovery of a radio-loud NLS1 source, SDSS J103024.95+551622.7 (hereafter J1030+5516) with arcsec-scale structure similar to those of Fanaroff-Riley type II radio galaxies (FR II, Fanaroff and Riley 1974). The projected linear size, ∼110 kpc is among the largest values in radio-loud NLS1 sources. Using low-resolution (∼5 ) radio data of J1030+5516 from the Faint Images of the Radio Sky at Twenty-Centimeters (FIRST) survey (Becker et al. 1995), Rakshit et al. (2018) argue that the inclination angle of the jet in the core region is <12 • with respect to the line of sight.
Here we present sub-arcsec resolution archival Very Large Array (VLA) A-configuration data, and milli-arcsec (mas) resolution VLBI data of J1030+5516, which support the claims of Rakshit et al. (2018). Additionally, we analysed more than 10 yr of Fermi Large Area Telescope (LAT; Acero et al. 2015, and references therein) data to constrain the high-energy properties of the source. We also reevaluated the mid-infrared light curves covering more than 7 yr, obtained with the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010) satellite.
In the following, we assume a flat ΛCDM cosmological model with H 0 = 70 km s −1 Mpc −1 and Ω m = 0.27. At the redshift of J1030+5516, z = 0.435, 1 angular size corresponds to a projected linear size of 5.65 kpc.
2 Observing data 2.1 Archival VLA radio data J1030+5516 was observed at 5 GHz with the VLA in its most extended A configuration on 1992 October 20 (project code: AF233). The raw data were obtained from the US National Radio Astronomy Observatory (NRAO) archive. 1 The 1 https://archive.nrao.edu/. on-source integration time was 1 min, the total bandwidth was 100 MHz. Phases and amplitudes were calibrated in the NRAO Astronomical Image Processing System (AIPS, Greisen 2003) in a standard way. The flux density scale was set using the amplitude calibrator source 3C286. We used the DIFMAP (Shepherd et al. 1994) software for imaging and for fitting Gaussian brightness distribution model components directly to the interferometer visibility data.

Archival VLBA radio data
J1030+5516 was observed with the Very Long Baseline Array (VLBA) on 2016 May 28 at 4.3 and 7.6 GHz (project code: BP192, PI: L. Petrov) in the framework of the widefield VLBA calibrator survey (L. Petrov 2019, in preparation). Nine (Brewster, Fort Davis, Hancock, Kitt Peak, Los Alamos, North Liberty, Owens Valley, Pie Town, St. Croix) and eight (all the above but Hancock) antennas of the array were used at the lower and higher frequency, respectively. The bandwidth was 256 MHz and the integration time was nearly 1 min at both frequencies. The calibrated visibilities were obtained from the Astrogeo website. 2

Archival Fermi/LAT γ -ray data
We analysed the archival Fermi/LAT data of J1030+5516. We looked for γ -ray signals in the available 10.3-year data which cover the time range between 2008 August 4 and 2018 November 26. We derived Fermi/LAT upper limits using the routines included in the fermipy package (Wood et al. 2017). We selected a 15 • circular region around the position of J1030+5516, and an energy range of 0.1-100 GeV. We used Pass 8, SOURCE type photons, with P8R2_SOURCE_V6 responses. The fit included a diffuse galactic foreground (gll_iem_v6) and an isotropic component (iso_P8R2_SOURCE_V6_v06).

WISE data
The mid-infrared WISE satellite scanned the whole sky in four bands at 3.4, 4.6, 12, and 22 µm (referred as W1, W2, W3, and W4) during its original mission phase in 2010 (Wright et al. 2010). Afterwards, the satellite measurements are continued within the framework of the NEOWISE (Near-Earth Object WISE) project (Mainzer et al. 2014). After four months of NEOWISE observations, the satellite was hibernated for 34 months. Then the NEOWISE Reactivation Mission continued. In this currently on-going phase, observations are conducted only at the two shorter wavelength bands, since the cooling material required for W3 and W4 receivers has been depleted. The WISE satellite observes the same regions of the sky in every ∼180 days.
We downloaded the WISE single exposure data 3 up until 2017 November, and followed the procedure as in Gabányi et al. (2018a). We adapted the guidelines in the Explanatory Supplement Series 4 to discard bad quality data points. None of the measurements were affected by the South-Atlantic Anomaly, or scattered light from the Moon. We used only those measurements for which the frame image quality score ('qi_fact') was 1.0, since values less than 1.0 mark data where residual light system motion may degrade the flux measurements. Ten per cent of data were discarded because of this effect. The contamination and confusion flags ('cc_flags') did not indicate problems for any of the measurements.
The final light curve contains 177 points grouped into ten mission phases both in W1 and W2 bands. Each mission phase lasted usually for 2 days and contains 12-16 points, except for one mission phase, which is the combination of two ∼1-day long observations separated by 2 days and contains 26 data points.

VLA data
The 5-GHz VLA-A observation revealed three distinct features, the core, the southwest (SW) and the northeast (NE) lobes ( Fig. 1). This structure is in good agreement with the lower-resolution FIRST image presented by Rakshit et al. (2018). Three circular Gaussian model components can adequately describe the radio structure. Their parameters are given in Table 1. The distance between the components fitted to the hotspots in the SW and NE lobes is 20. 43 ± 0. 01, which agrees well with the 20. 5 source size at 1.4 GHz, as derived by Rakshit et al. (2018) from the FIRST image.
Using the overall spectral index value α = −0.65 (defined as S ∝ ν α , where S is the flux density and ν the frequency) and the integrated flux density S 1.4 = 155 mJy derived from the 1.4-GHz FIRST data by Rakshit et al. (2018), the flux density at 5 GHz can be expected as ∼ 68 mJy. The sum of flux densities of the fitted components (Table 1), S 5 = 27.8 mJy, is well below this value, indicating that a significant amount of diffuse radio emission was resolved out in the 5-GHz VLA-A observation.

VLBA data
At both frequencies, a single unresolved component was detected with the VLBA. There was no additional radioemitting feature down to 0.8 mJy beam −1 within the undistorted field of view with a radius of 0. 3 at 4.3 GHz, and  Table 1 are labeled down to 1.1 mJy beam −1 within the undistorted field of view with a radius of 0. 2 at 7.6 GHz. We used the DIFMAP (Shepherd et al. 1994) software to fit the visibilities with a circular Gaussian brightness distribution. At 4.3 GHz, the radio emission can be best described by a point source model. During model-fitting using a Gaussian brightness distribution as a starting model, the FWHM size of the feature converged to an unrealistically small value (∼10 −6 mas), indicating the unresolved nature of the detected component. The flux density of the point source model component is S 4.3 = 10.6 ± 0.3 mJy.
At 7.6 GHz, a stable fit could be reached with using a single circular Gaussian component. Its parameters are: flux density S 7.6 = 15.2 ± 0.5 mJy, FWHM diameter θ = 0.31 ± 0.09 mas. However, this size is still smaller than the minimum resolvable angular size of the interferometer array, 0.4 mas, calculated following the formula of Kovalev et al. (2005). Therefore J1030+5516 remained unresolved with VLBI at 7.6 GHz as well.
Using the 0.4 mas upper limit to the source size, we can calculate a lower limit to the brightness temperature: where z is the redshift, S is the flux density in Jy, ν is the observing frequency in GHz, and θ is the FWHM size of the Gaussian component in mas. Thus, the lower limit to the brightness temperature of J1030+5516 is T B ≥ 3 × 10 9 K. This value is below the equipartition brightness temperature limit (∼5 × 10 10 K; Readhead 1994) by an order of magnitude. However, since it is a lower limit only, it does not exclude the possibility of Doppler boosting caused by relativistic beaming in the jet. More sensitive and higherresolution VLBI data would be needed to place tighter constraints on the angular size of the compact central source in J1030+5516 and thus on the brightness temperature. Using these simultaneous dual-frequency VLBA observations, we can derive the two-point spectral index α 7.6 4.3 = 0.6 ± 0.2 of the core, indicating an inverted radio spectrum of the most compact component. This is not in contrast with the spectral index derived by Rakshit et al. (2018) (α = −0.65 ± 0.04), since that value was obtained using low-resolution radio observations. Those measure the total flux density of J1030+5516 which is dominated by the steep-spectrum lobes and diffuse emission, completely resolved out on the long baselines of the VLBA.

Fermi/LAT data
No γ -ray emitting source was detected in the available Fermi/LAT data at the position of J1030+5516. The upper limits of the γ -ray fluxes in six energy ranges are given in Table 2.

WISE mid-infrared lightcurve
We converted the WISE magnitudes to flux densities, following the description in the Explanatory Supplement Se- ries. 5 To investigate the variability, we calculated the reduced χ 2 for all the data in a given band, and for each mission phase, separately. For every single mission phase, the reduced χ 2 2, indicating no variability on a few-day long time scale. The reduced χ 2 values for all measurements are 3.3 at 3.4 µm and 2.8 at 4.6 µm, showing a hint of variability at the shorter wavelength. Rakshit et al. (2018) also analysed the infrared light curve of J1030+5516 measured with the WISE satellite. At the time of their publication only eight mission phases were available, until 2016 November. They found that the object is not variable at 3.4 µm and 4.6 µm. We also calculated the reduced χ 2 using the first eight mission phases used by Rakshit et al. (2018). We found that the reduced χ 2 are lower for both bands, 3.1 at 3.4 µm and 2.4 at 4.6 µm.
In Fig. 2, we plot the weighted average flux densities and the standard deviations for each mission phase in both bands.

Discussion
The obtained lower limit of the brightness temperature of the mas-scale radio emitting core of J1030+5516 agrees well with the values derived by Gu et al. (2015) for a sample of 14 radio-loud NLS1 sources, 10 8.4 K < T B < 10 11.4 K. These low values compared to powerful blazar jets are explained by intrinsically low jet power by Gu et al. (2015). The brightness temperature measured in J1030+5516 is also similar to another NLS1 source with a radio structure extended to ∼150 kpc, SDSS J110006.07+442144.3 (Gabányi et al. 2018b). Richards and Lister (2015) investigated three radio-loud NLS1 galaxies with kpc-scale radio structures and found mildly relativistic jets. Comparing the mas-scale structures, J1030+5516 is more compact and its brightness temperature is 3-20 times larger than the three NLS1 sources studied by Richards and Lister (2015).
The flux density of the core component measured at 5 GHz with the VLA at arcsec scale (9.3 ± 0.4 mJy) is below the value measured at mas-scale resolution with the VLBA at a slightly lower frequency of 4.3 GHz (10.6 ± 0.3 mJy). As this cannot be explained by resolution effect, it is more likely related to variability of the radio flux density.
In the standard picture of expanding radio galaxies, the lobe located closer to the observer is seen farther away from the host galaxy in projection (Longair and Riley 1979), thus the arm-length ratio of the brighter to the fainter lobe is larger than one. In J1030+5516, the distance between the brighter NE lobe and the central compact radio feature, C is d NE = 13. 3 ± 0. 1, while the fainter SW lobe is at d SW = 7. 35 ± 0. 04. Therefore, it is most likely that the NE feature is at the approaching side of the source.
The arm-length ratio of the approaching to the receding lobes can be used to estimate the inclination angle (i) of the source, assuming there is no significant difference in the medium surrounding the jets on the two sides of the host galaxy. Using the equation of Taylor and Vermeulen (1997), the arm-length ratio can be given as From the VLA data of J1030+5516, the arm-length ratio is 1.81, thus β cos i = 0.29, which gives a lower limit for the jet speed β > 0.3, and an upper limit for the inclination angle i < 73 • . Similar equation describes the flux density ratio of the jet and counter-jet side (Taylor and Vermeulen 1997): where k equals 2 for continuous jet, and 3 for discrete jet components. Rakshit et al. (2018) using the FIRST data obtained a flux density ratio of 4.35, from which using the spectral index α = −0.65 they derived β > 0.2, and i < 79 • , which agree with our values obtained from the armlength ratio. The higher resolution 5-GHz VLA data give a flux density ratio of the approaching and receding lobe of 2.85 ± 1.8, which agrees within the errors with the value Rakshit et al. (2018) used. (Since we have no information on the spectral index of these radio features, we did not use it further to calculate the β cos i value.) The Doppler factor is defined as, where γ = 1/ 1 − β 2 is the Lorentz factor. Using the value β cos i = 0.29, and the lower limit on β, 0.29, one can obtain an upper limit on the Doppler factor, δ < 1.34. Doi et al. (2012) estimated the Doppler factor of a few radio-loud NLS1 sources by comparing the observed and intrinsic core powers. Fitting the VLA data we obtained the core flux density at 5 GHz. Assuming the spectral index of the core derived from the VLBA data (α 7.6 4.3 = 0.6), the observed 5-GHz power of the core is 3.6 · 10 24 W Hz −1 . The intrinsic core power can be estimated using the empirical correlation found for radio galaxies, log P core 5 GHz = (0.62 ± 0.04) log P total 408 MHz + (7.6 ± 1.1), where P total 408 MHz is the source's total power at 408 MHz (Giovannini et al. 2001). The closest frequency where the flux density of J1030+5516 was measured is 365 MHz within the framework of the Texas survey (Douglas et al. 1996). We used that value (0.474 ± 0.038 Jy) and the spectral index of the whole source derived by Rakshit et al. (2018) to calculate P total 408 MHz . The corresponding intrinsic 5-GHz core power is ∼ 10 24 W Hz −1 . Following Doi et al. (2012), if the difference between the observed and intrinsic core power is caused by relativistic beaming, their ratio can be given as δ 3−α . This implies a δ = 1.7. However, the uncertainty of the estimation of intrinsic core power allows for lower δ values and as high as 12. Rakshit et al. (2018) also used the above argument to estimate a Doppler factor. Instead of the core power, they used the intrinsic and observed core dominance parameter and obtained δ = 3.3. The derived δ values do not contradict the result of the VLBA observation which gave a lower limit of the brightness temperature. On the other hand, δ > 1.34 cannot be accommodated with the β cos i = 0.29 derived from the kpc-scale structure. Thus, if we accept the higher Doppler factors derived from the power of the core component, either the jet direction, i, or the jet speed, β, or both of them change significantly between the kpc and pc scale. The fact that only a single, compact radio-emitting feature was detected at mas-scale resolution indicates that the jet is not oriented close to the plane of the sky at pc scales.
The 1.4-GHz radio power of the source calculated from the flux density detected in the FIRST survey (∼150 mJy) is P total 1.4 GHz = 9 · 10 25 W Hz −1 . An and Baan (2012) studied the evolutionary sequence of symmetric extragalactic radio sources. On their radio power versus projected linear size diagram, J1030+5516 is among the large symmetric objects and the low-power FR II radio galaxies.
The projected linear size of the source is D ∼ 115 kpc. Assuming a constant expansion velocity and using the limit β cos i = 0.29 derived from the kpc-scale radio structure, the kinematic age of the source can be estimated as where D 0 is the full (deprojected) size of the radio source.
Using the upper limit on the inclination angle (73 • ), the lower limit on the age of the kpc-scale radio structure is 2 · 10 5 yr. The inclination angle of J1030+5516 should be below ∼32 • to obtain an age 10 6 yr and its age would reach 10 7 yr, if the inclination angle would be 5 • , which is not consistent with a nearly symmetric kpc-scale structure. Thus, J1030+5516 seems to be younger than typical FRII radio galaxies whose lifetimes are estimated to be (10 6 -10 7 ) yr (O'Dea et al. 2009). Rakshit et al. (2019) investigated the infrared properties of 520 NLS1 sources using WISE data. They found that more than 50% of the sources classified as variable in the AllWISE Source Catalog ) fall within a specific region in the WISE color-color diagram, the 'WISE Gamma-ray strip' (WGS). The WGS is defined by D' Abrusco et al. (2012) and Massaro et al. (2012) as a distinct region where blazars are located. The infrared colors of J1030+5516 place this source within the WGS as well. Based upon this and its radio-loudness, D' Abrusco et al. (2014) included J1030+5516 in the list of γ -ray emitting blazar candidates. They defined three classes of blazar candidate sources based on their decreasing likeliness of being blazars. J1030+5516 fell into the second class. This motivated us to analyse the available Fermi/LAT data to look for evidence of γ -ray emission. However, no γ -ray emission was found at the position of J1030+5516. Rakshit et al. (2018) reported that J1030+5516 did not show variability in infrared at the wavelengths 3.4 µm and 4.6 µm measured with the WISE satellite. Since their publication, data from two additional epochs of WISE measurements became public. Our analysis of ten mission phases of the WISE data indicate a slight hint of long-term variability at the shorter wavelength due to the two additional epochs. There is no sign of short time scale (few day long) flux density changes neither in W1 nor in W2 bands.
Contrary to the finding of Rakshit et al. (2018), Graham et al. (2015) detected optical variability in J1030+5516. Graham et al. (2015) studied the optical light curves measured by the Catalina Real-time Transient Survey (Drake et al. 2009) to look for periodic variability in quasars, which they interpreted as induced by a closely separated binary supermassive black hole in those sources. They list J1030+5516 among the binary candidates with a period of 1515 days, a separation of 0.006 pc and a rest-frame merger time of 2.2 · 10 5 yr (assuming a mass ratio of 0.5). However, Vaughan et al. (2016) called for cautious approach when only a few cycles are used to assess periodic variability. In any case, whether periodic or not, J1030+5516 seems to show some variability in optical and infrared bands on time scales of years, as suggested by Graham et al. (2015) and perhaps also by our light curve compiled from WISE data (Fig. 2).

Summary and conclusion
To reveal its structural properties at high resolution, we analysed archival radio interferometric observations of J1030+5516, a recently discovered rare radio-loud NLS1 source (Rakshit et al. 2018) having kpc-scale extended structure reminiscent of FR II radio galaxies. The compact central core and the nearly symmetric hotspots on its two sides seen in our 5-GHz VLA A-configuration image (Fig. 1) are consistent with the lower-resolution 1.4-GHz VLA FIRST survey image. A large fraction, more than a half of the radio emission originates from diffuse structures related to the two lobes, and is resolved out with the VLA at 5 GHz. At even higher, mas-scale resolution, the 4.6 and 7.6 GHz VLBA observations show a single, unresolved core feature, with a brightness temperature T B > 3 · 10 9 K. The observed high-resolution radio morphology and the derived parameters of J1030+5516 make it similar to the few known NLS1 sources with large-scale extended radio emission (Richards and Lister 2015;Gabányi et al. 2018b).
Using the arm-length ratio, we derived a lower limit for the kpc-scale jet speed (β > 0.3), and a corresponding upper limit of the inclination angle, i < 73 • , which are slightly more stringent than the ones given by Rakshit et al. (2018). The lack of complex radio morphology at mas scale indicates that the inner jet of J1030+5516 is not oriented close to the plane of the sky at pc scales.
The WISE infrared colors of the source place it in the WGS. However, our analysis of all available Fermi/LAT data did not reveal any γ -ray emission at the position of J1030+5516. There is indication of at least some variability on a multi-year time scale, both in the mid-infrared based on the most complete WISE light curve (Fig. 2), and in the optical according to monitoring observations analysed by Graham et al. (2015). Based on two epochs of radio observations, we also found indication of flux density variability.
Our results support the conclusion of Rakshit et al. (2018) that the NLS1 galaxy J1030+5516 is likely a low-power and young version of an FR II radio galaxy with a double-lobed structure. There is no indication of relativistic boosting in J1030+5516, neither from radio imaging, nor from variability or γ -ray data, however it cannot be completely ruled out either.
Acknowledgements Open access funding provided by MTA Research Centre for Astronomy and Earth Sciences (MTA CSFK). This project received support from the Hungarian Research, Development and Innovation Office (OTKA NN110333). K.É.G. was supported by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences. P.V. acknowledges support from Fermi grant NNM11AA01A. We used in our work VLBA data from project BP192F7 provided by the National Radio Astronomy Observatory that is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. We thank Leonid Petrov for providing results of data analysis prior to publication. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

Conflict of Interest
The authors declare that they have no conflict of interest.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creative commons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.