© The Authors 2024

1. Introduction

Quasars are unique laboratories for the extreme physics governing active galactic nucleus (AGN) accretion and feedback and they are important drivers for galaxy evolution and enrichment (e.g., Sijacki et al. 2007; McNamara & Nulsen 2012; Moustakas et al. 2019). However, since the central source in a quasar can exhibit a visual luminosity comparable to the entire host galaxy that it resides within, the point spread function (PSF) of a quasar’s central source (when seen with a telescope with finite mirror size) often dominates the light at inner approximate kiloparsec (∼kpc) scales. Many features of quasars’ circumnuclear regions (e.g., Ramos Almeida & Ricci 2017), such as inflows, dusty tori, winds, and jets, have visual and infrared (IR) components that appear overwhelmed as a result (e.g., Ford et al. 1994, 2014). Moreover, the dynamics and morphology of the host galaxy close-in to the nuclear region are likewise “swamped out” in visible and near-infrared (NIR) observations. While advances in radio interferometry have allowed us to glimpse the event horizons of two supermassive black holes at a ∼0.1 mas scale and their surrounding environments (i.e., M87: Event Horizon Telescope Collaboration 2019, Sagittarius A*: Event Horizon Telescope Collaboration 2022; Lu et al. 2023), many processes critical to feeding and feedback in extremely luminous quasars will only be understood once we obtain high-contrast, high-resolution imaging in the visible and NIR (e.g., Martel et al. 2003; Gratadour et al. 2015; Brandl et al. 2008; Moustakas et al. 2019; Rouan et al. 2019; Grosset et al. 2021; Ding et al. 2023).

In studies of quasar hosts, IR and radio interferometry are capable instruments used to image the inner few parsecs of the circumnuclear region. This allows us to take look at the broadline region, inner radius of the torus, and jets on this scale (e.g., Kishimoto et al. 2011; Lister et al. 2013; GRAVITY Collaboration et al. 2018). However, the circumnuclear disk, infalling material, and jet activity in the narrowline region are best observed at an intermediate scale in the regions further out. To study these regions, observations of naturally dust-obscured quasars (i.e., “natural coronagraphs,” e.g., Jaffe et al. 1993; van der Marel & van den Bosch 1998) only preferentially sample these structures in disk-like hosts or during epochs of peak dust production early on in mergers, giving an incomplete picture of quasar evolution (Urrutia et al. 2008; Schawinski et al. 2012; Del Moro et al. 2017). However, the Space Telescope Imaging Spectrograph (STIS) coronagraph on board the Hubble Space Telescope (HST) can fill the gap at > 0.5 kpc to ∼kpc scales and larger to study the morphology of host galaxies. The visual imaging of the sub-kpc structures of the dust, jet, and host galaxies of quasars with the STIS coronagraph provides a unique opportunity to complement ground- and space-based IR high-contrast imaging with Keck and JWST.

The prototypical quasar 3C 273 was first identified based on its redshift of 0.158 by Schmidt (1963). With the central quasar dominating the signals across visible to radio wavelengths, thus overwhelming the host galaxy, the study of the latter makes it necessary to first properly remove the quasar light. In visible wavelengths, Martel et al. (2003) placed the 3C 273 behind a coronagraph using HST/ACS, removed the quasar light using reference star images, and revealed the host galaxy in visible light exterior to an angular radius of $\sim 1\stackrel{″}{.}5$. With non-coronagraphic imaging after empirical PSF removal, these non-asymmetric signals persist after host galaxy modeling (Zhang et al. 2019, Fig. 10 therein). Using (sub)-millimeter observations, Komugi et al. (2022) subtracted point source models using ALMA to reveal the surroundings where the millimeter continuum emission is colocated with the extended emission line region in [O III] observed with VLT/MUSE in Husemann et al. (2019). These efforts have unveiled complex structure for the 3C 273 host from multiple aspects, calling for dedicated imaging that would better reveal and characterize the 3C 273 host with available state-of-the-art instruments. In this work, we report our coronagraphic high-contrast imaging observations using HST/STIS, which reached an inner working angle (IWA) of $\sim 0\stackrel{″}{.}2$ in visible light to reveal the 3C 273 host galaxy.

2. Observations and data reduction

The HST/STIS coronagraph offers broadband imaging in visible through NIR light (0.2 μm–1.15 μm; Medallon & Welty 2023). Its narrowest occulter, BAR5, can image the surroundings of bright sources exterior to $\sim 0\stackrel{″}{.}2$ from the center (Schneider et al. 2017; Debes et al. 2019). To reveal the surroundings of bright central sources using STIS, a careful selection of reference stars is needed to avoid color mismatch given it is in the broadband (Debes et al. 2019); otherwise, a non-matching PSF can induce spurious signals (Ren et al. 2017, Fig. 8 therein). A reference star should be of similar color, magnitude and in close proximity to a science target to maximize the success in coronagraphic imaging (Debes et al. 2019). To explore the inner regions of the 3C 273 host galaxy down to $\sim 0\stackrel{″}{.}2$, or ∼0.5 kpc¹, we observed it with two reference stars using HST/STIS in HST GO-16715 (PI: B. Ren). The two reference stars serve as the empirical PSF templates to remove the 3C 273 quasar light.

2.1. PSF reference star selection

On the lower Earth orbit, HST is affected by instrumental effects such as breathing caused by changes in the Solar angle, radiation, and Earth’s shadow. These effects on STIS observations could be empirically captured and removed when well-chosen PSF stars are close in angle to the science target. To effectively remove the PSF from the central source, Debes et al. (2019) recommended close match in magnitude and color in the B and V band. With the high-sensitivity space-based measurements in visible light from the Gaia mission (DR3: Gaia Collaboration 2023), Walker et al. (2021) showed that magnitude and color match in Gaia filters may also provide reasonably good PSFs for STIS observations.

Using Gaia DR3, we selected two reference stars, TYC 287-284-1 (hereafter “PSF1”) and TYC 292-743-1 (hereafter “PSF2”), to serve as the PSF templates for 3C 273 (hereafter “science target”). More details are given in Appendix A. For 3C 273, its Gaia DR3 magnitude G = 12.84, color Bp − Rp = 0.494, and G − Rp = 0.348. Within $4\stackrel{°}{.}7$ from 3C 273, the chosen PSF1 has Gaia DR3 magnitude G = 12.77, color Bp − Rp = 0.515, and G − Rp = 0.327. Within $2\stackrel{°}{.}7$ from 3C 273, the chosen PSF2 has magnitude G = 11.488, color Bp − Rp = 0.536 and G − Rp = 0.340. By choosing two PSF reference stars with such faintness that do not have existing IR excess measurements, with IR excess being indicative of circumstellar disks (e.g., Cotten & Song 2016), we can reduce the probability that a reference star might host a circumstellar disk that negatively impacts the PSF removal for 3C 273. In addition, this strategy reduces the possibility of a negative impact of background objects, which are beyond the detection limits of exiting instruments, on the PSF removal of 3C 273.

2.2. Observations

Coronagraphic imaging with HST/STIS relies on the blockage of light in central regions using its physical occulters. In STIS, two nearly perpendicular occulting locations, BAR5 and WedgeA0.6, offer IWAs of $\sim 0\stackrel{″}{.}2$ and $\sim 0\stackrel{″}{.}3$, respectively (Medallon & Welty 2023). To enable a full angular coverage of extended structures, however, the two locations are near to the edges of the field of view of STIS, making it unrealistic to roll the telescope ∼90° to achieve a nearly 360° coverage given the scheduling limits of HST using four consecutive orbits (e.g., Fig. 3 of Debes et al. 2019)², see Appendix A. Therefore, we scheduled two sets of observations in GO-16715, with each set composed of four contiguous “back-to-back” HST orbits (see Table A.1 for the observation log). By observing 3C 273 at two different epochs spanning about two months in HST Cycle 29, the relative roll is $84\stackrel{°}{.}032$ between the central visits of the two epochs to approach a full angular coverage.

In GO-16715, we observed only one object in an orbital visit, with a sequence of “target-target-PSF-target” in a four-orbit observation set. This ensures that the telescope thermal distribution is stabilized when a PSF star is observed. The observations were in CCD Gain = 4 mode to permit high dynamic range imaging (e.g., Debes et al. 2019). For the three orbits on the target, we rolled the telescope by either ±15° (UT 2022-01-08) or ±5° (UT 2022-03-26)³ to approach a ∼360° angular coverage for it (see Fig. B.1 for the coverage map of 3C 273).

Within one science target orbit, we observed 3C 273 using both the BAR5 and the WedgeA0.6 occulting location with three readouts each. With three readouts, each having a 315 s exposure at an occulting location, the STIS flat-fielded files can identify cosmic rays for random noise removal. There are a total of 36 readouts from six orbits. For the reference stars, on the one hand, the Gaia DR3 G magnitude of PSF1 is similar as that of the science target, with PSF1 being 0.07 mag brighter. The observation strategy of PSF1 is identical to the target: to reach similar detector counts, each readout of PSF1 is 294 s. There are six readouts in total for PSF1. On the other hand, the Gaia DR3 G magnitude of PSF2 is 1.352 mag brighter than that of the science target. To reach similar detector counts as 3C 273, each readout of PSF2 is 125 s. Given the available time in one HST orbit, the relatively shorter readout time permit the dithering of the telescope: the on-sky step is 0.25 STIS pixel to reduce the impact from the non-repeatability of HST pointing on STIS results (e.g., Ren et al. 2017; Debes et al. 2019). At each occulting location, we dithered the telescope twice with each dithering location having two 125 s readouts, totaling six readouts per occulting location. There are a total of 12 readouts for PSF2.

2.3. Data reduction

2.3.1. Pre-processing

In long readouts (315 s for the target), the observation data quality with STIS may be compromised due to different noise sources (e.g., cosmic ray, shot noise, charge transfer inefficiency: CTI). To address this in data reduction, we first used the stis_cti package⁴ that applies the Anderson & Bedin (2010) correction to remove CTI effects for STIS CCD⁵. We then used the data quality map in the CTI-corrected flat-fielded files and followed Ren et al. (2017) to perform a median bad pixel replacement for the data that have been marked as a bad pixel (e.g., pixels with dark rate more than 5σ the median dark level, bad pixel in reference file, and pixels identified in cosmic ray rejection) around its 3 × 3-pixel neighbors. To correct for the geometric distortion in 2-dimensional STIS CCD images (Hodge et al. 1998), we then rectified the CTI-corrected images using the x2d function from stistools⁶.

In STIS coronagraphic imaging, the central star is blocked by the STIS occulters, therefore, we used the two diffraction spikes to align the images (e.g., Schneider et al. 2014; Ren et al. 2017). We used the “X marks the spot” method (e.g., Schneider et al. 2014): first, we fit Gaussian profiles to each column or row of the observation to identify the peak of the diffraction spikes and then fit the lines to obtain the intersection point as the location of the star. In this way, we can obtain both the locations and their associated uncertainties from the observed data.

2.3.2. Post-processing: PSF removal

For each target image, we minimized the standard deviation for the regions containing diffraction spikes after PSF subtraction to obtain a residual image. Specifically, we used the algorithmic mask from Debes et al. (2017) to identify the regions that are blocked by the STIS occulters in the observations. At each occulting location, we scaled the median of all the PSF readouts, and subtracted it from the target readouts to inspect the residual images. We derotated each residual map to north-up and east-left using its corresponding ORIENTAT header and used the element-wise median of all 36 residual maps (18 from each occulting location) as the final image for 3C 273 host galaxy. We multiplied the final image in units of counts s⁻¹ pixel⁻¹ by PHOTFLAM = 4.22 × 10⁻¹⁹, which is the inverse sensitivity parameter recorded in the FITS file headers, to obtain the calibrated final image in units of erg s⁻¹ cm⁻² Å⁻¹, see Fig. 1. We also calculated the standard deviation for the regions that do not host signals in the final image to obtain the signal-to-noise ratio (S/N) map (see Fig. B.3).

Fig. 1. 3C 273 host galaxy and surroundings in visible light seen with the HST/STIS coronagraph. The surface brightness is in log scale, and the central dark regions are not accessible due to coronagraph blockage. The data used to create this figure are available at the CDS.

To inspect the PSF removal quality dependence on different PSF references, we experimented using a single PSF star for all the images and we did observe a decrease in the S/N after combining the two PSF stars. In fact, in the observations for PSF1, we identified at least two faint background sources within the STIS field of view (e.g., Fig. A.2). Nevertheless, this was anticipated in our observation planning, and we minimized such risks by using multiple PSF stars in Table A.1.

3. Results

3.1. 3C 273 host

3.1.1. Isophote fitting

We present in Fig. 2 the host galaxy within 5″ from the quasar. At the largest spatial scale in Fig. 1, we can confirm the existence of two components: the outer component (OC) and the inner component (IC), as previously identified in isophotic contours in Martel et al. (2003) with HST/ACS. The OC is more extended on the north-east side than the IC in STIS wavelengths, with both centers offset from the central quasar, as reported in Martel et al. (2003). With STIS coronagraphy in Fig. 2, we confirm the existence of the E1, jet component (JC), and inner jet (IJ) components identified in Martel et al. (2003).

Fig. 2. Host galaxy of 3C 273 within 5″. (a) contains original data. (b) is the isophote model. (c) and (d) are isophote-removed data. We have newly identified a symmetric core component (CC) within ∼1″ (marked with dash dotted circle), a core blob (CB) component at ∼1″ to the west of the quasar, a core jet (CJ) component, three smaller scale blobs (marked with dotted circles) at ∼4σ levels in comparison with their surrounding host galaxy, and filamentary structures at ∼3σ levels. We recover the Martel et al. (2003) findings including jet component (JC), inner jet (IJ), and E1 component. Note: we masked out a PSF feature (e.g., Grady et al. 2003) in (c) and (d) a shaded ellipse to the east of the center. The data used to create this figure are available at the CDS.

To reveal small-scale structures, we performed isophote modeling to enhance the visibility of structures of smaller spatial scale. We used the isophote package from photuils (Bradley et al. 2023), which implemented the Jedrzejewski (1987) method to iteratively fit elliptical isophotes to the 3C 273 host image. We present the isophote model in Fig. 2b and the model-removed residuals in Fig. 2c, where we annotated the structures. To better reveal these structures, we reperformed isophote fitting while excluding them to reduce the fitting bias. The final isophote model and residuals are shown in Fig. 2.

3.1.2. New close-in host structures

We identified a core blob (CB) component at ∼1″ to the west of the quasar in Fig. 2c. We also identified a core jet (CJ) component along the direction of the large scale jet. In addition, we detected small-scale blobs spanning from 2″ to 4″ from the quasar. These blobs are marked with dotted circles in Fig. 2. We also detected a more symmetric core component (CC) for the host galaxy of 3C 273 within ≈1″ in Fig. 2a. The surface brightness⁷ is (5 ± 2)×10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹ for CB, (14 ± 5)×10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹ for CC, and (3 ± 2)×10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹ for CJ. In comparison, the background is (1.3 ± 1.3)×10⁻²² erg s⁻¹ cm⁻² Å⁻¹, which is two orders of magnitude fainter than these new close-in structures. We present the measurements in Table B.1.

We marked the credible new structures in Fig. 2c. First, the CB component is not a random or color-mismatch residual, in the sense that it does not have a radially shaped elongation that suggests a color mismatch or telescope breathing (e.g., Schneider et al. 2014). Second, its globular morphology does not resemble the spurious residuals in the existing STIS data archive for circumstellar structures in Ren et al. (2017). Third, the CJ component is not only along the direction for the jet, but also of elliptical shape in the isophote model in Fig. 2b. Finally, the blobs (i.e., b₁, b₂, b₃) do not arise from random noise and they are ∼4σ above their surrounding host galaxy (blob surface brightness: ∼4 × 10⁻²¹ erg s⁻¹ cm⁻² Å⁻¹). They persist in the individual reduction results from different telescope rolls in Fig. B.4 and they are not behind the diffraction spikes regions, where random noise may manifest positive residuals. By fitting the surface brightness of CC with a bivariate normal distribution while ignoring these small-scale structures, we obtained a semi-major axis of $a_{\mathrm{CC}}=0\stackrel{″}{.}578\pm 0\stackrel{″}{.}005$ and a semi-minor axis of $b_{\mathrm{CC}}=0\stackrel{″}{.}477\pm 0\stackrel{″}{.}004$, along with a position angle of $56\stackrel{°}{.}2\pm 1\stackrel{°}{.}5$ for the major axis. The fitted center is located at $\mathrm{\Delta }\mathrm{RA}=0\stackrel{″}{.}031\pm 0\stackrel{″}{.}004$ and $\mathrm{\Delta }\mathrm{Dec}=-0\stackrel{″}{.}024\pm 0\stackrel{″}{.}003$, or $0\stackrel{″}{.}039\pm 0\stackrel{″}{.}004$ from the quasar, which is less offset from the quasar than IC and OC ($0\stackrel{″}{.}65$–$1\stackrel{″}{.}40$; Martel et al. 2003). The ellipticity is η_CC = (a_CC − b_CC)/a_CC = 0.175 ± 0.010, which is smaller than the η ≳ 0.3 measurements in regions further out (Martel et al. 2003, Fig. 4 therein), suggesting that the host galaxy is more symmetric when it is closer to the quasar.

We detected evidence of extended filamentary structures to the northeast, east and west of the galactic nucleus. These structures are visible in the residuals image (filament surface brightness: ∼2 × 10⁻²¹ erg s⁻¹ cm⁻² Å⁻¹) in Figs. 2c and d, obtained after computing and subtracting an elliptical isophote model of the underlying early-type galaxy. Candidate filaments extending ∼5–10 kpc to the Northeast are reminiscent of the emission line nebulae seen, for instance, in the brightest cluster galaxy (BCG) in the Phoenix cluster (McDonald et al. 2018). They are believed to be multiphase gas that condenses out of the circumgalactic medium and fuels further AGN feedback (Gaspari et al. 2015, 2018). Structures seen to the Eest several kpc from the galactic nucleus and several kpc to the west are also potentially consistent with this scenario, although either deeper or (ideally) spectroscopic follow-up observations will be necessary to characterize these structures. Further multi-band photometric or spectroscopic follow-up observations with JWST will characterize these structures and determine the role they play in the lifecycle of AGN feedback in 3C 273.

3.2. Jet motion

3C 273 has been imaged with STIS, with the unocculted 50CCD imaging configuration, on 2000 April 3 in HST GO-8233 (PI: S. Baum). We did not remove the PSF for this observation using the coronagraphic archive from Ren et al. (2017), since unocculted PSFs do not resemble the coronagraphic ones (Grady et al. 2003). Nevertheless, the jet is visible in both the 2006 and our 2022 observations, establishing a 7950–8027 day separation (or 22 years) for the apparent motion measurement.

We aligned the rectified non-coronagraphic observations of 3C 273 with STIS using the “X marks the spot” method. In comparison with the coronagraphic observations, there are no data quality extensions in the flat-fielded files in HST GO-8233; thus, we removed the cosmic ray noise using the Astro-SCRAPPY code (McCully et al. 2018), which implements the van Dokkum (2001) approach in astropy (Astropy Collaboration 2013). In comparison with Meyer et al. (2016), where two HST instruments (WFPC2/PC and ACS/WFC) were used to obtain jet motion, our study with an identical instrument permits motion measurements with less offsets from different instruments.

To measure the jet motion between the two epochs using two images with different quality, we adopted the concept of dummy variables to simultaneously fit elliptical morphology and offset. Specifically, for one jet component, we fit identical bivariate normal distribution to its data in two epochs while allowing for translation and rotation (see Appendix C for details).

We also applied the same procedure to measure the offset for background sources to correct for motion biases due to different observations. The field rotation between the two epochs is $0\stackrel{°}{.}088\pm 0.004$, or $0\stackrel{°}{.}0040\pm 0\stackrel{°}{.}0002$ yr⁻¹. This rotation rate is faster than the result reported in the Ward-Duong et al. (2022) study of $0\stackrel{°}{.}0031\pm 0\stackrel{°}{.}0001$ yr⁻¹, potentially due to our combination of coronagraphic images after rotation along the center of the quasar. However, our result confirms the trend of STIS long-term rotational evolution of the STIS CCD images.

By adopting a conversion factor of $\frac{8.9856c}{{\mathrm{mas}\mathrm{yr}}^{-1}}$ as in Meyer et al. (2016), where c is the speed of light, we calculated the apparent proper motion of the jet components along the jet axis, which is at a position angle of $-139\stackrel{°}{.}0\pm 1\stackrel{°}{.}9$ east of north based on the jet components. We present the motion rates in Fig. 3 and Table C.1 following the annotations in Marshall et al. (2001). The motion of the jet components along the jet direction is consistent with zero at < 2σ levels. We witness a potential trend that the motion is faster when it is further out. By fitting a linear relationship between the distance to the quasar r in arcsec and motion v in c, and excluding the nearby galaxy components (i.e., In1, In2, and Ex1; e.g., Meyer et al. 2016), we have $\frac{\mathit{v}}{c}=0.9_{-0.5}^{+0.5}\left(\frac{r}{1^{″}}\right)-{15}_{-8}^{+8}$.

Fig. 3. Motion of different components v, along the jet direction r in units of speed of light c. The best-fit of $\frac{\mathit{v}}{c}=0.9_{-0.5}^{+0.5}\left(\frac{r}{1^{″}}\right)-{15}_{-8}^{+8}$ in Sect. 3.2 suggests faster outwards motion when the jet components are further from the quasar. The gray lines are 200 random samples from the best-fit parameters and their covariance.

4. Summary

In this study, we take advantage the coronagraphic imaging capabilities of HST/STIS in IWA and instrumentation stability (Debes et al. 2019) to validate STIS coronagraphy using extragalactic observations. By applying STIS coronagraph to the iconic quasar 3C 273 and capturing its coronagraphic PSF using two color-matching close-in reference stars, we have revealed its surrounding regions, including its host galaxy, from $\sim 0\stackrel{″}{.}2$ to ∼60″.

We detected a more symmetric core component, CC, for the host galaxy of 3C 273, in addition to confirming the existing large-scale asymmetric components IC and OC that were previously identified in HST/ACS coronagraphy from Martel et al. (2003). With the STIS coronagraphic observations, we also identified a core blob (CB) component, as well as other point-source-like objects, after removing isophotes from the host galaxy. The nature of the newly identified components, as well as the point-source-like objects, would require observations from other telescopes for further study.

Using a 22 yr timeline, we have constrained the apparent motion for the jet components of 3C 273. With an identical instrument between the two observations, we find that the jet components in the largest scale jet are likely to have faster motion when they are further from the quasar. We also confirm the long-term rotational trend of the STIS CCD images.

To characterize the observed host components for 3C 273, follow-up efforts in the NIR will help constrain the nature of them in multi-wavelength color and/or spectroscopy. From the ground, with the visible magnitude of quasars beyond the capability of most of coronagraphic imaging systems that perform adaptive optics correction in visible light, the use of pyramid wavefront sensing in NIR wavelengths could be necessary (Keck/NIRC2: Bond et al. 2020; VLT/ERIS: Davies et al. 2018; Kravchenko et al. 2022). From space, JWST coronagraphic and non-coronagraphic imaging efforts will also access such hosts, along with with the coronagraphic imaging observations needed for probing deeper into the quasar surroundings.

Given it is the only operating space-based coronagraph in visible light, HST/STIS offers an IWA of $\sim 0\stackrel{″}{.}2$ to image the surrounding environment around bright central sources. In addition to the carefully selected PSF stars for this study, contributing more reference star images taken at the smallest IWA positions after the commissioning of the BAR5 occulter for STIS (Schneider et al. 2014) to compile a library of PSF images (Ren et al. 2017) could help to reduce the impact of color differences and telescope status. Moreover, to reveal faint and extended structures that could reach the sensitivity limit of the STIS CCD, we recommend correcting its periodic readout variations (temporal variation: Jansen et al. 2003; Jansen 2013) before image rectification. Moving forward, an ongoing calibration program, HST GO-17135 for STIS coronagraphy, has adopted a carefully designed dithering strategy might provide smaller IWAs, which is expected to potentially push towards imaging the components that are behind the current coronagraphic occulters for stars and quasars towards $\sim 0\stackrel{″}{.}1$ (i.e., ∼0.3 kpc for 3C 273). Similarly, for JWST, such attempts are also necessary given the supported sizes of the coronagraphs are significantly larger than $0\stackrel{″}{.}1$ (e.g., JWST GO-3087). With smaller IWAs for both telescopes, we can both confirm the existence of closest-in components and constrain their physical properties from multi-band imaging. In high-energy observations, we can better characterize such structures, as well as binary active galactic nuclei (e.g., Pfeifle et al. 2023) when advanced coronagraphic instruments become available in the future.

Acknowledgments

We thank the anonymous referee for constructive suggestions that improved this Letter. B.B.R. would like to acknowledge the pioneering leadership of the late M. Schmidt, whose office at the Cahill Center in Caltech was assigned to B.B.R. in 2019, for inspiring future generations with his groundbreaking discovery of 3C 273. We thank Dean Hines, Anand Sivaramakrishnan, and Hsiang-Chih Hwang for discussions in the initial preparation of the observation proposal. We thank Paul Kalas, Tom Esposito, Yuguang Chen, and Zhi-Xiang Zhang for discussions. Based on observations with the NASA/ESA Hubble Space Telescope obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for program number HST GO-16715 was provided through a grant from the STScI under NASA contract NAS5-26555. This project has received funding from the European Union’s Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 101103114. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (PROTOPLANETS, grant agreement No. 101002188). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of Photutils, an Astropy package for detection and photometry of astronomical sources (Bradley et al. 2023).

References

Appendix A: HST/STIS coronagraphic imaging

HST/STIS coronagraphic imaging relies on well-chosen PSF templates for data reduction (e.g., Debes et al. 2019). We list the target and the PSFs in this study, along with their information and observation log, in Table A.1. We present the coronagraphic imaging occulting locations, as well as the coronagraphic data before and after PSF-removal, in Fig. A.1.

Fig. A.1. HST/STIS coronagraphic imaging can reveal faint structures that are overwhelmed by the PSF of central sources. (a) STIS coronagraphic occulter from Debes et al. (2017), with empirical occulting locations in Ren et al. (2017). (b) and (d) are example 3C 273 exposures with a few CTI trail examples marked using triangles from the BAR5 and WEDGEA0.6 occulters, respectively. (c) and (e) are the PSF-removed results on rectified CTI-corrected images for (b) and (d).

The coronagraphic imaging mode of STIS offers four major occulters: the two perpendicular wedges (WedgeA and WedgeB), BAR10, and BAR5 which offers the smallest inner working angle of $\sim 0\stackrel{″}{.}2$ (e.g., Medallon & Welty 2023). There are 13 default occulting locations, as shown in Fig. A.1 (a) for their empirical centers from Ren et al. (2017). In addition, STIS users can issue instrument offsets from these centers for customized imaging (i.e., “POSTARG”). Among the 13 occulting locations, WedgeA1.0 (IWA $\approx 0\stackrel{″}{.}5$) and WedgeA0.6 (IWA $\approx 0\stackrel{″}{.}3$), have mostly been used in existing locations (Ren et al. 2017), thus offering the largest PSF template archives for PSF removal.

BAR5 is oriented ∼12° clockwise from a horizontal occulter, and we can image the surrounding environments of central sources down to $\sim 0\stackrel{″}{.}2$. In fact, due to occulter deformation pre-launch, the BAR5 occulter was not enabled for HST General Observers science until Schneider et al. (2017). In comparison, WedgeA is positioned vertically with WedgeA0.6 offering a $\sim 0\stackrel{″}{.}3$ IWA. With BAR5 additionally blocking nearly horizontal regions and WedgeA0.6 blocking vertical regions due to the existence of the coronagraphic occulters, the combination of WedgeA0.6 and BAR5 is thus necessary to offer a nearly 360° angular coverage (with an exception for the diagonal and off-diagonal diffraction spikes) with the smallest IWAs. In addition, WedgeA0.6 is located 68 pixel or $3\stackrel{″}{.}45$ from the bottom of the coronagraphic edge, while, on the other hand, BAR5 is located 54 pixel or $2\stackrel{″}{.}74$ from the right coronagraphic edge; thus multiple telescope rolls can help in imaging structures beyond these angular radii.

To demonstrate the improvement in PSF removal with coronagraphic imaging, we present in Fig. A.1(b) and (d) two exposures on the detector frame at different roll angles using BAR5 and WedgeA0.6, respectively. With CTI effect corrected, Fig. A.1(c) and (e) are the corresponding 3C 273 surroundings after PSF removal. The count rates can be reduced by a factor of ∼5 at $\sim 1\stackrel{″}{.}5$ after PSF removal, thus extracting the surroundings that are overwhelmed by the PSF. In addition, the count rate reduction is performed on the PSF halo of the central source behind the coronagraph, making another step forward in actual suppression of the central source light using a combination of instrumentation and data reduction. To validate the PSF removal results for 3C 273, we reduced the PSF1 exposures using the PSF2 exposures for a comparison by following an identical PSF removal process as 3C 273. In the PSF-removed residuals for PSF1 in Fig. A.2, no significant signals resembling the 3C 273 residuals exist.

Fig. A.2. PSF removal for PSF1 using PSF2. (a) and (c) are rectified CTI-removed PSF1 exposures under BAR5 and WedgeA0.6, respectively. (b) and (d) are the corresponding PSF-removed residuals for PSF1 for (a) and (c). In comparison with the 3C 273 residuals in Fig. A.1 with identical color bars, there are no significant residuals for PSF1 that resemble 3C 273 residuals.

To obtain the final image in Fig. 1, multiple telescope rolls are necessary in obtaining a ∼360° angular coverage to resolve the constraints from the occulting locations being located near the edge of the STIS field of view and diffraction spikes. In addition, STIS allows subarray readouts in obtaining smaller field of view (e.g., Schneider et al. 2018) to increase observation efficiency in reducing readout time in its electronics. We did not request the subarray readout here, since the purpose is to image the 3C 273 surroundings to the largest spatial extent.

In the post-processing of the STIS coronagraphic imaging data in Sect. 2.3.2, we used the median exposures from two stars as the empirical PSF template for 3C 273. Alternatively, the usage of archival data may better capture the PSFs of central sources to extract extended structures (e.g., Soummer et al. 2014; Ren et al. 2018; Sanghi et al. 2022; Xie et al. 2023). To explore this, we updated the coronagraphic PSF archive at the BAR5 location for STIS from Ren et al. (2017). With both the principal-component-analysis-based PSF modeling approach (Soummer et al. 2012) and the non-negative matrix factorization method (Ren et al. 2018) applied to STIS coronagraphic imaging (e.g., Ren et al. 2017; ren18; walker21), we did not observe any significant improvement in the data reduction quality. This is due to the fact that BAR5 was supported relatively late (Schneider et al. 2017; Debes et al. 2019) since the installation of STIS on HST in 1997, and its PSF diversity is not comparable to that of the other STIS occulters (e.g., WedgeA0.6, WedgeA1.0: Ren et al. 2017). Nevertheless, given that it offers the narrowest IWA in visible wavelengths (Debes et al. 2019), the increase of popularity in BAR5 (e.g., Schneider et al. 2018; Walker et al. 2021; Ren et al. 2023; Debes et al. 2023) can help provide better PSF templates in the future. Its potential in PSF modeling can be maximized when a PSF reference does not match that of the target in complicated scenarios (e.g., source-variability-induced color and thus PSF change in STIS: Stark et al. 2023).

Appendix B: Auxiliary images

B.1. Phase II planning and coverage map

To enable HST users in phase II for realistic observation coverage planning with the APT for STIS coronagraphy, we have developed the VISIT-STIS-Coron visibility tool (Ren 2022). In the visibility tool, the parameters are consistent with the APT conventions, including occulting location, telescope parameters (e.g., ORIENT,⁸ POSTARG), etc. In addition, the occulting locations are empirically measured using the entire STIS coronagraphic archive in Ren et al. (2017), which can ensure maximum telescope pointing repeatability between observation planning and execution. For STIS coronagraphic data reduction, VISIT-STIS-Coron also includes a coronagraphic mask available in FITS format created in Debes et al. (2017).

Using on-sky observations, we created a coverage map for GO-16715 in Fig. B.1, which also takes into account of the diffraction spikes of the central quasar from HST optics. Despite a nearly 360° azimuthal coverage, the northeast region of the quasar surroundings has less exposures, which was not originally planned in the Phase II, but instead due to an update³ in the HST guide star catalog in 2022 between the two sets of 3C 273 visits.

Fig. B.1. On-sky coverage of 3C 273 surroundings using the STIS coronagraph from 36 individual readouts in six orbits. The colored values denote the total readout count for a specific location. The STIS coronagraphic occulters (e.g., Debes et al. 2019) are annotatated with arrows, and regions with zero coverage due to occulting are colored white. See Fig. 1 for the corresponding quasar surroundings with a same field of view.

B.2. S/N map

We calculated the pixel-wise S/N map for Fig. 1 as follows. First we identified four background regions that do not contain astrophysical signals by eye (see the full-frame result in Fig. B.2). We then calculated the standard deviation of the selected regions to estimate the noise. The noise level does not change significantly when we use only a few regions for analysis or change their locations and sizes. Finally, we divided the data in Fig. 1 by that calculated value of standard deviation for these background regions. We present the S/N map in Fig. B.3.

Fig. B.2. Full-frame result of 3C 273 surroundings from GO-16715. The regions within the four dashed circles are used to estimate the noise.

Fig. B.3. Pixel-wise S/N map, obtained through dividing Fig. 1 by the pixel-wise standard deviation of the astrophysical-signal-free regions identified by eye.

To investigate the impact of the selection of background regions, we changed the locations and sizes selected regions for noise estimation and the resulting S/N maps do not have significant change from Fig. B.3. In fact, due to the positions of the BAR5 and WedgeA0.6 occulters, the telescope roll angles, and our readout in the full frame, the final full-frame result has a field of view of ∼100″ × 100″. In comparison, the 3C 273 host resides within a 40″ × 40″ region in Fig. 1 which is cropped from the full-frame result. As a result, while there are foreground stars (2MASS J12290318+0203185 and 2MASS J12290786+0203359, which are located at 644 ± 7 pc and 203.4 ± 1.9 pc from the Sun in Gaia DR3, respectively. The two are not shown in Fig. 1, but they are available in the corresponding full-frame image in Fig. B.2) and background galaxies in the full-frame result, which can reach a large fraction of the areas in a 120″ × 120″ field; the majority of the full-frame result contains background regions that can be used to quantify the noise in the reduction.

B.3. Small-scale features

In Fig. B.4, we present the small-scale features identified in Fig. 2 at different telescope orientations and occulting locations. In Table B.1, we present the surface brightness for the newly identified features, as well as those in Martel et al. (2003) from Fig. 2.

Fig. B.4. Small-scale features in isophote-removed residuals in Fig. 2 can persist in multiple telescope orientations. (a), (c), and (e) are for BAR5. (b), (d), and (f) are for WedgeA0.6.

Appendix C: Jet motion measurement

We assume that the surface brightness distribution of an elliptical component follows a bivariate normal distribution,

$$\begin{array}{c}\hfill S({\mathit{r}}_0)\sim {\mathcal{N}}_2({\mathit{\mu }}_{\mathbf{0}},{\mathbf{\Sigma }}_{\mathbf{0}}),\end{array}$$(C.1)

where 𝒩_k denotes a normal distribution with a dimension of k ∈ ℤ and r₀ = [x, y]^T ∈ ℝ^2 × 1 denotes the on-sky location, while μ₀ ∈ ℝ^2 × 1 and Σ₀ ∈ ℝ^2 × 2 are the expectation and covariance matrix for the distribution, respectively. We can perform matrix translation then rotation to obtain a new surface brightness distribution, S′.

C.1. Elliptical component motion

To enable the translation and rotation of the surface brightness distribution in Eq. (C.1), we define $\mathit{r}={[{\mathit{r}}_0^{\top },1]}^{\top }$ to be a 3 × 1 column matrix, where ^⊤ denotes matrix transpose. We additionally define its corresponding expectation and covariance matrices to be:

$$\begin{array}{c}\hfill \mathit{\mu }={\mathit{\mu }}_0\oplus 1,\end{array}$$(C.2)

and

$$\begin{array}{c}\hfill \mathbf{\Sigma }={\mathbf{\Sigma }}_0\oplus 0,\end{array}$$(C.3)

respectively, where ⊕ is matrix direct sum. With these, we can rewrite Eq. (C.1) in a 3D form,

$$\begin{array}{c}\hfill S(\mathit{r})\sim {\mathcal{N}}_3(\mathit{\mu },\mathbf{\Sigma }).\end{array}$$(C.4)

For a translation matrix T, we have

$$\begin{array}{c}\hfill T=\left[\begin{array}{ccc}1& 0& t_x\\ 0& 1& t_y\\ 0& 0& 1\\ \end{array}\right],\end{array}$$

where t_x ∈ ℝ and t_y ∈ ℝ denote the translation along the x-direction and y-direction, respectively. For a rotation matrix R, we have

$$\begin{array}{c}\hfill R=\left[\begin{array}{ccc}cos{t}_{\theta }& -sin{t}_{\theta }& 0\\ sin{t}_{\theta }& cos{t}_{\theta }& 0\\ 0& 0& 1\\ \end{array}\right],\end{array}$$

where t_θ ∈ [ − π, π) denotes the clockwise rotation about the origin. A translation then rotation a surface brightness distribution following Eq. (C.4) then follows

$$\begin{array}{c}\hfill S^{\prime }(\mathit{r})=RTS(\mathit{r})\sim {\mathcal{N}}_3(RT\mathit{\mu },RT\mathbf{\Sigma }{T}^{\top }{R}^{\top }).\end{array}$$

Given that the last row and column of Eq. (C.3) are composed of 0, we can rewrite the above distribution as

$$\begin{array}{c}\hfill S^{\prime }(\mathit{r})\sim {\mathcal{N}}_3(RT\mathit{\mu },R\mathbf{\Sigma }{R}^{\top }).\end{array}$$(C.5)

C.2. Motion quantification

For a bivariate normal distribution at two epochs, we assume that its morphology does not change significantly. In this way, we can use Eqs. (C.4) and (C.5), which share identical morphological parameters μ and Σ, to obtain the relative motion between two epochs. To maximize the information from both datasets, we use the statistical concept of dummy variables [see, e.g.,][for an application in high-contrast imaging science]ren20 to fit Eqs. (C.4) and (C.5) simultaneously.

First, we rewrite the distribution of the first two entries in Eqs. (C.4) in a general form in Cartesian coordinates, with the distribution centered at (x₀, y₀) and rotated θ₀ ∈ [ − π, π) along the counterclockwise rotation from a rectangular bivariate normal distribution. We have

$$\begin{array}{c}\hfill S(x,y)=Aexp\{-\frac{1}{2}[a{(x-x_0)}^2+2b(x-x_0)(y-y_0)+c{(y-y_0)}^2]\},\end{array}$$(C.6)

where A ∈ ℝ is the amplitude, and

$$\begin{array}{c}\hfill \{\begin{array}{cc}a& =\frac{{cos}^2{\theta }_0}{{\sigma }_x^2}+\frac{{sin}^2{\theta }_0}{{\sigma }_y^2}\\ b& =-\frac{sin2{\theta }_0}{2{\sigma }_x^2}+\frac{sin2{\theta }_0}{2{\sigma }_y^2}\\ c& =\frac{{sin}^2{\theta }_0}{{\sigma }_x^2}+\frac{{cos}^2{\theta }_0}{{\sigma }_y^2}\end{array}.\end{array}$$(C.7)

Second, we introduce a dummy variable D ∈ {0, 1} to Eq. (C.6) to denote the offset along x- and y-axis, and rotation from θ₀ using the following substitution:

$$\begin{array}{c}\hfill \{\begin{array}{c}{x}_0\equiv x_0+t_{x}D\\ y_0\equiv y_0+t_{y}D\\ {\theta }_0\equiv {\theta }_0+t_{\theta }D\end{array},\end{array}$$(C.8)

to obtain the general form for the first two entries in Eq. (C.5) when D = 1.

Third, by combining Eqs. (C.6)–(C.8), we can obtain the general form for the first two entries in Eq. (C.4), when the dummy variable is D = 0, and the two in Eq. (C.5), when D = 1. To obtain the morphological and offset parameters in two epochs, we use (x, y, D) as the independent variables and surface brightness, S(x, y, D), as the dependent variable.

For each elliptical component, using the curve_fit function from scipy (Virtanen et al. 2020), we can obtain the best-fit parameters for the elliptical morphology (i.e., A, x₀, y₀, σ_x, σ_y, θ₀) and the motion (i.e., t_x, t_y, t_θ), as well as their covariance matrix. In practice, we also had a local background brightness value that was added to Eq. (C.6) to quantify the background difference in two epochs. We adopt the square root of the diagonal values of the covariance matrix from curve_fit as the uncertainty of the elliptical and motion parameters. In Fig. C.1, we present a demonstration of the fitting process.

Fig. C.1. Ellipse fitting to obtain the morphology and motion for jet component A2. The (a) 2000 and (b) 2022 data are normalized to have a peak count of one in each panel. Then, (c) and (d) are the best-fit models, which follow identical bivariate normal distribution but with location offset, for (a) and (b), respectively. Finally, (e) and (f) are the residuals after subtracting the models from the corresponding data.

C.3. Instrument offset calibration

To address any potential bias in our alignment and rotation of two images in 2000 and 2022, we used background objects to calibrate the translation and rotation offsets between two epochs (see Fig. C.2). Assuming the background sources in the two images do not move, we can first offset then rotate one image to another to calibrate the global offset.

Fig. C.2. Components (marked by squares) and background objects (marked with diamonds) used in motion measurement and instrument offset calibration, respectively.

For a background object, we can follow Appendix C.2 to obtain its best-fit and uncertainty values for the motion parameters: t_x, t_y, and t_θ. To quantify the global offset in ${\widehat{t}}_x$, ${\widehat{t}}_y$, and ${\widehat{t}}_{\theta }$, we assume that the motion parameters are dependent on the centers of the ellipses (x₀, y₀). Using the odr function which performs orthogonal least-squared-fitting for data with both input and output uncertainties (Boggs et al. 1989) from scipy, we obtained the best-fit parameters for ${\widehat{t}}_x$, ${\widehat{t}}_y$, and ${\widehat{t}}_{\theta }$.

For a jet component, we followed Appendix C.2 to obtain its offset, then corrected the global offset measured from background objects. To obtain the jet component motion along the jet direction, we first obtained the position angle of the jet by performing least-square fit to the location of the jet components. We then projected the calibrated offset for a jet component, which is expressed as a bivariate normal distribution from the odr outputs, to the jet direction though matrix rotation, see Eq. (B8) in Shuai et al. (2022) for a similar approach.

C.4. Measured motion

We present the measured motion rates from Sect. 3.2 in Table C.1.

All Tables

All Figures

	Fig. 1. 3C 273 host galaxy and surroundings in visible light seen with the HST/STIS coronagraph. The surface brightness is in log scale, and the central dark regions are not accessible due to coronagraph blockage. The data used to create this figure are available at the CDS.
In the text

Fig. 2. Host galaxy of 3C 273 within 5″. (a) contains original data. (b) is the isophote model. (c) and (d) are isophote-removed data. We have newly identified a symmetric core component (CC) within ∼1″ (marked with dash dotted circle), a core blob (CB) component at ∼1″ to the west of the quasar, a core jet (CJ) component, three smaller scale blobs (marked with dotted circles) at ∼4σ levels in comparison with their surrounding host galaxy, and filamentary structures at ∼3σ levels. We recover the Martel et al. (2003) findings including jet component (JC), inner jet (IJ), and E1 component. Note: we masked out a PSF feature (e.g., Grady et al. 2003) in (c) and (d) a shaded ellipse to the east of the center. The data used to create this figure are available at the CDS.

In the text

	Fig. 3. Motion of different components v, along the jet direction r in units of speed of light c. The best-fit of $\frac{\mathit{v}}{c}=0.9_{-0.5}^{+0.5}\left(\frac{r}{1^{″}}\right)-{15}_{-8}^{+8}$ in Sect. 3.2 suggests faster outwards motion when the jet components are further from the quasar. The gray lines are 200 random samples from the best-fit parameters and their covariance.
In the text

Fig. A.1. HST/STIS coronagraphic imaging can reveal faint structures that are overwhelmed by the PSF of central sources. (a) STIS coronagraphic occulter from Debes et al. (2017), with empirical occulting locations in Ren et al. (2017). (b) and (d) are example 3C 273 exposures with a few CTI trail examples marked using triangles from the BAR5 and WEDGEA0.6 occulters, respectively. (c) and (e) are the PSF-removed results on rectified CTI-corrected images for (b) and (d).

In the text

	Fig. A.2. PSF removal for PSF1 using PSF2. (a) and (c) are rectified CTI-removed PSF1 exposures under BAR5 and WedgeA0.6, respectively. (b) and (d) are the corresponding PSF-removed residuals for PSF1 for (a) and (c). In comparison with the 3C 273 residuals in Fig. A.1 with identical color bars, there are no significant residuals for PSF1 that resemble 3C 273 residuals.
In the text

Fig. B.1. On-sky coverage of 3C 273 surroundings using the STIS coronagraph from 36 individual readouts in six orbits. The colored values denote the total readout count for a specific location. The STIS coronagraphic occulters (e.g., Debes et al. 2019) are annotatated with arrows, and regions with zero coverage due to occulting are colored white. See Fig. 1 for the corresponding quasar surroundings with a same field of view.

In the text

	Fig. B.2. Full-frame result of 3C 273 surroundings from GO-16715. The regions within the four dashed circles are used to estimate the noise.
In the text

	Fig. B.3. Pixel-wise S/N map, obtained through dividing Fig. 1 by the pixel-wise standard deviation of the astrophysical-signal-free regions identified by eye.
In the text

	Fig. B.4. Small-scale features in isophote-removed residuals in Fig. 2 can persist in multiple telescope orientations. (a), (c), and (e) are for BAR5. (b), (d), and (f) are for WedgeA0.6.
In the text

Fig. C.1. Ellipse fitting to obtain the morphology and motion for jet component A2. The (a) 2000 and (b) 2022 data are normalized to have a peak count of one in each panel. Then, (c) and (d) are the best-fit models, which follow identical bivariate normal distribution but with location offset, for (a) and (b), respectively. Finally, (e) and (f) are the residuals after subtracting the models from the corresponding data.

In the text

	Fig. C.2. Components (marked by squares) and background objects (marked with diamonds) used in motion measurement and instrument offset calibration, respectively.
In the text