The first spectroscopic IR reverberation programme on Mrk 509

Near IR spectroscopic reverberation of Active Galactic Nuclei (AGN) potentially allows the infrared (IR) broad line region (BLR) to be reverberated alongside the disc and dust continua, while the spectra can also reveal details of dust astro-chemistry. Here, we describe results of a short pilot study (17 near-IR spectra over a 183 d period) for Mrk 509. The spectra give a luminosity-weighted dust radius of $\langle R_{\mathrm{d,lum}} \rangle = 186 \pm 4$ light-days for blackbody (large grain dust), consistent with previous (photometric) reverberation campaigns, whereas carbon and silicate dust give much larger radii. We develop a method of calibrating spectral data in objects where the narrow lines are extended beyond the slit width. We demonstrate this by showing our resultant photometric band light curves are consistent with previous results, with a hot dust lag at>40 d in the K band, clearly different from the accretion disc response at<20 d in the z band. We place this limit of 40 d by demonstrating clearly that the modest variability that we do detect in the H and K band does not reverberate on time-scales of less than 40 d. We also extract the Pa$\beta$ line light curve, and find a lag which is consistent with the optical BLR H$\beta$ line of $\sim$70-90 d. This is important as direct imaging of the near-IR BLR is now possible in a few objects, so we need to understand its relation to the better studied optical BLR.


INTRODUCTION
Active Galactic Nuclei (AGN) are powered by mass accretion onto a central Supermassive Black Hole (SMBH), and are the most luminous compact objects in the Universe.They exhibit a wide range of emission characteristics resulting in an AGN zoo of differing classifications.One of the fundamental distinctions between categories ★ E-mail: jake.a.mitchell@durham.ac.uk † Visiting Astronomer at the Infrared Telescope Facility, which is operated by the University of Hawaii under contract NNH14CK55B with the National Aeronautics and Space Administration.
is that which separates AGN with both narrow and broad emission lines seen in their optical spectra (Type 1) from those which only display narrow lines (Type 2).Historically unified schemes proposed a basic AGN anatomy consisting of a central SMBH, a thin accretion disc, separate compact broad and extended narrow line regions and an outer dusty toroidal structure, with each of these components existing on different physical size scales.It has been argued, through these unified schemes, that orientation effects resulting in the obscuration of the Broad Line Region (BLR) for certain viewing angles is the physical mechanism separating AGN into type 1 and type 2 sources (Lawrence 1987;Antonucci 1993;Urry & Padovani 1995;Netzer 2015).Originally the dichotomy between type 1 and 2 AGN was solely explained by these orientation effects.However the predominant view now is more holistic, with differences in emission characteristics attributed to numerous factors, such as mass accretion rate and black hole mass in addition to orientation (Jin et al. 2012;Kubota & Done 2018).Despite this slight paradigm shift, much evidence for the presence of a large dusty obscuring structure, still commonly referred to as the torus for ease, remains.Strong infrared continuum emission at  > 1 m, thought to be reprocessed optical/UV emission from the accretion disc, and observations of broad emission lines in the polarised spectra of some type-2 AGN (Antonucci & Miller 1985;Antonucci 1993;Cohen et al. 1999;Lumsden et al. 2001;Tran 2003) are prime examples of evidence for the torus.For a small subset of AGN, constraints on the spatial extent of the torus have been placed at (0.1-6) pc through near-IR, mid-IR (Tristram et al. 2009;Pott et al. 2010;Kishimoto et al. 2011;Burtscher et al. 2013;Kishimoto et al. 2013) and sub-mm interferometric observations (García-Burillo et al. 2016;Imanishi et al. 2018) along with high angular resolution mid-IR imaging (Packham et al. 2005;Ramos Almeida et al. 2011a).These findings are roughly consistent with a powerlaw size-luminosity relation  dust =  0.5 (Kishimoto et al. 2009(Kishimoto et al. , 2011)).Modelling the torus has also provided a useful probe of this structure.Through the fitting of infrared SED's, using the CLUMPY torus model (Alonso-Herrero et al. 2011), and spectroscopy, tight constraints have been placed on a range of torus properties.This and other studies indicate that the torus is of a clumpy geometry (Ramos Almeida et al. 2011b;Ichikawa et al. 2015).However for the general case the precise structure and geometry of the torus remains unknown, largely due to the size scales in question, rendering direct imaging impossible with current technology for the vast majority of sources.
Where interferomtery is not possible, the inner dusty structures of AGN can be probed with reverberation mapping techniques, which are predominantly used in mapping the spatial extent of the broad line region in AGN, but can also be used to measure the physical size of other components such as the accretion disc or hot dust (Clavel et al. 1989;Glass 1992;Sitko et al. 1993;Nelson 1996;Oknyanskĳ & Horne 2001;Glass 2004).This method assumes that the variable near-IR continuum arises from hot dust reprocessing of the variable optical/UV emission from the accretion disc and therefore the inner edge of the torus can be measured as the time lag between these two variable signals.This method is well tested, with 20 AGNs having been reverberation mapped in the near-IR through long term coordinated optical and near-IR monitoring campaigns (Clavel et al. 1989;Sitko et al. 1993;Nelson 1996;Oknyanskĳ & Horne 2001;Glass 2004;Minezaki et al. 2004;Suganuma et al. 2006;Koshida et al. 2014;Vazquez et al. 2015).Dust radii measured in this way are systematically smaller than radii found with near-IR interferometry (Oknyanskĳ & Horne 2001;Kishimoto et al. 2007;Nenkova et al. 2008).The same discrepancy is seen between reverberated radii and those determined by the dust sublimation temperature and AGN bolometric luminosity.Previous studies have attributed this offset to the dust geometry, advocating for a bowlshaped rather than a spherical (doughnut-like) structure, caused by the anisotropy of the accretion disc emission (Kawaguchi & Mori 2010, 2011).This shape would also be seen if the dust formed part of a radiatively accelerated outflow from the disc, or a disc wind (Hönig & Kishimoto 2017).
Reverberation mapping has most commonly been carried out with photometric data, largely due to the relative ease in obtaining high cadence lightcurves for multiple objects simultaneously.However photometric reverberation mapping has limitations that can be overcome by using spectroscopy instead.The lightcurves are not true reflections of the continuum flux variability as, due to fairly wide bandwidth filters, they are often contaminated by broad line emission, which superimpose their own intrinsic variability originating from the BLR, and show reverberation, but from different radii than the hot dust emission (Peterson et al. 2004;Koshida et al. 2014).With time-resolved spectroscopic data these line profiles can be reliably separated from the continuum emission, providing not only an uncontaminated continuum lightcurve, but also allowing for the assembly of simultaneous emission-line lightcurves.This allows us to investigate the connection between variations of emission-line flux and profile shape simultaneously with dust continuum variations.In the scheme that the BLR is limited by the dusty torus, we could even gain insights into the shape of the BLR, current evidence would suggest that the BLR, like the torus could be bowl shaped (Goad et al. 2012).Most BLR reverberation mapping is carried out on the optical lines, however, and so is most probably sampling a smaller size scale than the most extended part of the BLR.In this work we extend emission-line reverberation studies to infrared lines, namely Pa and Pa.In addition to line measurements, spectroscopy also makes the simultaneous measurement of dust temperature possible, making astro-chemical analysis possible through sublimation arguments.This technique is now viable due to the advent of short-cross dispersed near-IR spectrographs mounted on medium sized telescopes, and has been utilised in Landt et al. (2019) which focused on NGC 5548.
In this work we present our findings from a 7-month long spectroscopic near-IR reverberation campaign on the nearby, well studied AGN Mrk 509 conducted between 2019 May to November.In this paper, we focus mainly on the near-IR continuum emission, but also investigate the variability of Pa and Pa.We ultimately are able to place lower limits on the hot dust and Paschen lags and thereby demonstrate the feasibility of measuring both of these values with a better sampled data set.
We adopt here cosmological parameters  0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7, which give a luminosity distance to Mrk 509 of 150 Mpc and an angular scale at the source of 682 pc per arcsec.

THE SCIENCE TARGET
This study focuses on the Type-1 AGN Mrk 509.We selected this target for several reasons: (i) it has a previously measured, short and variable dust response time of (126 ± 11) days (Koshida et al. 2014), therefore within our 6 month baseline of observations we should be able to recover a hot dust lag measurement; (ii) it has clearly discernible broad emission lines in the Paschen series allowing us to separate narrow-and broad-line components; (iii) multiple reverberation campaigns have determined its black hole mass.Peterson et al. (2004) reported the black hole mass to be (1.43±0.12)×10 8M ⊙ .In addition Mrk 509 has a spatially resolved torus radius of (296 ± 31) lt-days from near-IR interferometry (Gravity Collaboration et al. 2020).
Mrk 509 (J2000 sky coordinates R.A. 20 ℎ 44  09.8  and Decl.−10 • 43 ′ 24.7 ′′ ) can be observed from both the northern and southern hemispheres.Its relatively low redshift of  = 0.0344, enables us to sample the hot dust component and several near-IR coronal lines using a cross-dispersed near-IR spectrum.This AGN inhabits a compact host galaxy, perhaps of type S0 (Kriss et al. 2011) and is relatively bright in the near-IR (2MASS  = 11.6 mag,  = 10.0 mag; Skrutskie et al. (2006)) ensuring good quality near-IR spectra are possible using reasonable exposure times on a 4 m class telescope.

THE OBSERVATIONS
Our aim is to measure a reverberation signal in the near-IR hot dust emission and the Paschen emission lines of Mrk 509, and to measure the dust temperature.We therefore monitor Mrk 509 using near-IR spectroscopy and also require well sampled optical photometry to trace the accretion disc variations.

The near-IR spectroscopy
Using the Spex spectrograph at the 3 m NASA Infrared Telescope Facility (IRTF) on Manaukea, Hawaii, we observed the source Mrk 509 between 2019 May and November (semesters 2019A and 2019B) (Rayner et al. 2003).Having requested 36 nights of observing to obtain a cadence of ∼1 week, we were scheduled 24, 2-2.5 hour, observing windows between 2019 April and November (resulting in successful observations between May and November 2019).We obtained a total of 17 near-IR spectra as we lost seven of the scheduled observations due to weather and site inaccessibility.The journal of observations is listed in Table 1.
All observations were carried out using the (0.7-2.55) m short cross-dispersed mode (SXD) equipped with the 0.3"x15" slit, oriented at the parallactic angle.Using this narrow slit allowed us to minimise the contamination from the host galaxy, whilst giving a sufficiently high spectral resolution, of =2000 or full width at half maximum (FWHM) ∼ 150 km s −1 , to discern clearly any narrow and broad emission lines and study their profiles.We observed the source using a standard ABBA nodding pattern.As the source is not extended we nodded along the slit rather than on sky, therefore obtaining an image of the target in both the A and B frames from which the background subtraction could be carried out.The usual on-source exposure time was 32×120s.On seven of the seventeen observing runs we obtained fewer exposures than this due to either weather or technical issues.We also note that on 2019 September 10 we achieved 38 exposures due to spare time.This exposure time was chosen in order to guarantee high enough (S/N) to reliably measure not only the broad, but also the narrow emission line profiles, which can be key in performing photometric corrections to the spectra, as discussed in Section 3.3.
We observed a standard star to correct for telluric absorption and for flux calibration.Therefore after the science target we observed a nearby (in airmass and position) A0 V star, HD 203893 which has accurate optical magnitudes (B = 6.89 and V = 6.84).This spectral type has a well defined spectral shape into the infrared and using the quoted optical magnitudes along with the Spextool reduction pipeline this shape was used to perform the IR telluric correction.Where possible we took flats and arcs after the science target.Photometric corrections to the spectra in order to investigate intrinsic variability is carried out through a separate process detailed in Section 3.3.
The data were reduced using Spextool (version 4.1) a software package created for use by Spex users and based in Interactive Data Language (IDL) (Cushing et al. 2004).This package was used on both the science target and telluric standard.Spextool carries out all of the procedures required in the production of fully reduced spectra.This includes the preparation of calibration frames, the subtraction of flat fields, the fitting and subtraction of local background, interactive spectral extraction from science frames using an optimally weighted extraction scheme (Horne 1986), telluric correction of the science spectra, merging of orders, and combination of the different frames into a single reduced science spectrum for each observing window.The final spectrum was corrected for Galactic extinction using a value of  V = 0.079, which we derived from the Galactic hydrogen column densities published by Dickey & Lockman (1990).In Figure 1, we show the spectrum from 2019 May 9 as a representative example.

The LCOGT observations
Between 2019 April 2 and November 12 we monitored Mrk 509 in the optical with the 1 m and 2 m robotic telescopes of the Las Cumbres Observatory (LCOGT; Brown et al. 2013) almost daily in four bands (, ,  and   ).This baseline of observations begins ∼1 month before our near-IR spectra and ends 2 days before the end of our near-IR campaign.We use here mainly the lightcurves in the  and   filters, which have a central wavelength of 4770 Å and 8700 Å and a width of 1500 Å and 1040 Å respectively.The Sinistro cameras on the 1 m telescopes have a field-of-view (FOV) of 26.5 ′ × 26.5 ′ and a plate scale of 0.389 ′′ per pixel.The Spectral cameras on the 2 m telescopes have a FOV of 10.5 ′ × 10.5 ′ and a plate scale of 0.304 ′′ per pixel.
The frames were first processed by LCOGT's BANZAI pipeline (McCully et al. 2018) in the usual way (bias and dark subtraction, flat-fielding and image correction) and were subsequently analysed with our custom-made pipeline, which has been described in detail by Hernández Santisteban et al. (2020).In short, after performing aperture photometry with a diameter of 7 ′′ and subtracting a background model, we produced stable light curves by constructing a curve of growth using the standard stars on each individual frame and measuring the correction factors required to bring all different light curves to a common flux level.This approach mimicked extraction from an azimuthally averaged point-spread function (PSF), but without actually performing a PSF extraction.A colour correction and correction for atmospheric extinction were applied before the photometric calibration.Finally, we used comparison stars in each field to perform an image zero-point calibration at each epoch.In Figure 2, we show the  and in Figure 6 the   light-curve, with the latter overlapping in wavelength with the near-IR spectrum.

The ZTF data
The LCO -band lightcurve has been supplemented with data taken from the Zwicky Transient Facility (ZTF) (Bellm 2014).These data were corrected to match the flux scale of the LCO observations by interpolating the LCO flux lightcurve to the ZTF observation times, calculating a mean flux difference, and then subtracting this from the ZTF lightcurves.This correction value was ∼0.7 mJy.The resulting ZTF lightcurve is displayed in Figure 2 alongside the LCO observations.The ZTF data fills large gaps present in the LCO dataset, and therefore provides a helpful constraint on the behaviour of Mrk 509 during these time periods.

The GROND observations
We obtained simultaneous optical and near-IR photometry in seven bands with the Gamma-Ray Burst Optical and Near-Infrared Detector (GROND; Greiner et al. 2008) mounted on the MPG 2.2 m telescope, La Silla, Chile.The observations were taken on 2019 August   26, simultaneous with one of our near-IR spectra.The optical channels have a FOV of 5.4 ′ × 5.4 ′ and a plate scale of 0.158 ′′ per pixel, whereas both the FOV and plate scale of the near-IR channels is larger (10 ′ × 10 ′ and 0.6 ′′ per pixel, respectively).The filter passbands are relatively wide (∼ 1000 − 3000 Å).Images were reduced with the IRAF-based pipelines developed by R. Decarli and G. De Rosa (Morganson et al. 2012), modified and extended for our purpose.After bias and dark subtraction and flat-fielding of each individual image, a scaled median sky image constructed from separate sky exposures was subtracted from each frame.The dither step size was 18 ′′ for both the science and sky images, with five dither positions obtained for each.
In Table 2, we list the optical and near-IR nuclear fluxes.They were obtained with GALFIT (Peng et al. 2002), a software package that models the object's surface brightness profile with a point spread function (PSF), for the unresolved AGN, and a host galaxy component.We modelled the host galaxy with only a bulge due to the relatively low S/N of the images.In addition, GALFIT simultaneously models the PSF fluxes of the reference stars, which we used to flux-calibrate the optical and near-IR images based on the known stellar fluxes as listed in the SDSS (Abazajian et al. 2009) and 2MASS (Skrutskie et al. 2006).We note that we had to omit the -band images due to lack of SDSS reference stars in this area and that the host galaxy parameters in the  and  bands had to be constrained (based on the other images) rather than left to vary. Figure 3 shows clearly that the GROND photometry matches the spectral shape and flux level of our near-IR spectrum, which pro-vides a useful check of our telluric and correction, and indicates the lack of host galaxy contribution to our near-IR spectra.

Photometric corrections to the spectra
For the purposes of this work it is essential to distinguish between intrinsic AGN variability and extrinsic sources of variability such as changing night to night observing conditions.Therefore to obtain meaningful lightcurves from the spectra it was necessary to perform a photometric correction to the spectra.Failing to do this could result in a contamination of the final dust response lightcurves Bentz et al. (2016).
Several software packages employ different approaches to perform the photometric re-scaling of spectra, calibrating out timedependent instrumental effects, e.g prepspec (written by Keith Horne and discussed in Shen et al. (2016)) and GW92 (van Groningen & Wanders 1992).This work uses the open source python package mapspec (Fausnaugh 2017) which is an updated and improved version of the GW92 package.mapspec , like GW92, aligns an assumed non variable line profile of a given observed spectrum to a reference spectrum that has been constructed, by applying a smoothing kernel, a flux re-scaling factor, and a wavelength shift.
In order to correct for time-dependent instrumental effects, a spectral feature which is intrinsically non variable over the baseline of observations is required.by a pc-to-kpc scale region and therefore constant on timescales much larger than the baseline of observations typically used in reverberation studies such as this (Fausnaugh 2017).However, this emission line is not visible in the near-IR, therefore we have selected the [S III] 9531 narrow emission line, which is analogous to the [O III] 5007 region due to its extent.
A complication in the use of the [S III] 9531 narrow emission line is that it is blended with the variable broad Pa.It was therefore important to separate these two profiles to allow for accurate flux scaling.In order to do this a mean spectrum of our 4 photometric nights was constructed.Using this mean spectrum we isolated and fit the broad Pa to the broad Pa profile in velocity space and then converted back to wavelength space.This fit is displayed in Figure 4.
A linear fit to the Pa profile within the [S III] window of (9520 -9560) Å was performed to calculate the slope of an artificial continuum which could be used to isolate the intrinsically nonvariable [S III] narrow line from the variable broad Pa upon which it sits.Figure 5  Using mapspec, a reference spectrum was constructed from 4 photometric spectra using the artificial continua as displayed in Figure 4 as the continuum level, thereby isolating the [S III] profile.This reference spectrum was used alongside all 17 spectra with artificial continua in order to calculate re-scaling factors.mapspec allows for three different smoothing methods, a Delta function, a Gaussian and Gauss-Hermite smoothing kernel.We have taken the final scaling values as the Gauss-Hermite smoothing value as this is the most stable of the three methods (Fausnaugh 2017).We performed a test using two different continua, the broad line profile and the artificial continua shown in Figure 4, resulting in an average difference in scale factor of 3.5% between the two continuum options for the Gauss-Hermite method compared to an average difference in scale factor of 10.4% for the delta function, therefore confirming this as the most robust smoothing method.
To quantify the uncertainties from the scaling procedure we have calculated scaling factors using three alternative methods.Firstly by calculating the integrated flux in the [S III] line using the artificial continua and scaling these areas to the mean area of all four photometric nights.Secondly, using the interactive fitting package qdp, we modelled the [S III] region with a linear continuum plus broad and narrow Gaussian components.The linear continuum and broad line components were subtracted and the flux enclosed by the narrow line was calculated and scaled to the mean of the photometric nights.Finally we used prepspec to calculate a set of flux correction factors.In Figure 5, data points representing the photometric nights (enclosed in red) show very little spread in correction factor and all lie within a few percent of unity, as expected.Across all observations, the corrections factors for each method of calculation described above agree within ≲ 5%. Figure 5 shows that on average mapspec is closest to the mean correction factor for each observation, but the distribution around the mean for the other three correction factors appears random, indicating that one method alone does not significantly increase the spread displayed in the bottom panel.
The two central panels of Figure 6 show the mapspec corrected LCO   -band spectral fluxes compared with LCO   -band photometry.Whilst an exact agreement with the photometry would not be necessary to confirm the validity of the correction factors, it is clear that there are two measurements which are very discrepant with respect to the contemporaneous photometric measurements.These large discrepancies reduce the confidence in the scaling procedure and the assumptions upon which it is based, perhaps indicating that there is some extension in the [S III] region causing variability in this line between observations not due to atmospheric effects.We therefore scale to the LCO   -band photometry, as shown in the bottom panel of Figure 6.To do this, we bin the photometric data into 3-day windows around each spectral data point and scale to the weighted mean of this value, we then propagate the weighted standard error of this bin forward as the scaling uncertainty, and in any bin where this is less than the average spread of the data in any 3-day period we adopt this value as the uncertainty on the correction.
To test for potential spatial extension of the [S III] emitting region, we extracted spectra from the Integral Field Unit (IFU) image detailed in Fischer et al. (2015) (private communication from R. Wilman and T. Fischer), for each position angles (PAs) of the 17 spectra.By running mapspec over this dataset we found that the distribution of the extended emission is such that the addition of this relatively small contribution in each of the PAs we used does not result in variations above a few percent.

Variability in the broad Paschen emission lines
As Mrk 509 is a Type 1 Seyfert it exhibits both broad and narrow emission lines in its optical and near-IR spectra.The Paschen series are a group of hydrogen recombination transitions that fall within the spectral range of the IRTF Spex instrument, and therefore are detected in our spectra.Figure 1 shows these broad emission lines, the most prominent of which are Pa, Pa and Pa.It is clear from their profile shapes that these emission lines are emitted from the BLR, which is observed to vary due to its illumination by the accretion disc (Peterson et al. 2004).Therefore it should be possible to measure variability in the Paschen lines from our spectra, and to search for a reverberation signal from our 'driving' optical LCO lightcurve and thereby place constraints on the location of the Paschen BLR.
To do this we have selected the Pa and Pa lines, opting to reject Pa due to the effect of telluric absorption on its blueward wing, which can be seen in Figure 1, but is far more prominent on nights of poorer quality.Pa is surrounded by a relatively clean continuum, however there is a slight complication as on top of its broad profile sits [S III], procedures for the removal of Pa from the The root mean square (rms) fractional variation ( var ) is a measure of fractional "excess variance", and therefore is a measure of intrinsic variability above that due to flux measurement errors. var is calculated as, following Rodríguez-Pascual et al. (1997), where ⟨⟩ is the sample mean of the  = 17 flux data   ; and  is the square root of the sample variance, calculated as, (3) Finally,  2 is the mean square uncertainty of the fluxes, calculated as, where   is the uncertainty on each flux measurement   .
Figure 7 shows the  var spectra (calculated as detailed in Equation 1) for the Pa region.Firstly we plot the  and  values, defined in Equations 4 and 3 respectively, as a function of wavelength.Secondly the  var values are plotted along with their 1 uncertainties.We compute and plot an upper limit for any spectral element for which the  value exceeds the  value.
We note that applying a percentage error to each flux element of the spectra leads to an apparent dip in the broad line region of the  var spectrum.This artefact is resultant from the increase in the absolute error value for spectral elements of higher flux values.This produces a bias in the weighting step as the  values increase (black line Figure 7) whereas the  values do not (blue line Figure 7).Therefore spectral elements with higher flux will be weighted more heavily when calculating  var thus resulting in a dip in the broad line profile  var values.
However, in the case of Pa the [S III] narrow line profile shows an  var amplitude of ∼ 6%.This indicates that when scaled to the photometry, the [S III] profile is variable.As the  var has been calculated for each spectral element, the [S III] profile calculation does include the continuum variability which would increase  var , but also includes the systematically increased errors which would act to reduce  var .Despite these competing effects, we still detect variability in the [S III] line profile above the level we detect for the continuum.
Whilst the narrow-line region has been found to vary over large timescales, it is unlikely that we are sensitive to this over the 183 days of our observations.The narrow-line region has been found to be extended over large physical size scales, meaning that we could see variation in the [S III] region due to our chosen slit width and observational effects such as seeing (Peterson et al. 2013).For reasons discussed in Section 3.3 we believe scaling to the photometry to be more reliable than assuming a non-variable [S III] component.

The spectral continuum components
We aim to compare response-and luminosity-weighted dust radii.For the dust time delay we assembled the lightcurves of both the irradiating flux and hot dust.For the calculation of the luminosityweighted dust radius we measured the dust temperature and estimated the total accretion disc luminosity that heats the dust.The large wavelength coverage of our cross-dispersed near-IR spectra gives in Mrk 509 roughly half of the hot dust SED, and the improved wavelength coverage of SpeX in the blue samples a considerable part of the accretion disc spectrum, which is expected to dominate the total continuum flux up to a rest-frame wavelength of ∼ 1 m (Landt et al. 2011a,b).We note that, since we used a relatively small spectral aperture of 0.3" which corresponds to a physical size of ∼205 pc at the source, the contribution from the host galaxy to the total observed continuum flux is expected to be small in this luminous AGN.
Following Landt et al. (2019), we decomposed the spectral continuum into two components.We first approximated the rest-frame wavelength range of ∼ < 1 m with an accretion disc spectrum, which we subsequently subtracted from the total spectrum.For its calculation we adopted a black hole mass of  BH = 1.09 × 10 8 solar masses, which is the value derived by optical reverberation campaigns (Peterson et al. 2004;Bentz & Katz 2015) using a geometrical scaling factor of  = 4.3 to convert the measured virial product (Grier et al. 2013).Furthermore, we assumed that the disc extends to  out = 10 4  g , where  g =   BH / 2 is the gravitational radius, with  the gravitational constant and  the speed of light, the value of  was resultant from each fit and an efficiency of  = 0.057 was assumed.We then fitted the resultant hot dust spectrum at wave-lengths > 1 m with a blackbody, representing emission by large dust grains, and with two blackbodies modified by a power-law of the form   () ∝   , approximating with  = −1 and  = −2 the emissivity of sub-micron silicate and carbon dust grains, respectively (Landt et al. 2019, see their Figure 8).Despite our fitting of only two components, there is evidence to suggest that at 1m diffuse BLR continuum emission could account for 25% of the total continuum emission with this increasing to 35% at 8000Å (Korista & Goad 2019).This could have an effect on our spectral decomposition, and on our spectral lightcurves, particularly those extracted at wavelengths <1m, see Section 4.1.Table 3 lists the relevant physical parameters extracted from the spectral decomposition.We obtain average temperatures of ⟨⟩ = 1365 ± 4 K, 1166±3 K and 1018±3 K for emissivity laws with  = 0, −1 and −2, respectively.
As in Landt et al. (2019), we calculated luminosity-weighted dust radii,  d,lum , from the best-fit dust temperatures assuming radiative equilibrium between the luminosity of the irradiating source and the dust: where  is the Stefan-Boltzmann constant and ⟨ em ⟩ is the Planck-averaged value of   ().We approximated  uv with the accretion disc luminosity and have used for the Planckaveraged emission efficiencies in the case of  = −1 a value of ⟨ em ⟩ = 0.0210 appropriate for silicates of  = 1259 K and  = 0.1 m (Laor & Draine 1993) and in the case of  = −2 a value of ⟨ em ⟩ = 0.0875 appropriate for graphite of  = 1000 K and  = 0.1 m (Draine 2016).The average luminosity-weighted dust radii are ⟨ d,lum ⟩ = 186 ± 4 light-days, 1763±36 light-days and 1132±23 light-days in the case of a blackbody, and small-grain silicate and carbon dust, respectively.

The observed spectral continuum light-curves
Lightcurves for four spectral continuum windows are constructed from fluxes extracted from near-IR spectra and displayed in Figure 8.In addition, we have the -band photometric light-curve for the accretion disc, shown in Figure 2. Ideally, we would like to measure the accretion disc flux at the shortest wavelengths, where it is least contaminated by the hot dust, however, the S/N ratio of the near-IR spectra decreases significantly towards short wavelengths.Therefore, we measured the accretion disc flux in the 60 Å wide rest-frame wavelength region of  = 9730 − 9790 Å, which lies between the two broad hydrogen emission lines Pa and Pa and is known to be line-free, and also in the wider (300 Å) rest-frame wavelength region of  = 8700−9000 Å, which is not contaminated by significant emission-line flux in Mrk 509 (see Figure 1).
We measured the flux in two line-free, 500 Å wide rest-frame wavelength regions, namely,  = 2.25 − 2.30 m, which is towards the red end of the  band and thus close to the peak of the blackbody spectrum, and  = 1.50 − 1.55 m, which is close to the middle of the  band.These bands were selected due to the absence of any emission line contamination and the lack of any telluric effects even in the worst quality spectra.These two wavebands have a significantly lower contribution from the accretion disc than the two shorter wavelengths selected.In Figure 8 the 30% rise in -band flux between MJD ∼57830 and ∼58800 has no clear counterpart in the H and K bands as opposed to the shorter-wavelength bands  The columns are: (1) Universal Time (UT) date of observation; (2) total accretion disc luminosity; for a blackbody emissivity (3) dust temperature; (4) total dust luminosity and (5) dust radius; for an emissivity law appropriate for silicate dust with small grain sizes of  ∼ < 0.1 m (6) dust temperature; (7) total dust luminosity and (8) dust radius; for an emissivity law appropriate for carbon dust with small grain sizes of  ∼ < 0.1 m (9) dust temperature; (10) total dust luminosity and (11) dust radius.
which possess this distinctive rise also seen in our -band lightcurve, thus demonstrating the reduction in accretion disc contribution well.
The reduction in accretion disc contribution at longer wavelengths is well documented in Mrk 509 and similar sources (Hernán-Caballero et al. 2015, 2017;Kuraszkiewicz et al. 2003).Hernán-Caballero et al. (2017) show that at H the contributions from the accretion disc and hot dust are roughly equal, whereas at K the hot dust emission significantly dominates over that of the accretion disc.
It is worth noting that the shorter wavelength variability that we sample could contain a significant contribution from diffuse BLR continuum emission, or a wind on the inner edge of the BLR, rather than simply being from the disc (Netzer 2021;Kara et al. 2021).
It is clear from Figure 8 that the level of detected variability decreases with increasing wavelength, with the most significant variability being detected in the -band.
From the (126±11) lt-day lag measured by Koshida et al. (2014) and the interferometrically determined radius of (296±31) lt-days (GRAVITY Collaboration: et al. 2020), along with our luminosity weighted dust radius of (186±4) lt-days, it appears that we should expect to measure a lag of (100-200) days.A time lag of this size would necessarily require a smoothing timescale of a similar order meaning that the hot dust emission originating from the torus would not display a high level of variability on short timescales.This could explain why, with our limited baseline and number of observations, 17 data points across 183 days, we detect little variability in the hot dust emission.We would only expect the already low level of variability seen in the H and K-band lightcurves to decrease if the accretion disc component were removed, and therefore we would almost certainly not be sensitive to variability in the hot dust emission alone.

Emission line light-curves
To construct lightcurves, we extract the Pa and Pa line fluxes.For Pa we first fit a linear continuum to redward and blueward clean continuum regions, shown as the red dot-dashed line in Figure 4. We subtract this linear continuum, and then integrate under the data using Simpson's method.We must subtract the [S III] component from each observation in order to obtain the flux enclosed by Pa.As detailed in Section 3.3 we opt to scale to the photometry due to major discrepancies in the scaling procedures relying on the [S III] being constant.Therefore the subtracted [S III] values vary by ∼ 6% in line with our calibration procedure.The flux extraction process was slightly more straightforward for the Pa profile: we fit a linear continuum to relatively clean blueward and redward continuum regions, again displayed by the dashed line in Figure 4.The selected redward continuum window is necessarily short in this case as this region is affected by telluric absorption.We then subtract this continuum and simply integrate under the line, using Simpson's method, to give the fluxes displayed in Figure 9.The Paschen line flux uncertainties shown in Figure 9 were calculated using the Monte-Carlo method, with 1000 iterations, selecting the 1 bounds of the resultant Gaussian distribution and then numerically propagating this with a 2% systematic from the continuum selection.
Figure 9 shows the lightcurves for Pa and Pa, it is clear that there is not a high level of variability.This is confirmed by the calculated  var values listed in Table 4.

Continuum and hot dust
JAVELIN is an open source reverberation mapping package (Zu Fluxes have been calibrated onto an absolute scale using mapspec scale factors and then a secondary photometric correction, therefore any variability is considered to be intrinsic.Utilising the LCO and ZTF -band data as a 'driving' lightcurve, and the four spectroscopic near-IR lightcurves shown in Figure 8 as 'response' lightcurves we run JAVELIN with 10 5 iterations and a (0-100) day limit on the time lag search range in the form of strong preferential priors.This lag range was selected in order to prevent a significant number of spectral data points being shifted into a seasonal gap in the 'driving' lightcurve.The results are displayed in Figure 10.
Figure 10 shows that the (0.870-0.900) m and (0.973-0.979) m bands clearly reverberate on timescales shorter than ∼25 ltdays, whereas the H and K-band lightcurves tend towards longer lags.This is to be expected given that the two shorter-wavelength bands are dominated by the accretion disc, which would reverberate on much shorter time scales than any hot dust emission which would dominate to a larger degree at H and even more so at K.
Due to the limited sampling of our data and ∼150 day seasonal gaps in the 'driving' -band lightcurve, there are strong aliasing problems, and we cannot therefore report any confident lag.However we can state that the H and K-band lightcurves do not reverberate at <40 days.This is consistent with previous findings.Koshida et al. (2014) measure a reverberated dust lag in the K band of (126.8±11)lt-days, and Gravity Collaboration et al. ( 2020) measure the spatially resolved torus at a radius of (296 ± 31) lt-days using near-IR interferometry.
We are not only limited by a lack of significant variability features over the majority of our campaign, but also the data gaps present in both the spectral dust, and photometric continuum lightcurves.The large seasonal data gap shown clearly in Figure 2 between MJD ∼58450 and ∼58550 limits our ability to measure lags from the hot dust.The -band lag measured by Koshida et al. (2014) of (126.8±11)lt-days would be extremely difficult to measure with our data-set as a lag of this length would place much of the spectral data into a seasonal data gap, and the length of the lag would smooth out what little variability we do sample in the hot dust response.In fact any lag of more than ∼65 days would place a significant proportion of our already limited data sample into this seasonal gap, which greatly limits the ability to measure a credible hot dust response with our current sampling.
In Appendix:A we present a wider range of positive and negative lags to demonstrate the prevalence of alias lags in all of the spectral lightcurves.

The broad emission line region
Literature values for the H lag in Mrk 509 place the Balmer BLR at 82.3 +6.3  −5.6 days and 79 +6 −7 days for  cent and  peak respectively, using Cross-Correlation Function (CCF) analysis (Peterson et al. 2004).Unless Mrk 509 has changed state significantly since these measurements were taken (1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996), we would expect to see the Paschen line lags to at least be this large.As with measuring reverberation signals in the hot dust of Mrk 509, this presents a problem due to the sparse sampling, short baseline of observations and data gaps present in our data.
As in the case of the hot dust lightcurves, we used JAVELIN to search for Paschen line reverberation signals, the results of which are displayed in Figure 11.
JAVELIN reports a lag of 62 +7 −6 days in Pa.This result is consistent (within 2) with previous measurements of the H lag (Peterson et al. 2004).Our reverberation result of 62 +7 −6 days, does show some smoothing, with a value of 21 +29 −19 days, however this is not as large as we would expect for a lag of this size.In addition, Figure 11 shows that with a lag of 62 days, 3 out of 17 data points are shifted into a data gap and do not align with any optical data points.Upon visual inspection, we see that the remaining 10 data points do not show a good agreement with the pattern of variability displayed by the optical lightcurve.Therefore we do not report this result with a high level of confidence.However, as with the  and -band spectral lightcurves, we are able to place a lower limit on the Pa lag of ∼40 days as we see no counterpart in the lightcurve of the distinctive rise in flux in the last 80 days of the campaign, marked in red in Figure 2.
We can place this same lower limit on Pa which shows the same distinct absence of the major rise in flux present in our 'driving' lightcurve.We measure the Pa at 91 +7 −37 days, shown in Figure 12.This lag places 6 of our 17 data points into a data gap and therefore is not reported with a high level of confidence despite being consistent with previous measurements of the H lag (Peterson et al. 2004).
JAVELIN results for Pa and Pa are displayed in separate plots as the analysis was carried out separately.Running JAVELIN separately resulted in a fit with higher level of smoothing in both cases, and therefore in some senses a more physical result.
In Appendix:A we present a wider range of positive and negative lags to demonstrate the prevalence of alias lags in both Pa and Pa.

SUMMARY AND CONCLUSIONS
We monitored the AGN Mrk 509 to document its infrared (hot dust) and optical (accretion disc) variations.Our IRTF spectra span 0.7 to 2.4 microns at 17 epochs over 200 days from 2019 May to Nov. The -band monitoring spans 5 seasons with sub-day cadence during the IR campaign.
Using spectral decomposition we measure luminosity weighted dust radii for different astro-chemical regimes to be ⟨ d,lum ⟩ = 186 ± 4 light-days, 1763±36 light-days and 1132±23 light-days in the case of a blackbody, small-grain silicate and carbon dust, respectively.Comparison with the previously measured photometric dust response time of (126±11) days (Koshida et al. 2014) and an interferometrically determined radius of (296 ± 31) lt-days (Gravity Collaboration et al. 2020) suggests that the dust emits like a black body and therefore is most likely to be composed of large carbonaceous grains.This is consistent with results for the Seyfert NGC 5548 (Landt et al. 2019).
Our -band lightcurve shows a distinctive rise by ∼30% in the last 80 days of the 2019 season, with no clear counterpart in the  and  lightcurves, which have intrinsic variations less than ∼4%.The little variability we do detect in H and K does not reverberate on timescales less than 40 days and therefore we find a lower limit on the dust response in Mrk 509 which is consistent with previously measured reverberation measurements and also interferometric measurements of the hot dust radius (Koshida et al. 2014;Gravity Collaboration et al. 2020).
We are also limited in our ability to recover a dust lag due to seasonal gaps in both the driving and response lightcurves.This combined with the poor cadence and low levels of variability in the dust response means that application of a large smoothing width would remove any significant variability, which is already low.
We measure lags of 62 +7 −6 days and 91 +7 −37 days for Pa and Pa respectively.The Pa lag is considered a more reliable measurement than that determined for Pa due to a more concentrated peak, fewer data points placed in seasonal gaps of the 'driving' lightcurve and a more reasonable smoothing result.The Pa result is consistent with literature values within 1 and the Pa result consistent within 2 (Peterson et al. 2004).However due to the seasonal gaps in our 'driving' lightcurve and poor sampling of the emission line lightcurves, we do not report either of these results with a high level of confidence.As with the  and -bands, we can place a lower limit on the lag of ∼40 days as we do not detect any lags below this value and do not recover the distinctive variability feature seen in our 'driving' lightcurve during the last ∼80 days of our spectral campaign.
A longer baseline of observation, and a higher cadence is needed to recover a reverberated signal from the hot dust and place constraints on the physical size scales of this emission region.

Figure 1 .
Figure 1.IRTF SpeX near-IR spectrum from 2019 May 9 shown as observed flux versus rest-frame wavelength.Emission lines listed in Table4ofLandt et al. (2008) are marked by dotted lines and labeled; black: permitted transitions, green: permitted Fe ii multiplets (not labeled), red: forbidden transitions and cyan: forbidden transitions of iron (those of [Fe ii] not labeled).

Figure 2 .
Figure 2. LCO -band observations of Mrk 509 with complementary corrected ZTF  band data.The most significant variability feature measured contemporaneously with our near-IR spectra is shaded in red.Archival LCOGT data (MJD<58200) is shown to extend the time baseline.Spectral epochs are marked by red arrows.

Figure 3 .
Figure 3. IRTF SpeX near-IR spectrum from 26 th August 2019, highlighted in Figure 6 (black), overlaid with the optical and near-IR nuclear fluxes from the GROND photometry obtained on the same day (red filled circles).
top panel shows the Pa and [S III] emission region from all 17 Mrk 509 spectra with the real continuum and the broad Pa profile, with the bottom panel displaying the [S III] scale factors calculated in four different ways.

Figure 4 .
Figure 4. Top Panel: Pa emission line from Mrk 509 spectrum taken on 26 August 2019.Red line represents a linearly fitted continuum, subtracted in order to calculate enclosed flux.The continuum has been fit using a blueward clean continuum window of (12500-12640) Å and the redward clean continuum window of (13000-13050) Å, these regions are highlighted in blue.The redward window is necessarily short as beyond this there is significant telluric absorption.Bottom Panel: Mean spectrum of Mrk 509 from all 4 nights with photometric conditions.Pa profile (black) scaled to match Pa (red), in order to isolate the constant [S III] 9531 component from any variable broad line emission.The slope of the Pa profile in the [S III] window is shown as the artificial continuum in grey.Both panels centre on the rest frame wavelengths of 12821.6Å and 9548.6 Å for Pa and Pa respectively and have a matched velocity range.

Figure 5 .
Figure 5. Top Panel: All 17 Mrk 509 spectra, showing the region around Pa and [S III].The four photometric spectra are coloured red, the others black.Second Panel: [S III] flux correction factors calculated using four different methods for each night of observation, photometric nights are enclosed by red boxes.The four methods used are as follows; mapspec using spectra with an appended artificial continuum (see Figure 4), interactive fitting routine qdp with a linear, narrow Gaussian and broad Gaussian component, numerically integrating under the [S III] line profile after the removal of the artificial continuum and prepspec.Third Panel: Percentage of the mean correction factor for each value obtained via the four methods described above, coloured similarly to the panel above.Bottom Panel: Percentage range of the four methods around their mean.

Figure 6 .
Figure 6.Top Panel: LCO   band photometry, uncorrected for host galaxy contribution in black, with all 17 spectral LCO   band (8180-9220 Å observed frame) uncorrected fluxes plotted in red.Second Panel: LCO   band photometry, uncorrected for host galaxy contribution in black, with all 17 spectral LCO   band ((8180-9220) Å observed frame) fluxes, corrected with mapspec scaling factors are plotted in red.Third Panel: LCO   band photometry, corrected for a constant host galaxy contribution in black, with all 17 spectral LCO   band (8180-9220 Å observed frame) fluxes, corrected with the mapspec scaling factors plotted in red.Spectral fluxes measured on MJD 58745.30and 58801.20 are encircled in red as they do not match well with the LCO   band photometry.The night of 26 th August 2019 is highlighted with a dashed red line.Bottom Panel:The mapspec corrected lightcurve has been further corrected to the LCO   band photometry by scaling to the weighted mean of all photometric data points within three days of each spectral data point.The errors displayed are defined by the spread of the data to which the spectral points have been corrected, or the mean spread of the photometry within any given 3 day period of the spread if the correction was less than this value.These uncertainties have been propagated to our corrected lightcurves along with the mapspec uncertainties displayed in the panel above.

Figure 7 .
Figure 7. ,  (top panels) and  var spectrum (bottom panels) for Pa (top plot) and Pa (bottom plot) calculated through the process described in Section 3.4. var values are shown in red with their 1 uncertainties calculated via Monte Carlo methods, represented by the red shaded region, upper limits on  var are shown as red triangles. values are displayed as the blue lines along with their 1 uncertainties calculated via Monte Carlo methods. values shown as the black solid line.

Figure 9 .
Figure 9. Pa and Pa emission line lightcurve, constructed using the procedure described in Section 4.2.

Figure 10 .
Figure 10.JAVELIN results for the four spectroscopic bands, (0.870-0.900) m, (0.973-0.979) m, (1.50-1.55)m, (2.25-2.30)m, displayed from top to bottom.Each run consists of 10 5 iterations, with the determined lags displayed on the right hand of each panel.The 17 red data points represent spectral data points, and the continuous grey curve is the smoothed, shifted and scaled JAVELIN model of the 'driving' -band lightcurve displayed in Figure 2. Lag limits placed at (0-100) days.

Figure 11 .Figure 12 .
Figure 11.JAVELIN results for the Pa broad emission line.Each run consists of 10 5 iterations, with the determined lags displayed on the right hand of each panel.The 17 red data points represent enclosed emission line fluxes, and the continuous grey curve is the smoothed, shifted and scaled JAVELIN model of the 'driving' -band lightcurve displayed in Figure 2. Lag limits placed at (0-100) days.

Table 3 .
Physical parameters for the calculation of luminosity-weighted dust radii