Regulation of accretion by its outflow in a symbiotic star: the 2016 outflow fast state of MWC 560

TL;DR

This study investigates how outflows influence accretion discs in the symbiotic star MWC 560 during its 2016 fast outflow state, revealing stable disc conditions and outflow-disc interactions similar to those in X-ray binaries.

Contribution

It provides multi-wavelength observations of the 2016 outflow fast state, showing stable accretion disc properties and the role of outflows in regulating accretion in a symbiotic binary.

Findings

Outflow power increased abruptly during the 2016 state.

High-velocity and low-velocity outflows coexist, producing soft X-ray emission.

Stable accretion disc persisted despite high accretion rates.

Abstract

How are accretion discs affected by their outflows? To address this question for white dwarfs accreting from cool giants, we performed optical, radio, X-ray, and ultraviolet observations of the outflow-driving symbiotic star MWC 560 (=V694 Mon) during its 2016 optical high state. We tracked multi-wavelength changes that signalled an abrupt increase in outflow power at the initiation of a months-long outflow fast state, just as the optical flux peaked: (1) an abrupt doubling of Balmer absorption velocities; (2) the onset of a $20$ $\mu$Jy/month increase in radio flux; and (3) an order-of-magnitude increase in soft X-ray flux. Juxtaposing to prior X-ray observations and their coeval optical spectra, we infer that both high-velocity and low-velocity optical outflow components must be simultaneously present to yield a large soft X-ray flux, which may originate in shocks where these fast and slow absorbers collide. Our optical and ultraviolet spectra indicate that the broad absorption-line gas was fast, stable, and dense ($\gtrsim10^{6.5}$ cm$^{-3}$) throughout the 2016 outflow fast state, steadily feeding a lower-density ($\lesssim10^{5.5}$ cm$^{-3}$) region of radio-emitting gas. Persistent optical and ultraviolet flickering indicate that the accretion disc remained intact. The stability of these properties in 2016 contrasts to their instability during MWC 560's 1990 outburst, even though the disc reached a similar accretion rate. We propose that the self-regulatory effect of a steady fast outflow from the disc in 2016 prevented a catastrophic ejection of the inner disc. This behaviour in a symbiotic binary resembles disc/outflow relationships governing accretion state changes in X-ray binaries.

Figures

Tables

Table 1. Table 1: VLA observations and results. All detections were unresolved point sources.

Date a	Config.b	Band	$\overline{\nu }$ (GHz) c	On-source time (min)	Flux density ($\mu$Jy)	${\sigma }_{rand}$ ($\mu$Jy) d	${\sigma }_{tot}$ ($\mu$Jy) e
2014 Oct 02	DnC	X	9.83	88	37	3	4
2016 Apr 04	C	X	9.84	68	85	4	6
2016 May 01	CnB	X	9.84	66	105f	19f	20
2016 May 24	B	X	9.86	67	136	4	8
2016 Jul 29	B	X	9.82	40	175	5	10
2016 Oct 17	A	X	9.85	40	144	4	9
2016 Nov 09	A	X	9.85	45	138	4	8
2017 Jan 18	A	X	9.80	42	120	4	7
2016 May 24	B	S	3.10	67	128	5	8
2016 Jul 29	B	S	3.08	38	163	7	11
2016 Oct 17	A	S	3.07	38	162	7	11
2016 Nov 09	A	S	3.08	37	135	7	10
2017 Jan 18	A	S	3.11	37	112	6	8
2016 Jul 29	B	Ka	33.07	22	183	16	24
2016 Oct 12	A	Ka	33.07	21	153	10	19
2016 Nov 09	A	Ka	33.09	26	163	10	19
2017 Jan 18	A	Ka	32.96	24	128	11	17

Table 2. Table 2: Radio spectral indices. Spectral index is defined such that F ∝ ν α proportional-to 𝐹 superscript 𝜈 𝛼 F\propto\nu^{\alpha} between the two given bands; approximate central frequencies of the bands are listed for convenience. σ α subscript 𝜎 𝛼 \sigma_{\alpha} is the uncertainty in the spectral index α 𝛼 \alpha propagated from the total uncertainty σ t ​ o ​ t subscript 𝜎 𝑡 𝑜 𝑡 \sigma_{tot} in flux from Table 1 .

Bands	Epoch	$\alpha$	${\sigma }_{\alpha }$
SX (3.1 to 9.8 GHz)	2016 May 24	0.05	0.07
	2016 Jul 29	0.06	0.07
	2016 Oct 12	-0.10	0.07
	2016 Nov 09	0.02	0.08
	2017 Jan 18	0.06	0.08
SKa (3.1 to 33 GHz)	2016 Jul 29	0.05	0.06
	2016 Oct 12,17a	-0.03	0.06
	2016 Nov 09	0.08	0.06
	2017 Jan 18	0.05	0.06
XKa (9.8 to 33 GHz)	2016 Jul 29	0.04	0.12
	2016 Oct 12,17a	0.05	0.11
	2016 Nov 09	0.14	0.11
	2017 Jan 18	0.05	0.12

Table 3. Table 3: X-ray model for Chandra 2016 March 8-9 epoch. Both components with temperature kT and normalization norm were modelled using the CIAO Sherpa clones of the XSPEC apec diffuse collisionally-ionized plasma model with solar abundance ( angr ) and redshift 0, each absorbed by their own wabs model with hydrogen column nH.

Component	Parameter	Best fit	-2$\sigma$	+2$\sigma$
soft	kT (keV)	0.77	-0.66	+0.22
	norm	2.0E-05	-1.2E-05	+0.036
	nH (1022 cm-2)	0.46	-0.33	+0.64
hard	kT (keV)	11.26	frozen
	norm	8.9E-05	-6.7E-05	+1.39E-4
	nH (1022 cm-2)	13.48	-12.74	+19.39

Table 4. Table 4: Physical properties of the outflow fast state (Jan-Jul 2016).

Density of absorption-line region:	$\gtrsim$106.5 cm-3
Column density of outflow:	$\gtrsim$1023 cm-2
Maximum outflow velocity:	2500–3000 km s-1
Strong-shock velocity for soft X-rays:	300–900 km s-1
Radio emission mechanism:	optically-thin thermal
Density of radio-emitting region:	$\lesssim$105.5 (2.5 kpc/d)1/2 cm-3
Outflow radius for radio emissions:	$\gtrsim$60 (d/2.5 kpc) au
Mass-outflow rate for radio emissions:*	$\gtrsim$10-6 (d/2.5 kpc)5/2 M⊙ yr-1
Inner accretion disc:	intact with persistent optical/UV flickering
Accretion-disc bolometric luminosity peak:	1800 (d/2.5 kpc)2 L⊙
Accretion rate peak:$†$	6$\times$10-7 (d/2.5 kpc)2 (R / 0.01 R⊙) (0.9 M⊙/ M) M⊙ yr-1

Table 5. Table 5: Swift UV spectra. Observation metadata and quality flags for Swift UVOT UV Grism spectra.

Start time (UT)	Exposure time (s)	Roll angle ($deg$)	Quality flags*
2016 Mar 02 07:47	785.2	251	ZOc:2300-2376,2590-2650
2016 Mar 05 12:09	892.5	254	ZOc:2790-2890,2880-3200
2016 Mar 09 12:00	307.8	251	ZOc:2300,2376,2590-2650
2016 Mar 12 08:39	463.3	254	ZOc:2790-2890,2880-3200”
2016 Mar 17 22:25	892.5	256	ZOc:2285-2385
2016 Mar 20 22:19	999.8	260	ZOc:2720-2820
2016 Mar 23 22:03	470.1	260	ZOc:2740-2810
2016 Mar 26 17:00	1008.6	264	ZOc:1886
2016 Mar 29 10:25	946.6	268	ZOc:1775,1905-2005,3230-3330; bright FO nearby
2016 Apr 01 15:00	892.5	264	ZOc:1870,FOc:$>$3700
2016 Apr 04 14:49	892.5	268	ZOc:1720,1930-2000,2180-2255,2300-2417,2547-2600,3240-3325
2016 Apr 20 09:03	892.5	281	ZOc:1756,2400-2440w,2540-2600,2935-3110
2016 Apr 27 05:06	892.5	283	ZOc:1880-1947,2006-2090,2510-2640
2016 May 04 14:16	463.3	284	ZOc:2012-2070,2525-2622,3025-3127; bright FO nearby
2016 May 11 04:02	892.5	286	ZOc:1970-2010,2250-2300; fluxes underestimated: loss of lock
2016 May 18 17:54	892.5	297	ZOc:2620-2720
2016 Jun 01 23:17	845.9	307	ZOc:1990

Table 6. Table 6: Optical spectra. ECH = Echelle spectrograph on the Asiago 0.61m or 1.82m telescopes (R=18000-20000). B&C = B&C spectrograph on Asiago 1.22m telescope (R=1300). LT = Liverpool Telescope (R ∼ similar-to \sim 2500). UCS8 = home-built spectrograph in L’Aquila (R ∼ similar-to \sim 1300). PSO = Shelyak LHires III spectrograph in Tata (resolution varies and is specified in table). For more details, see § 2.1 .

Date	Time (UT)	Exposure (s)	Telescope+Instrument
2007 Oct 30	03:01	4500	0.61m+ECH
2012 Jan 27	21:40	900	1.22m+B&C
2012 Feb 08	22:10	600	1.82m+ECH
2012 Mar 31	18:52	480	1.22m+B&C
2013 Jan 26	22:13	900	1.22m+B&C
2014 Feb 09	21:22	480	1.22m+B&C
2014 Mar 14	20:08	600	1.82m+ECH
2014 Oct 30	02:30	900	1.22+B&C
2015 Mar 08	20:24	900	1.82m+ECH
2015 Mar 10	19:23	600	1.22m+B&C
2015 Oct 02	02:56	3620	PSO (R$\sim$450)
2015 Oct 04	02:41	4433	PSO (R$\sim$14000)
2015 Dec 31	00:37	4959	PSO (R$\sim$3300)
2016 Jan 21	21:23	300	1.22m+B&C
2016 Feb 05	17:49	4500	0.61m+ECH
2016 Feb 05	19:07	480	1.22m+B&C
2016 Feb 05	22:55	2800	USC8
2016 Feb 11	19:44	2800	USC8
2016 Feb 20	20:20	900	1.82m+ECH
2016 Feb 20	21:08	480	1.22m+B&C
2016 Feb 20	17:58	2800	USC8
2016 Feb 23	18:25	3600	0.61m+ECH
2016 Feb 25	18:32	3600	0.61m+ECH
2016 Mar 01	17:58	3600	0.61m+ECH
2016 Mar 02	19:52	3600	USC8
2016 Mar 03	18:54	4500	0.61m+ECH
2016 Mar 04	18:38	3200	USC8
2016 Mar 06	20:31	4500	0.61m+ECH
2016 Mar 06	18:10	3200	USC8
2016 Mar 07	18:15	3600	0.61m+ECH
2016 Mar 07	20:14	60	LT
2016 Mar 07	23:46	60	LT
2016 Mar 08	18:22	4500	0.61m+ECH
2016 Mar 08	20:25	3600	0.61m+ECH

Table 7. Table B2: – continued

2016 Mar 08	22:40	60	LT
2016 Mar 08	00:03	60	LT
2016 Mar 08	20:11	60	LT
2016 Mar 08	21:22	60	LT
2016 Mar 08	18:23	3600	USC8
2016 Mar 09	20:14	60	LT
2016 Mar 09	22:58	60	LT
2016 Mar 10	18:34	4500	0.61m+ECH
2016 Mar 10	18:15	480	1.22m+B&C
2016 Mar 11	18:43	5400	USC8
2016 Mar 13	18:13	360	1.22m+B&C
2016 Mar 14	21:32	1800	0.61m+ECH
2016 Mar 17	19:01	5400	0.61m+ECH
2016 Mar 18	19:03	5400	0.61m+ECH
2016 Mar 18	18:40	5400	USC8
2016 Mar 19	19:08	900	1.82m+ECH
2016 Mar 19	18:58	360	1.22m+B&C
2016 Mar 19	18:53	4800	USC8
2016 Mar 21	18:58	5400	0.61m+ECH
2016 Mar 23	18:42	4500	0.61m+ECH
2016 Mar 24	18:34	3600	0.61m+ECH
2016 Apr 03	19:13	300	1.22m+B&C
2016 Apr 06	19:10	3600	0.61m+ECH
2016 Apr 10	18:53	3600	0.61m+ECH
2016 Apr 12	18:54	480	1.22m+B&C
2016 Apr 14	19:39	3600	0.61m+ECH
2016 Apr 19	19:09	3600	0.61m+ECH
2016 Apr 19	19:10	360	1.22m+B&C
2016 Apr 27	19:24	3600	0.61m+ECH
2016 May 06	19:27	300	1.22m+B&C
2016 Oct 05	03:07	4500	0.61m+ECH
2016 Oct 13	03:41	900	1.82m+ECH
2016 Oct 23	03:22	1200	PSO (R$\sim$ 2500)
2016 Dec 04	01:15	600	PSO (R$\sim$450)
2016 Dec 11	01:38	3600	PSO (R$\sim$450)
2016 Dec 16	01:48	1200	1.82m+ECH
2017 Jan 08	20:43	600	PSO (R$\sim$450)
2017 Jan 20	20:43	1800	PSO (R$\sim$750)
2017 Jan 21	20:43	2400	PSO (R$\sim$750)
2017 Jan 22	20:43	600	PSO (R$\sim$750)

Equations

$$F_{\nu }\approx (T_e/1200K)(\nu /GHz{)}^2({\theta }_s/arcsec{)}^2\tau$$

$$\tau \approx 3.28\times 10^{-7}(T_e/10^{4}K{)}^{-1.35}(\nu /GHz{)}^{-2.1}(n_e^{2}s/c{m}^{-6}pc)$$

Full text

††pubyear: 2019††pagerange: Regulation of accretion by its outflow in a symbiotic star: the 2016 outflow fast state of MWC 560–G

1 Introduction

White dwarf (WD) symbiotic stars, hereafter symbiotics, are interacting binaries in which a WD accretes from a cool giant. They have larger accretion discs than their counterparts with main-sequence-like donors (cataclysmic variables: CVs), and several have spatially resolved jets. They are candidate progenitors for single-degenerate supernovae Ia (e.g., Cao et al., 2015; Munari & Renzini, 1992) and possible ancestors of double-degenerate supernovae Ia, raising the stakes for understanding their accretion discs, WD masses, and evolution. Accretion disc outflows may be fundamental to that investigation.

MWC 560 ($\equiv$V694 Mon) is a remarkable symbiotic in which the accretion disc drives a powerful outflow, producing broad, blue-shifted, variable absorption lines that extend up to thousands of km s*-1* from atomic transitions in the infrared (IR), optical, near-ultraviolet (NUV), and far-ultraviolet (FUV): He i, H i, Al iii, Mg ii, Fe ii, Cr ii, Si ii,C ii, Ca ii, Mg i, Na i, O i, C iv, Si iv, and N v (Bond et al., 1984; Tomov et al., 1990; Iijima, 2001; Michalitsianos et al., 1991; Tomov et al., 1992; Meier et al., 1996; Schmid et al., 2001; Goranskij et al., 2011; Lucy et al., 2018). It may be the prototype of a population of broad absorption line symbiotics (Lucy et al., 2018), and a useful nanoscale-mass analogy to quasars both in terms of its emission lines (Zamanov & Marziani, 2002) and its absorption lines (Lucy et al., 2018). The distance to the system is probably about 2.5 kpc (Appendix A). The red giant (RG) is mid-M and not a Mira; the infrared spectrum is S-type in that it is dominated by the RG rather than by dust (Meier et al., 1996). The WD mass is probably at least 0.9 M*⊙* (Zamanov et al. 2011a, Stute & Sahai 2009). The outflow may be a highly-collimated, baryon-loaded jet (Schmid et al., 2001) or a polar wind with a wider opening angle (Lucy et al., 2018). Zamanov et al. (2011a) have previously proposed that switching between different mass outflow regimes in MWC 560 is related in some way to switching between different regimes of mass inflow in the inner accretion disc, but the nature of that relationship is undetermined.

Here we describe coordinated multi-wavelength observations of MWC 560 conducted during 2016, prompted by the January 2016 peak of a year-long rise in optical flux (which may have been predicted by the system’s flux periodicities; Leibowitz & Formiggini 2015). Our 2016 observations are supplemented by data collected throughout the preceding decade. In § 2, we describe our observations and data reduction in the optical (§ 2.1), radio (§ 2.2), X-ray (§ 2.3), and NUV (§ 2.4) wavebands. We present our results in § 3, examining the optical absorption (§ 3.1), the radio flux (§ 3.2), the X-ray spectrum (§ 3.3), the NUV absorption (§ 3.4), and the optical and NUV flickering (§ 3.5). We cohere our results into a physical narrative in § 4, describing the state of the outflow (§ 4.1) and of the accretion disc (§ 4.2), then demonstrating a regulatory relationship between the disc and its outflow over the course of MWC 560’s history (§ 4.3). We summarize our conclusions in § 5.

Throughout this paper, we use heliocentric velocities, and scale to a distance of 2.5 kpc and a WD mass of 0.9 M*⊙* unless otherwise noted.

2 Observations and Data Reduction

2.1 Optical spectroscopy

We obtained 74 optical spectra of MWC 560111Including one spectrum presented previously in Munari et al. (2016)., listed in Appendix B and delineated below.

We obtained 30 echelle spectra with two telescopes at Asiago: the ANS Collaboration 0.61m operated in Varese by Schiaparelly Observatory, and the 1.82m operated by the National Institute of Astrophysics (INAF). The Varese 0.61m telescope used an Astrolight Instruments mark.III Multi Mode Spectrograph with a SBIG ST10XME CCD camera, covering 4225–8900 Å in 27 orders with a spectral resolution of R=18000, The Asiago 1.82m telescope used the REOSC Echelle spectrograph with an Andor DW436-BV and an E2V CCD42-40 AIMO back illuminated CCD, covering 3600–7300Å in 32 orders with R=20000.

The wavelength calibration for each echelle spectrum was obtained by exposing on a Thorium lamp before and after the science spectrum. A pruned list of about 800 unblended Thorium lines evenly distributed along and among the Echelle orders was fitted to the observed spectra, providing a wavelength solution with an rms of 0.005 and 0.009 Å ($\equiv$0.3 and 0.6 km s*-1*) for the Asiago and Varese spectrographs, respectively.

We also obtained 16 spectra with R$\approx$1300 (3300–8050 Å) at Asiago with the 1.22m telescope + B&C spectrograph, operated in Asiago by the University of Padova, using an ANDOR iDus DU440A with a back-illuminated E2V 42-10 sensor.

At all Asiago telescopes, the long slit, with a 2 arcsec slit width, was aligned with the parallactic angle. Spectrophotometric standards at similar airmass were observed immediately before and after MWC 560. The spectra were similarly reduced within IRAF, carefully involving all steps connected with bias correction, darks and flats, long-slit sky subtraction, and wavelength and flux calibration.

We obtained 10 spectra with R$\approx$1300 (3750-7350 Å) using a home-built slit-spectrograph mounted on an 8-in C8 Schmidt-Cassegrain telescope located in L’Aquila, Italy. The slit width was set to 3 arcsec and a fixed East-West orientation. The wavelength solution was obtain by fitting to an Ar-Ne comparison lamp, and observation of spectrophotometric standard stars was used to flux the spectra. Due to the consistency of the data reduction process for these spectra and the close match to other spectra when they were taken at similar times, we will hereafter group these data with the 1.22m telescope + B&C spectrograph spectra, which have the same resolving power.

We obtained 10 H$\alpha$ spectra of varying resolutions with a Shelyak LHires III spectrograph, using 150–2400 gr mm*-1* gratings and an ATIK 414 EXm camera, on 25 and 30 cm Newtonian telescopes located in Tata, Hungary. Slit widths varied from 15 to 35 microns. These spectra were reduced in ISIS following the standard procedure. The response calibrator star was fit with a 3rd order polynomial for spectra with R$\sim$3000 and above, and a manually selected spline for lower-resolution spectra only. The response calibration does not exhibit any features that would give rise to results discussed in this paper.

We obtained 8 spectra with the 2m fully robotic Liverpool Telescope (LT; Steele et al. 2004) using the FRODOSpec instrument (Barnsley et al., 2012). FRODOSpec was operated in its “low resolution” mode, obtaining spectra covering 3900–5700 Å and 5800–9400 Å at a resolution of $R\sim 2500$. Each spectra epoch consisted of $3\times 60$ s exposures. Wavelength calibration was performed by comparison to a Xe arc lamp and the data were reduced using the pipeline described in Barnsley et al. (2012).

2.2 Radio

In 2016, we observed for a total of 18 hours on the VLA through Projects 16A-448 and 16A-490, divided between 5 epochs in the S band (3.1 GHz), 7 epochs in the X band (9.8 GHz) and 4 epochs in the Ka band (33.1 GHz). Additionally, on 2014 October 2, we observed MWC 560 in the X band using about 2 hours of VLA Project 14B-394. More observational details are listed in Table 1.

We observed 3C147 for flux/gain calibration in all bands, QSO J0730-116 for phase calibration in S and X, J0724-0715 for phase calibration in Ka, and J0319+4130 for polarization leakage calibration in the later epochs. Source/phase-calibrator switching was performed at the recommended rates for each band and antennae configuration; pointing observations and slews to the satellite free zone were used when appropriate.

We used the default correlator setups for broadband continuum observations in all bands, including 8-bit samplers in S-band and 3-bit samplers in X and Ka bands, full polarization products, and effective bandpass sizes (after data reduction) of 2.05 (S), 4.1 (X), and about 8.1 (Ka) GHz.

We reduced the data following the standard analysis procedure using the Common Astronomy Software Applications (CASA) and the Astronomical Image Processing System (AIPS) software packages. The initial flagging and calibration used the VLA Scripted Calibration Pipeline v1.3.8 in CASA v4.6.0. For the mixed set-up observations, after initial flagging, the data for different bands were selected and separated into multiple files by running the CASA task split, before running the pipeline for each band separately. The pipeline output was inspected critically, and, in some cases, further manual flagging (and then pipeline re-calibration) was also performed. After flagging and calibration, the visibility data for the target field were converted using the CASA task exportuvfits to standard FITS format for imaging and self-calibration, which was then performed in AIPS (version 31DEC16) using the AIPS tasks IMAGR and CALIB. CLEANing boxes (including strong sources in outlier fields) were used for faster convergence in imaging and deconvolution. The “phase-only” self-calibration was done using the CLEAN components iteratively. The final images thus obtained were used to estimate the flux densities and associated statistical uncertainties using JMFIT in AIPS.

As the target source in all our observations was unresolved, in JMFIT we set DOWIDTH = -1 (which adopts the restoring CLEAN beam shape for the source) to get robust estimates of the flux densities. For the S and X band data, we also checked the rough consistency of flux densities for other nearby sources in the field. The random uncertainty output by AIPS was propagated in quadrature with an estimated systematic uncertainty of $5\%$ (S and X bands) or $10\%$ (Ka band) of the measured flux, following Weston et al. (2016), to obtain an estimated total uncertainty $\sigma_{tot}$.

The X band observation on 2016 May 1 had some unidentified weather-related problem222The problem was likely a combination of bad weather, bad ionosphere, and low inclination. Gusty winds approached the recommended X band limit of 15 m s*-1*, and towards the end of the observation 7 antennae started auto-stowing in sequence for up to 28 min due to high winds. API RMS phase reached at least 26.4$\degree$, close to the 30$\degree$ limit. The target elevation ranged from 25 to 11$\degree$, close to the recommended limit., and we were unable to estimate the flux density of the target robustly for this epoch. However, comparing the flux density of three other sources in the field with that of the adjacent epochs, finding them to be too dim by a factor $\approx$3.6, and scaling both the target flux and the observed RMS noise value by 3.6, we obtained an estimate with a large uncertainty. Reassuringly, the estimated value is intermediate between the flux densities observed on April 4 and May 24.

2.3 X-rays

Chandra observed the system in Cycle 17 for 24.76 ks on 2016 March 8.285 UT and 24.76 ks on 2016 March 9.079 UT, using chip S3 on the ACIS-S array in VFAINT mode without a grating.

We reprocessed the data following standard CIAO (v4.8, CALDB v4.7.2) procedures. Chandra detected MWC 560 as a point source always within 1 pixel (0.5 arcsec) of the expected coordinates on each dimension. We used dmstat to obtain the centroid source position (to accommodate uncertainty in the aspect solution and the expected position; for both observations, the calculated position was 7h25m51s.322,-7$\degree$44’07”.99), then extracted PSF-corrected spectra with specextract (with weight=no correctpsf=yes) using a constant circular source aperture radius of 1.93 arcsec (the 90% enclosed count fraction radius at 6.4 keV for the Chandra PSF, which is smaller at lower energies) and a circular background annulus with inner and outer radii of 18 and 30 arcsec. As a check, we also extracted the Chandra radius with a 2.5 arcsec source aperture radius for comparison. The exposure time of good data was 49.3 ks. We extracted light curves from Chandra for the same source and background regions using dmextract with 1, 5, and 12.5 ks bins.

We observed MWC 560 with the Neil Gehrels Swift Observatory for 48 ks, and reduced data from the Swift X-ray Telescope (XRT; Burrows et al. 2005) using the online product-building tool333www.swift.ac.uk/user_objects/. The source was only marginally detected at less than 3$\sigma$, so we did not attempt a centroid. We used the default grade range, and downloaded spectra corresponding to the average of all Swift data (2016 March 2 – June 1) and of the last month (2016 May 5 – June 1).

For comparison purposes, we obtained from the XMM-Newton archive an unpublished X-ray spectrum observed on 2013 April 12 (PI: Stute), reduced following standard procedures in the Science Analysis Software (SAS).

2.4 Ultraviolet

2.4.1 UV spectroscopy

We obtained 17 low-resolution (R$\sim 150$) UV Grism spectra (1700–2900 Å) on Swift UVOT (Roming et al., 2005) at roughly regular intervals between 2016 March 2 – June 1, using a total of 15 ks on Swift; the dates are listed in Appendix B.

These data were reduced using the Swift UVOT Grism UVOTPY package (version 2.1.3, Kuin 2014, an implementation of the calibration from Kuin et al., 2015). The spectra were extracted with default parameters and the locations of zeroth order contamination were determined by inspecting the images. These contaminated regions are approximations; features (especially those near contaminated regions) that did not recur between spectra obtained at different roll angles may be spurious. The contaminated regions and a few whole-spectrum quality flags are tabulated in Appendix B. The wavelength scale was shifted to match the IUE spectrum LWP19113 from 1990 November 2 (Michalitsianos et al., 1991), which we obtained from the IUE archive. We focused on the 2100-2900Å  range for comparison to older observations and to avoid auto-flagged short-wavelength data and long-wavelength second-order contamination.

We omit from our analysis three spectra noted in Appendix B as having nearby bright first order flux or as suffering from underestimated flux due to loss of lock, leaving 14 good spectra.

2.4.2 UV photometry

We observed MWC 560 with the UVOT UVM2 ($\lambda_{cen}$=2246Å, FWHM=513Å) photometric filter on Swift for 33.6 ks in event mode from 2016 Mar–Jun, using a 5x5 arcmin hardware window to minimize coincidence loss.

Reductions were performed in HEAsoft. The data were screened using a development version of uvotscreen (the stable version at the time, HEAsoft v6.21, could apply the orbit file incorrectly) using aoexpr="ELV > 10. && SAA == 0" and allowing quality flags of only 0 and 256. The photometry was extracted following the method described in Oates et al. (2009), and we checked the images at each time step. A development version of uvotsource, patched to eliminate data from debris-shadowed regions on the detector444https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/swift/docs/uvot/uvotcaldb_sss_01.pdf, was employed using a circular aperture of usually 5 arcsec radius (with larger radii during periods of mediocre tracking). We then binned the extracted data into 60 second intervals, propagating the errors appropriately. The same procedures were performed on a reference star, TYC 5396-570-1.

Our method yielded 21.6 ks of good live exposure on MWC 560. As a check, we also performed a more standard reduction using a stricter attitude filtering expression555aoexpr="ANG_DIST < 100. && ELV > 10. && SAA == 0 && SAFEHOLD==0 && SAC_ADERR<0.2 && STLOCKFL==1 && STAST_LOSSFCN<1.0e-9 && TIME{1} - TIME < 1.1 && TIME - TIME{-1} < 1.1", manual removal of additional intervals with loss of lock or mediocre tracking, and a strict 5 arcsec radius aperture. An additional few thousand seconds of observation were lost using this check method, but all long-term and short-term variability trends described in this paper were still observed.

3 Data Analysis and Results

Fig. 1 shows that at the peak of a year-long rise in optical flux tracked by Munari et al. (2016) and the AAVSO (Kafka, 2017), the 2016 outflow fast state began abruptly: Balmer absorption-line velocities suddenly doubled in January 2016 (§ 3.1), radio emissions began rising by 20$\mu$Jy/month (§ 3.2), and soft X-rays strengthened by $\sim 10$ times (§ 3.3). The fast state ended by the time that radio flux began gradually descending in July 2016 as Balmer absorption velocities decreased; still, Balmer absorption maximum velocities remained higher than 1500 km s*-1* even at the end of 2016, and as usual, portions of their profile were consistently saturated up to at least H$\gamma$. The NUV iron curtain absorption remained optically thick throughout our observations (§ 3.4). Optical/NUV flickering were also persistent throughout (§ 3.5).

3.1 Optical absorption

Sometime after the last day of 2015 and before 2016 January 21, the maximum velocity of MWC 560’s broad, blue-shifted Balmer absorption lines doubled from below 1500 km s*-1* into the 2500–3000 km s*-1* range. They stayed in this fast state until at least mid-April. The timeline is illustrated in Fig. 1. A 1-month velocity burst in 2015 January–February and the longer outflow fast state in 2016 stand out in this figure with large equivalent widths of high-velocity Balmer absorption.

H$\alpha$ velocity profiles in Fig. 2 show that 2016 January through mid-April saw consistently higher velocities than the preceding four years of quiescence, except for the 1-month burst of high-velocities in 2015 January–February. In the months following 2016 April, the Balmer velocities gradually decreased, leading to profiles (not plotted) intermediate between typical quiescent and fast state profiles. Large sections of the line profiles remained saturated in high-resolution spectra throughout, for Balmer transitions up to at least H$\gamma$ (and H$\delta$ when included in the wavelength coverage). At no point during 2016 did maximum Balmer velocities drop below 1500 km s*-1* in high-resolution spectra.

The appearance of high-velocity Balmer absorption in 2016 January was abrupt, occurring entirely during a 3-week gap in observations. The depth of the absorption line relative to the continuum decreased between 2015 Mar 8 and 2015 Oct 4, and between 2015 Oct 4 and 2015 Dec 31, which we attribute only to a rise in the underlying emission line. The absorption optical depth and velocity did not change in this period, even as late as 2015 December 31.

We supplemented our data in Fig. 1 and Fig. 2 with spectra from other sources, but only in time periods which were poorly sampled by our data (before and after the 2016 outflow fast state). 47 spectra were obtained from the Astronomical Ring for Access to Spectroscopy (ARAS) public database666http://www.astrosurf.com/aras/Aras_DataBase/Symbiotics/V694Mon.htm, comprising all spectra before or after the 2016 outflow fast state and covering H$\alpha$ with a large enough wavelength range for continuum placement. A further 5 spectra were obtained at the Special Astrophysical Observatory (SAO) by E. Barsukova and V. Goranskij (2017, private communication). We focus on H$\alpha$ because it allows the best time domain coverage in available data, by far, and a continuum that is easy to place. Individually examining the higher-order Balmer lines, of which we have poorer temporal coverage but in which the RG contributes less to the continuum, we found the same variability patterns as for H$\alpha$.

To obtain the optical absorption plots and equivalent width calculations discussed in this section, we first smoothed all spectra to a common resolving power R=450, then normalized the flux to the continuum in all spectra as measured by the median of the flux -4000 – -3000 km s*-1* blueward of the rest wavelength of H$\alpha$ (with rare exceptions on narrow-band spectra that do not include this wavelength range, in which case a continuum region was manually obtained). Then we measured the equivalent width between -1500 and -3000 km s*-1*, excluding wavelengths with flux higher than the continuum. Lower-velocity absorption below -1500 km s*-1* was ignored in order to minimize contamination by the variable broad H$\alpha$ emission line. The most variable part of the H$\alpha$ absorption line throughout the 2012–2016 period was on the high-velocity end—and the Balmer absorption lines were always contiguous with their associated emission lines—so this method accurately reflects changes in the whole line. We closely examined the spectra individually, validating that this method did not exaggerate the abruptness of line variability.

In Appendix C, we demonstrate a correlation between high-velocity Balmer absorption strength and optical/NUV flux on week time-scales throughout the 2016 outflow fast state, as both varied together in a narrow range around their maxima. High-resolution velocity profiles of H$\alpha$, H$\beta$, and Fe ii in the fast state are also presented in that appendix, along with further details on the spectral smoothing used in Fig. 1 and Fig. 2.

3.2 Radio rise

The onset of the 2016 outflow fast state coincided with a rapid, roughly 20$\mu$Jy/month rise in flat-spectrum radio emissions up to a maximum of 175 $\pm$ 10 $\mu$Jy at 9.8 GHz on 2016 July 29, about 5 times brighter than the flux density observed in the only prior radio detection of this system on 2014 October 2777A non-detection was incorrectly reported for the 2014 October 2 observation in Weston (2016) and Lucy et al. (2016), and corrected by the erratum Lucy et al. (2017).. We plot the 9.8 GHz flux density measurements as a function of time in Fig. 1. The flux densities observed at 3.1, 9.8, and 33.1 GHz are tabulated in Table 1.

The radio spectrum was always flat between all observed bands, with spectral index $\alpha$ ($F\propto\nu^{\alpha}$, F=flux and $\nu$=frequency) between -0.1 and 0.14 (modulo uncertainties $\sim 0.1$); Table 2 lists these measurements. In all cases except the SKa and XKa indices for October, the observations used to calculate $\alpha$ were obtained nearly simultaneously, within 3 hours of each other. In all cases, the in-band spectral indices are also consistent with being flat within large uncertainties.

No intra-observation variability was apparent from the visibilities; in particular, we checked for X-band variability within our brightest observation (2016 July 29) by modelling all other sources in the field, subtracting the model using uvsub in CASA, and examining the residuals in the visibilities. We also imaged that observation in four equal time sections. The flux of MWC 560 in this observation was constant within uncertainties (which were 3–4 times larger in the shorter images than for the whole observation).

We checked for polarized emission for the 2016 July 29 and October 17 observations using S and X band data, and found none. MWC 560 was not detected in Stokes Q and U, in which the noise was only slightly higher than in Stokes I.

3.3 Strengthened soft X-ray flux

Our Chandra spectrum on 2016 March 8–9 showed the emergence of a soft X-ray component about $10$ times stronger than during the 2007 and 2013 X-ray epochs. The hard X-ray component remained as weak as it was in 2013. The source aperture contained $254$ total source counts summed over the pair of exposures, amounting to $5.15\pm 0.3(1\sigma)\times 10^{-3}$ source counts s*-1*. Our data and best fit model are plotted in Fig. 3, and the best fit parameters are listed in Table 3.

The soft component was weak in 2007 and 2013 and strong in 2016, while the hard component was strong in 2007 and weak in 2013 and 2016. This result is obvious from Fig. 4, which plots the 2007 September 27 best fit model reported by Stute & Sahai (2009), our background-subtracted reduction of the 2013 April 12 XMM PN data, and our 2016 best fit model.

Following the MWC 560 models used by Stute & Sahai (2009) and the general symbiotic star models used by Luna et al. (2013), we obtained our 2016 model by fitting both the hard and soft components with separately-absorbed (photoelectric absorption with Wisconsin cross-sections), collisionally-ionized, optically-thin diffuse plasma emission spectra. These models were inspired by Patterson & Raymond (1985) for emission of a hard component by the extended halo of a WD accretion disc’s boundary layer (BL), and Muerset et al. (1997) for emission of a soft ($\beta$-type) component by colliding winds; the soft component is almost certainly not supersoft. In CIAO’s Sherpa, we grouped counts from our co-added Chandra spectra into bins of 20 source-region counts, and fit the spectrum over the range for which the Chandra energy scale is calibrated (0.277–9.886 keV) with a (wabssoft$\times$apecsoft)$+$(wabshard$\times$apechard) model. The small flux of the hard component in our spectrum required us to fix its temperature to the 11.26 keV best fit obtained during its brightest state, the 2007 epoch (Stute & Sahai, 2009).

The soft component, which has a temperature somewhere between 0.1 and 1 keV and remained roughly constant between 2007 and 2013, had an observed, absorbed flux at least 7 times brighter in 2016 (best fit 13 times) than in the prior epochs. This strengthening of the soft X-ray flux is plotted in Fig. 1. The 2016 observed flux from the soft component model was $1.7_{-0.2}^{+0.3}\times 10^{-14}$ erg s*-1* cm*-2*, which also acts as a lower limit on the intrinsic, unabsorbed flux. The upper limit on the unabsorbed flux is poorly constrained888This is because the normalization of the soft component flux is degenerate with the absorbing hydrogen column. Absorption mainly affects the peak and low-energy flank of the observed soft component flux, so at high temperatures in the allowed parameter space, the high-energy flank places an upper limit to the normalization. At lower temperatures, however, where the model’s high-energy flank is weak relative to the rest of the profile, flux normalization and absorbing column can be freely increased together without much affecting the modelled spectrum. at $\lesssim 3\times 10^{-11}$ erg s*-1* cm*-2*. These correspond to luminosities of $1.3^{+0.2}_{-0.2}\times 10^{31}$ erg s*-1* (d/2.5kpc${}^{2})$ and $\lesssim 2\times 10^{34}$ erg s*-1* (d/2.5kpc${}^{2})$, respectively.

The hard component, which by 2013 had dimmed to $0.3_{-0.3}^{+0.2}$ times its 2007 observed flux (i.e., dimmed to $\leq$1$\times 10^{-13}$ erg s*-1* cm*-2*), remained dim in 2016. No statistically significant change occurred in this component between 2013 and 2016.

In Appendix D, we demonstrate that there is no evidence for X-ray variability during the 2016 outflow fast state; however, the source was not bright enough for strong constraints. In Appendix E, we report the details of our fitting routines and the calculation of inter-epoch variability constraints. In Appendix F, we demonstrate that optical loading definitely does not affect our Chandra observations or our main conclusions, and probably does not affect our Swift XRT observations.

3.4 Persistent ultraviolet absorption

Our NUV spectra, which span 2016 March 29 to 2016 June 1, show that the “iron curtain” of overlapping absorption troughs from Fe ii and similar ions (Michalitsianos et al., 1991; Lucy et al., 2018) remained optically thick throughout the 2016 outflow fast state. These overlapping broad lines absorb on the continuum and leave behind a pseudo-continuum that looks like, but is not, an emission spectrum.

Fig. 5 shows that our 14 usable Swift UV grism spectra keep a roughly constant shape, with the differences between them best explained as a simple scaling factor applied uniformly to the whole NUV spectrum; the variability is in the underlying continuum, not in the absorption. Comparing to MWC 560’s spectral morphologies from prior epochs reported in Bond et al. (1984) and Michalitsianos et al. (1991) and also categorized in Lucy et al. (2018), the 2016 spectra are a reasonable match to the 1990 March 14 spectral morphology, and an even better match to the 1991 September 29 spectral morphology due to additional absorption between 2400-2550Å, which we ascribe to Fe*+* lines with an upper level of 2.7 eV.

We further infer that the underlying continuum being absorbed upon at the NUV flux maxima in 2016, March 2 and April 4, is very similar to the spectrum observed by Michalitsianos et al. (1991) on 1990 April 29, which is largely unabsorbed (excepting Mg ii). We retrieved this spectrum (International Ultraviolet Observer LWP17832, R$\approx$200-350 comparable to Swift) from the Mikulski Archive for Space Telescopes (MAST), and include it for reference in Fig. 5.

3.5 Persistent optical+NUV flickering

The rapid flickering in the optical and NUV that almost always characterizes the MWC 560 system (e.g., Bond et al., 1984; Tomov et al., 1996; Zamanov et al., 2011a, c) persisted throughout the 2016 outflow fast state. Densely-sampled AAVSO V-band observations covered February through April, our Swift UVM2 observations covered late-February into June, and strong variability on short time-scales was reliably observed throughout.

V-band magnitudes typically varied by at least 0.1 mag per 15–30 minutes, and the light curves strongly resembled past periods of strong flickering (comparable, e.g., to the bottom panels of Figure 1 in Zamanov et al., 2011a). For example, our Fig. 6 reproduces a characteristic light curve from the AAVSO (Kafka, 2017) on 2016 March 1–2, continuous over 11 hours thanks to volunteers observing at widely separated geographic locations. To verify this result, we confirmed that there was no systematic relationship between flux and airmass in these data.

The expected NUV counterpart to MWC 560’s optical flickering is observationally confirmed here for the first time; the most robust examples are shown in Fig. 7. Swift UVM2 ($\lambda_{cen}$=2246Å, FWHM=513Å) magnitude varied by up to 0.4 magnitudes peak-to-peak over time-scales as short as 10 minutes, consistent with flickering. Longer time-scale NUV flux variability is discussed and plotted in Appendix C.

It is virtually certain that the rapid UV variability in Fig. 7, and most of the remaining UV variability included in Appendix C, is real and not a systematic error. The reference star, TYC 5396-570-1 (AB UVM2 $\approx$15.8 magnitudes), was much less variable than MWC 560 on all time-scales; less than 2% of its light curve deviated by more than 3$\sigma$ (0.15 magnitudes) from the reference star’s median magnitude. In particular, the MWC 560 light curves in Fig. 7 were accompanied by stability in the reference star, with less than 20% of the data deviating from the local median by more than 1$\sigma$, and no noticeable evidence of trends in the direction of the MWC 560 variability. Furthermore, the UVOT instrument is known to be very stable, and has been shown by the Swift team to vary by less than 0.1 magnitudes in the brightness measured for a WD calibration source placed at one location on the detector over the course of weeks.

4 Discussion

At the peak of a year-long rise in accretion rate through the visual-emitting disc, the dense accretion disc outflow abruptly jumped in power in January 2016. The dense outflow remained fast and stable for several months (§ 4.1.1), steadily feeding a lower-density region of radio-emitting gas (§ 4.1.2), and shocks at the collision of high-velocity and low-velocity components of the dense outflow began radiating a large soft X-ray flux (§ 4.1.3). The inner accretion disc remained intact throughout (§ 4.2). Properties of the system during this outflow fast state, calculated throughout these sections, are listed in Table 4. The stability of the outflow and the disc in 2016 stands in marked contrast to the instability and disruption of both outflow and disc in 1990, even though the 2016 and 1990 brightening events reached a similar peak accretion rate. We ascribe 2016’s stability to the regulatory and stabilizing effect of the outflow on the accretion disc during both a short high-velocity burst in early-2015 and during the 2016 outflow fast state (§ 4.3).

4.1 The abrupt onset of a stable fast state

4.1.1 The dense, stable, fast absorption-line outflow

The best explanation for the jump in Balmer absorption maximum velocities between 2015 December 31 and 2016 January 21 (§ 3.1) is the sudden appearance of high-velocity material in MWC 560’s dense outflow. The appearance of high-velocity Balmer absorption occurred abruptly only towards the peak of a year-long, smooth rise in optical flux, without any evidence of a commensurately abrupt change in photoionizing flux, so it is unlikely to have been due to a photoionization effect. The new fast flow remained until at least mid-April, sometime after which it began very gradually slowing down throughout the year.

After the abrupt velocity jump, the outflow was remarkably stable and predictable throughout its 2016 fast state. High-velocity Balmer absorption, which was always present during the fast state, varied slightly on week time-scales in time with the varying optical and NUV flux; this correlation between optical depth and broadband flux was likely attributable to a photoionization effect (Appendix C).

As usual in MWC 560, the absorption-line outflow was dense. The saturation to black of large portions of the H$\gamma$ absorption trough (§ 3.1; $\tau\gtrapprox 1.6$) required at minimum a density n${}_{H}\gtrsim 10^{6.5}$ cm*-3* and a column density N${}_{H}\gtrsim 10^{23}$ cm*-2*, for any ionization parameter and assuming turbulent velocities $\gtrsim 100$ km s*-1* (Hamann et al., 2019; Williams et al., 2017). The persistence of the NUV iron curtain (§ 3.4) is also consistent with continually high densities, although it is not as restrictive: n${}_{H}\gtrsim 10^{5}$ cm*-3* for excited states up to 1.1 eV above the ground state (Lucy et al., 2014), which we certainly observe in MWC 560. We also tentatively identified Fe ii transitions from 2.7 eV above the ground state; by comparison, only some quasars with broad Balmer absorption also feature these extremely high-excitation lines, so they may impose stricter lower limits on the density in photoionization modelling (Aoki, 2010).

4.1.2 Steadily feeding a lower-density radio emission region

The radio data provide further evidence that January 2016 marked a jump in outflow power and the onset of a stable fast state. In this section, we argue that the radio emissions originated as optically-thin thermal emission in a region of low-density ($\lesssim 10^{5.5}$ cm*-3*) gas being steadily fed by the fast, denser absorption-line outflow.

The flat-spectrum radio flux density doubled in the 4 months between 2016 April 4 and July 29, growing even as the optical/NUV flux varied up and down. The slope of this radio rise is consistent with linear growth starting at its quiescent-state value in January 2016 (inferred to be 37 $\mu$Jy per the 2014 observation), roughly corresponding to the onset of fast optical absorption. Further supporting a fast-outflow origin for the radio emissions, the H$\beta$ absorption profile for 2014 October 30, obtained less than a month after the quiescent-state radio flux was measured to be weak, lacks the high-velocity absorption component that appeared in January 2016.

We favour thermal emission as the radio source, in a growing region steadily fed by the stable absorption-line outflow. The radio spectral index $\alpha$ is in the range -0.17 to 0.25 (which includes measurement uncertainties at the edges) and categorically rules out the $\alpha\geqslant 0.6$ often observed from many symbiotic star nebulae due to optically-thick thermal bremsstrahlung emission with a radially-dependent photosphere from an asymmetrically-ionized wind (Seaquist et al., 1993; Seaquist & Taylor, 1990; Seaquist et al., 1984; Reynolds, 1986; Wright & Barlow, 1975; Panagia & Felli, 1975), as well as the $\alpha\approx-0.7$ of optically thin synchrotron emission. Indeed, while symbiotics—especially S-type symbiotics like MWC 560—seem to generally have $\alpha$>0.6 in quiescence (Seaquist et al. 1993; Seaquist & Taylor 1990), flatter indices tend to emerge from jets (e.g., Weston, 2016). An outflow origin is therefore more plausible than any alternative source for the radio emission.

Furthermore, an optically thick / self-absorbed synchrotron emission mechanism, such as is favoured by Coppejans et al. (2016), Russell et al. (2016), and Körding et al. (2008) for their dwarf nova radio outbursts, is not consistent with the brightness temperatures for our data unless the emitting source is extremely small. In particular, the brightness temperature at 10 GHz (where there is no doubt the spectrum is flat) ranges from 0.6–1.3$\times 10^{7}$ (1 au/s)2 (d/2.5kpc)2 K where s is the average of the outflow’s major and minor axes, and d is the distance to the system. The kinetic temperature of synchrotron-emitting electrons at this frequency is 1.0$\times 10^{13}$ (B / 10*-4* Gauss)-1/2 K where B is the magnetic field strength (e.g., Condon & Ransom, 2016). The brightness temperature must be on the order of the kinetic temperature of the emitting electrons for synchrotron self-absorption to be significant and the spectrum to flatten (e.g., Williams, 1963). Given the speed of the outflow (travelling roughly 1 au away from the disc per day), it would be difficult to contrive a situation where the emitting region is compact or collimated sufficiently to fulfil this condition. Fulfilling the condition becomes even more difficult (by a factor of about 20) if the spectrum really is flat out to 33 GHz, as it appears to be.

Given our conclusion that the emission is thermal, the flat spectrum suggests an almost fully optically thin emitting region, which places an upper limit on the density and a very rough lower limit on the mass outflow rate. For such a gas, the Rayleigh-Jeans Planck approximation yields

[TABLE]

where $F_{\nu}$ is flux density in mJy, $\theta_{s}$ is the angular size of the emitting gas in the plane of the sky, the electron temperature $T_{e}$ is set to the gas temperature, and the optical depth is

[TABLE]

where s is the size of the emitting gas along the line of sight, for the simple case of uniform gas density equal to the electron density $n_{e}$. For $\tau\lesssim 0.1$, $T_{e}\lesssim 10^{4}$ K, $F_{\nu}=175\leavevmode\nobreak\ \mu$Jy, and $\nu=$ 10 GHz, we obtain an emitting region size scale in the plane of the sky of $\gtrsim 6\times 10^{-4}$ (d/2.5 kpc) pc (i.e., $\gtrsim 120$ au; as a consistency check, this scale could be reached by a bipolar outflow in $\gtrsim 50$ days at 2000 km s*-1*). Assuming that the size scale in the plane of the sky is comparable to the size scale along the line of sight, we obtain a density $\lesssim 10^{5.5}$ (2.5 kpc/d)1/2 cm*-3*. As a check on plausible mass outflow rates, we can divide the radio emitting region’s mass by the time-scale over which it was formed. If it took about 200 days for the radio emission to reach 175$\mu$Jy (about the time between the emergence of high-velocity optical absorption and the 2016 July 29 radio observation), then a region with this volume and density could be filled by a mass outflow rate $\gtrsim 10^{-6}$ (d/2.5 kpc)5/2 M*⊙* yr*-1*, modulo the simplifying assumptions made for this toy model. However, if there is clumping in the outflow, the mass outflow rate could be lower.

4.1.3 Soft X-rays: shocks at the collision

of fast and slow absorbers

LABEL:fignew suggests that strengthened soft X-rays during the 2016 outflow fast state may have originated in a shock at the collision of the new high-velocity absorbers and the pre-existing low-velocity absorbers. The optical spectrum most coeval to the 2007 X-ray epoch exhibited only high-velocity Balmer absorption, lacking absorption below $\sim$1000 km s*-1*. Spectra obtained during the 2012-2014 quiescent-outflow periods, including one low-resolution spectrum obtained within 2.5 months of the 2013 X-ray epoch (an excellent match to the echelle spectrum used in LABEL:fignew, when the latter is smoothed to lower resolution), exhibited only lower-velocity absorption. Soft X-rays were weak in both 2007 and 2013. It was only in the 2016 outflow fast state, when high-velocity and low-velocity absorbers coexisted, that the soft X-ray component brightened by about an order of magnitude or more.

Our fit to the soft X-ray component observed on 2016 March 8-9 constrained the diffuse, collisionally-ionized gas component temperature to between 0.1 and 1.0 keV. This corresponds to a strong-shock velocity of 300 to 900 km s*-1*, consistent with the differential velocities between the high-velocity and low-velocity absorbers. The soft X-rays emission measure, n2V$\sim 10^{54}$ (d/2.5 kpc)2 cm*-3*, is easy to achieve; assuming a density $10^{7}$ cm*-3* and a spherical toy model, this corresponds to an emission-region radius $\sim 10^{-6}$ pc (i.e., $\sim 0.1$ au), which can be traversed by a 2000 km s*-1* outflow in mere hours. In contrast, the soft X-ray component is likely too hard to be supersoft emission from nuclear burning, and too soft to be boundary-layer emission. We can certainly rule out softening of a previously observed hard X-ray component; in the 2013 epoch there was no comparably bright hard component, and even at the hard component’s maximum in 2007 it was dimmer than the 2016 soft component.

4.2 Intact inner disc

The inner accretion disc of MWC 560 remained intact throughout 2016, and the increase in luminosity was powered by accretion alone without nuclear burning. Variability in optical/NUV photometry corresponds to variability in the accretion rate through the optical/NUV-emitting parts of the disc.

Persistent optical flickering was observed throughout the 2016 outflow fast state (§ 3.5) and in every preceding monitoring project except in 1990 (e.g., Bond et al. 1984; Tomov et al. 1996; Zamanov et al. 2011c; Zamanov et al. 2011b; Zamanov et al. 2011a). Rapid flickering most likely comes from instabilities in an inner accretion disc near the boundary layer, where the instability time-scales are short. While disc-less models are presumably conceivable, to our knowledge there is no evidence for disc-less flickering from accreting WDs (see review in Sokoloski et al., 2001). Hot spots can produce some flickering where the accretion stream hits a disc, but high frequency and large amplitude flickering from accreting WDs is generally caused by the inner disc or BL (see the introduction to Bruch, 2015, for a review). We also confirmed the expected NUV flickering in MWC 560 for the first time; the NUV spectral morphology remained constant, so the variability we observed must have originated in a variable NUV continuum. Luna et al. (2013) have shown that flickering is suppressed by a flood of invariant light if there is nuclear burning on the WD surface, so MWC 560 must be powered by accretion alone. Indeed, the short-term and long-term stability of the absorption NUV and optical absorption spectra (§ 3.1,§ 3.4,§ 4.1.1) further disfavors a thermonuclear nova interpretation.

We estimate a 2016 peak (V$\approx$8.8) accretion disc bolometric luminosity of about 1800 (d/2.5 kpc)2 L*⊙, and a quiescent (V$\sim$10.5) accretion disc bolometric luminosity of about 300 $(d/2.5kpc)^{2}$ L⊙. These correspond to maximum and quiescent accretion rates of 6 $\times 10^{-7}$ and 1 $\times\leavevmode\nobreak\ 10^{-7}$ (d/2.5 kpc)2 (R / 0.01 R⊙) (0.9 M⊙/ M) M⊙* yr*-1*, respectively (assuming $\dot{M}$=LR/GM). The reddening, E(B-V)=0.15 $\pm$ 0.05 (Schmid et al. 2001, and our Appendix A), introduces a $\pm$25% uncertainty in the accretion rates and luminosities, in addition to any uncertainty in the distance or in disentangling the RG flux from the accretion disc. In Appendix A, we review and update the Meier et al. (1996) and Schmid et al. (2001) case for a distance to MWC 560 of 2.5 kpc, including new supporting evidence. In Appendix G, we report our method for calculating the accretion rates and luminosities.

The estimated 2012–2014 quiescent accretion rate is very high for a symbiotic star with flickering (i.e., without WD surface burning). Comparison to some theoretical expectations for WD surface burning as a function of WD mass and accretion rate, as in figure 1 of Wolf et al. (2013), yields interesting results. For WD masses less than 0.9 M*⊙, MWC 560’s quiescent accretion rate would be expected to lead to stable burning, inconsistent with observations—but the large hard X-ray component temperature observed in MWC 560 by Stute & Sahai (2009) does imply a high white dwarf mass, so this is not a problem. At 0.9 M⊙, stable burning is narrowly avoided, but a recurrent nova is expected every century or so; in other words, MWC 560 may be due for an imminent nova. For WD masses higher than 0.9 M⊙, the nova recurrence time quickly becomes too short, inconsistent MWC 560’s 9-decade-long optical light curve (Leibowitz & Formiggini, 2015). If the WD mass is higher than 0.9 M⊙* (as perhaps suggested by the comparable boundary layer cooling flow temperatures between MWC 560 and RT Cru / T CrB; Stute & Sahai 2009), or a nova does not occur sometime soon, that may support competing theoretical pictures (e.g., Starrfield et al. 2012a, b) or suggest an important role for outflows in preventing novae. However, some caution is warranted due to the uncertainties introduced by flux modelling, reddening, and especially distance. A more careful assessment of MWC 560’s possibly liminal location in WD mass / accretion rate parameter space may be possible after future GAIA data releases with binary solutions and smaller uncertainties.

4.3 Self-regulation of the disc by its outflow

Differences between brightening events in 1990 and 2016 suggest that the 1-month outflow burst in January–February 2015 and the sustained outflow fast state in January–July 2016 (Fig. 1) may have helped stabilize the accretion system and helped keep the inner disc intact during an increase in accretion rate through the disc. Dramatic disruptions in the outflow and the disc in 1990 did not repeat in 2016, even though both events reached similar peak accretion rates.

4.3.1 Historical overview

In the 9 decades of observations preceding 2016 (Luthardt, 1991; Leibowitz & Formiggini, 2015; Munari et al., 2016), two other optical brightening events came close to the maximum optical luminosity of 2016. The 2010-2011 event reached at least V$\approx$9.6, and was accompanied by absorption velocities as high as -4200 km s*-1* (Goranskij et al., 2011). The 1989–1990 event reached at least V$\approx$9.2, and was accompanied by rapidly variable broad absorption lines with velocities as fast as -6000 km s*-1* (Tomov et al., 1990, 1992). 1989–1990 may have included a step-function-like shift in the light curve, separate from the postulated orbital and RG periodicities that may explain the timing of the three events (Leibowitz & Formiggini, 2015). At least, 1990 is consistent with a secular brightening over the last century (Munari et al., 2016); the three optically-brightest events have all occurred during or after 1990, and the average optical brightness after 1990 has been about triple the average optical brightness before 1990 (Leibowitz & Formiggini, 2015).

4.3.2 1990 vs. 2016

In one respect, 1990 and 2016 were similar: using the same distance and reddening, and a similar flux model, Schmid et al. (2001) obtained about the same peak accretion rate and luminosity for the 1990 event as we did for 2016 (§ 4.2). But otherwise, the 1990 and 2016 events differed starkly.

(1) In 1990, the Balmer absorbers were rapidly variable, appearing and vanishing in matters of days, detached from their associated emission lines, and ejected at velocities up to -6000 km s*-1* (Tomov et al., 1990, 1992); then, later in 1990, the outflow slowed to a very unusually small $\lesssim 900$ km s*-1* until about a year after the outburst (Zamanov et al., 2011a). In 2016, the outflow was much more stable (§ 3.1: Balmer lines varied only slightly (slowly and predictably in time with the optical/NUV flux, likely a photoionization effect; Appendix C), never vanished, stayed contiguous with their emission, and never went faster than $\approx 3000$km s*-1*; then, later in 2016, the Balmer lines slowed, but never fell below $\approx 1500$km s*-1*.

(2) Over the course of the 1990 outburst, the veil of Fe ii absorption in the UV lifted, clearing entirely by the end of April 1990 (Lucy et al., 2018; Maran et al., 1991); as discussed in Lucy et al. (2018), this likely represented a temporary switch from persistent outflow to discrete mass ejection, resembling in the latter state the P Cygni phase of some novae. In 2016, this iron curtain remained optically thick throughout and the spectral morphology did not vary significantly (§ 3.4), consistent with the stable and persistent outflow that we saw in the optical.

(3) Fast optical flickering was partially suppressed during the 1990 outburst, and thoroughly suppressed by the next observing season later in that year (Zamanov et al., 2011a). U band fast photometry collected during the 1990 outburst peak and discussed in Zamanov et al. (2011a) show comparatively smooth light curves. In 6 hours of U band photometry by T. Tomov over 4 days during the 1990 high state, provided to us by R. Zamanov (2017, private communication), there were slow 0.1 mag oscillations on hour time-scales, and on shorter time-scales one incident of 0.2 mag flickering and a few incidents of 0.1 flickering. Afterwards, at the same time as the Balmer velocities were unusually slow, essentially no flickering at all was observed on time-scales of minutes (Zamanov et al., 2011a). In contrast, flickering persisted throughout 2016 (§ 3.5); there was no evidence at all for suppression (though poor sampling after mid-2016 warrants a little caution).

In brief, it appears that in 1990, MWC 560 underwent a dramatic disruption of the flicker-producing inner disc (Zamanov et al., 2011a), perhaps as a result of ejections of mass evacuating the inner disc, and that this disruption interrupted the launching mechanism for the high-velocity outflow until the inner disc was rebuilt a year later. While no dramatic colour changes were observed in 1990, reprocessing of the accretion disc light by an optically thick, low velocity wind, like that proposed for MWC 560 by Panferov et al. (1997), could make any changes in accretion disc size and temperature profile difficult to observe in the UV and optical. Disc evacuation has also been proposed to have occurred in the symbiotic star CH Cyg, which is believed to have once been evacuated by a jet, causing its flickering to temporarily cease (Sokoloski & Kenyon, 2003a, b). Such incidents also have precedent in X-ray binaries, as recently demonstrated by V404 Cyg (Muñoz-Darias et al., 2016). Neilsen et al. (2011) also suggested that accretion discs in X-ray binary systems like GRS 1915+105 may suppress themselves by carrying away the mass of the inner-disc that launches them.

In 2016, none of that happened. The outflow power increased abruptly as the system reached its peak accretion rate through the visual-emitting disc, but the new outflow was stable, the inner disc was intact, and neither the outflow nor the inner disc were destroyed by the event. Fig. 8 shows a schematic diagram of MWC 560 before, during, and after the 1990 and 2016 accretion rate peaks, illustrating the difference between these two events.

4.3.3 The outflow as a gatekeeper

Two interrelated processes may explain the stability of the disc and the outflow in 2016: the timely evacuation in the outflow of mass from the accretion disc, and a long-term trend towards higher average accretion rates over the course of the last century. We propose that the outflow helped maintain equilibrium in the disc, slowing down changes in the disc in 2015, finally halting changes in the disc in 2016, and thereby preventing the inner disc’s evacuation.

There was a month-long expulsion of high-velocity material in January–February 2015 (Fig. 1), just before the optical flux began to slowly rise. Perhaps an early-2015 disc instability or an increase in the RG mass transfer rate, which would not necessarily immediately be detectable in the optical depending on where in the disc it manifested, was slowed and stabilized by a wind burst carrying away excess mass and angular momentum. Although optical rises reliably predict increased outflow velocity throughout the history of MWC 560 (Tomov et al., 1990; Iijima, 2001; Goranskij et al., 2011) in a clear cause-and-effect relationship, such outflow velocity jumps have previously occurred in optical flux quiescence too (Iijima, 2001); the system clearly varies in ways that we sometimes cannot detect, and the outflow could play a role in suppressing the propagation of these changes through the disc.

The outflow might suppress any heating or cooling wave that reaches its launching radius in the disc. After February 2015, the accretion rate through the visual-emitting disc continued to rise (Fig. 1), but the outflow remained slow and static until January 2016, when it jumped in power and subsequently persisted steadily for months as the optical light curve flattened (§ 4.1). We speculate that the delay between the start of the optical rise and the initiation of a persistent fast outflow could be due to feedback from the outflow burst and the multiple directions in which heating and cooling waves can travel through an accretion disc. The propagation of heating and cooling waves across the disc in a dwarf-nova-type outburst can involve mass being transported along the angular momentum gradient of the disc in both directions (e.g., Lasota, 2001; Schreiber et al., 2003). It may be that the fast outflow launching radius was roughly the last radius to experience the outburst responsible for the 2016 accretion rate peak. For example, an outside-in outburst might only affect a centrally-launched wind after the increased accretion rate finally propagated to the centre of the disc, at which point a fast outflow could launch while the disc stabilized. The 2016 optical rise occurred near the confluence of two long-term periodicities in MWC 560’s light curve (Leibowitz & Formiggini, 2015). A rise in the external accretion rate (i.e., rate of mass flow into the disc) as periastron approached during a period of high RG mass loss could induce multiple instabilities in the disc; those at the outflow launching radius would be suppressed by an outflow burst (and by the difficulty of propagating to higher surface densities and angular momenta; e.g., Schreiber et al., 2003), while others would have time to propagate through the disc before restarting the outflow, themselves getting suppressed by the 2016 outflow fast state.

What physical differences between 2016 and 1990 allowed the formation of a stable fast outflow in 2016, instead of driving catastrophic discrete evacuations of the inner disc as in 1990? The optical light curve in 2014–2016 was about the same as in 1988–1990, but the 2006–2013 period was typically at least twice as bright as the 1980–1987 period (Munari et al., 2016; Doroshenko et al., 1993; Luthardt, 1991), consistent with a secular brightening throughout the last century (Munari et al., 2016) or with a step-function tripling of the average optical luminosity in 1989–1990 (Leibowitz & Formiggini, 2015). So it may be that the structure of the disc (and therefore its outflow-driving mechanism) was acclimatized, on a time-scale of at least several years, to a much higher accretion rate through the disc in 2016 than in 1990. The large discs of symbiotic stars plausibly allow for viscous time-scales of at least several years to govern the nature of the as-yet-undetermined outflow-driving mechanism. Alternatively, the difference between 2016 and 1990 could involve changes in the boundary layer that are not always traced by the optical flux; the disappearance of hard X-rays between 2007 and 2013, and their continued absence in our 2016 observations, might indicate that the optical depth of the boundary layer had been increasing over at least the last decade. Future observations could test whether an optically thick boundary layer is necessary to initiate a stable outflow without evacuating the inner disc in discrete mass ejections.

Neilsen et al. (2011) suggests that, in the X-ray binary microquasar GRS 1915+105, a high-mass (Ṁwind $\geq$ 15Ṁacc) wind may be “effectively acting as a gatekeeper or a valve for the external accretion rate, and facilitating or inhibiting state transitions” (in that case, via Shields et al. 1986 oscillations). Our 2016 observations provide strong new evidence that the role of MWC 560’s outflow is similar to the role of GRS 1915+105’s outflow, a comparison first made in broader terms by Zamanov et al. (2011a). During the 2016 outflow fast state of MWC 560, the mass outflow rate estimated from our radio observations (perhaps $\gtrsim 10^{-6}$ (d/2.5 kpc)5/2 M*⊙* yr*-1* if the outflow has uniform density; § 4.1.2) may have been commensurate with the accretion rate through the disc (6 $\times$ 10*-7* (d/2.5 kpc)2 (R / 0.01 R*⊙) (0.9 M⊙/ M) M⊙* yr*-1*; § 4.2, Appendix G), although neither is sufficiently well constrained for this comparison to be more than a consistency check. The possible evacuation of the inner disc of MWC 560 during the 1990 outburst (Zamanov et al., 2011a), accompanied by the phenomenology of discrete mass ejection in 1990 (Lucy et al., 2018), suggests that MWC 560’s outflow may indeed be capable of acting as a valve, facilitating an accretion disc state transition in 1990 and inhibiting it in 2016.

5 Conclusions

If accretion disc outflows are common in symbiotic stars (e.g., Lucy et al., 2018; Sokoloski, 2003; Muerset et al., 1997), then the effect of outflows on their accretion discs could be fundamental to understanding their physics and evolution. The broad absorption line symbiotic star MWC 560 = V694 Mon is an important laboratory for this complex relationship. We conducted coordinated radio, optical, NUV, and X-ray spectroscopic and photometric observations of MWC 560 during an optical flux maximum in 2016, and we detected an outflow fast state from January 2016 into likely the end of July 2016. The variability in each band is plotted in Fig. 1 for the years 2012 through 2016.

1. (1)

The maximum velocity of a dense, Balmer line-absorbing outflow from MWC 560’s accretion disc abruptly doubled over the course of at most 3 weeks in January 2016, near the 2016 February 7 accretion rate peak and at the completion of a year-long rise in the accretion rate through the visual-emitting disc. The abrupt change in Balmer absorption velocity was almost certainly not due to photoionization effects, although a subsequent correlation on week time-scales between high-velocity Balmer opacity and optical/NUV flux can probably be attributed to photoionization. High velocities were stable and sustained through at least mid-April of the same year, before beginning a slow decline; even through to the end of 2016, velocities in high-resolution spectra never dropped below 1500 km s${{}^{-}1}$. The density of the Balmer-absorbing gas was $\gtrsim 10^{6.5}$ cm*-3*. (§ 3.1, § 4.1.1) 2. (2)

Radio emissions confirm an increase in outflow power at the onset of higher Balmer absorption velocities. Flat-spectrum radio emissions detected at 3.1, 9.8, and 33.1 GHz began to rise linearly at a rate of about 20 $\mu$Jy/month, even as the optical/NUV flux varied up and down. The slope of this radio rise suggests that it started at its 37 $\mu$Jy quiescent value (observed in 2014) in January 2016, around the same time as the high-velocity Balmer absorption appeared. Radio emissions reached a maximum of 175 $\pm$ 10 $\mu$Jy at 9.8 GHz on 2016 July 29 before beginning a slower decline. The emission mechanism was thermal and optically-thin, originating in gas with a density $\lesssim 10^{5.5}$ cm*-3*. We propose that this lower-density region was steadily fed by the denser Balmer absorption-line fast outflow. (§ 3.2, § 4.1.2) 3. (3)

Also during the 2016 outflow fast state, soft X-rays were observed to be brighter by an order of magnitude relative to both prior X-ray epochs (2013 and 2007). The plasma temperature of this component, constrained to between kT $\approx$ 0.1 and 1 keV, was consistent with strong-shock velocities of 300-900 km s*-1*, in turn consistent with differential velocities in the optical Balmer absorption lines. Meanwhile, the hard X-ray component in both 2016 and 2013 was about a third as bright as it had been in 2007. (§ 3.3, § 4.1.3, Fig. 4) 4. (4)

Balmer velocity profiles observed close to each of the three X-ray epochs suggest that soft X-rays are weak when only high-velocity absorption is present (2007) and when only low-velocity absorption is present (2013); only in 2016, when both high-velocity and low-velocity Balmer absorption was observed, was the soft X-ray flux strong. We propose that the soft X-ray component originates in a shock where these new fast absorbers and pre-existing slow absorbers in the absorption-line outflow collide. (§ 4.1.3, LABEL:fignew). 5. (5)

The iron curtain of overlapping absorption in the NUV remained optically thick throughout the 2016 outflow fast state, without substantial variability in its spectral morphology, further supporting the presence of a stable and sustained outflow. (§ 3.4, § 4.1.1) 6. (6)

Optical/NUV flickering shows that the inner accretion disc remained intact in 2016; the slow optical brightening from 2015 into 2016, which led to the outflow power jump, was due to an increase in the rate of accretion through the disc. Flickering by at least a typical 0.1 mag per 15–30 minutes in V band, up to 0.4 magnitudes over 10 minutes in the NUV, persisted throughout the 2016 outflow fast state—similar to that observed in quiescence and high states from the end of 1991 through 2015. This flickering requires that MWC 560’s luminosity be powered by accretion alone without WD surface burning. The high accretion rate, even in the relative quiescence of 2012–2014, may suggest that a nova should be expected in MWC 560 within the next century. (§ 3.5, § 4.2) 7. (7)

The peak accretion rate in 2016 was about the same as the peak accretion rate in the 1990 outburst of the same system. Despite this, the two events differed dramatically: In 1990, rapidly variable mass ejections up to 6000 km s*-1* appeared to evacuate the inner accretion disc, leading to the cessation of flickering and the suppression of the outflow to velocities below 900 km s*-1* for up to a year, and the temporary disappearance of iron curtain absorption in the UV. Throughout 2016, the outflow was stable and the inner disc remained intact; strong and rapid flickering continued, Balmer velocities reached up to 3000 km s*-1* and slightly varied in time with the optical/NUV flux, and the iron curtain remained optically thick. (§ 4.3) 8. (8)

We propose that the outflow sometimes inhibits structural changes in the accretion disc and sometimes facilitates them. In 1989–1990, the outflow evacuated the inner disc in a phase of discrete mass ejection. But in 2015–2016, a 1-month outflow velocity burst in January–February 2015 and the longer, stable 2016 outflow fast state may have prevented a catastrophic evacuation of the inner disc, by carrying away excess accreting mass and suppressing heating and cooling waves as they reached the outflow’s launching radius. The difference between the disc/outflow relationship in 2016 versus 1990, which reached similar peak accretion rates, might be due to a secular increase in the decade-averaged accretion rate throughout the last century. The complex, self-regulatory relationship between a symbiotic star accretion disc and its outflow resembles X-ray binary behaviour. (§ 4.3, Fig. 8)

Acknowledgements

With thanks to Neil Gehrels.

We thank Fred Hamann for helping us understand the relationship between Balmer absorption and photoionization, Elena Barsukova and Vitaly Goranskij for generously sending additional spectra, and Josh Peek and Kirill Tchernyshyov for assistance with their ISM velocity map. We acknowledge Steve Shore for wisely encouraging ARAS observations in early 2015.

We thank volunteer observers Teofilo Arranz, Gary Walker, and Geoffrey Stone, whose data were submitted to the AAVSO and are reproduced in Fig. 6. We are similarly indebted to all the other AAVSO, ARAS, and independent observers across the world whose work made this paper possible: H. Adler, S. Aguirre, T. Atwood, D. Barrett, P. Berardi, B. Billiaert, D. Blane, E. Blown, J. Bortle, D. Boyd, J. Briol, C. Buil, F. Campos, A. Capetillo Blanco, R. Carstens, J. Castellani, W. Clark, T. Colombo, A. Debackere, X. Domingo Martinez, S. Dvorak, J. Edlin, J. Foster, R. Fournier, L. Franco, J. Garlitz, A. Garofide, A. Glez-Herrera, K. Graham, C. Gualdoni, J. Guarro Flo, J. Foster, F. Guenther, C. Hadhazi, F. Hambsch, B. Harris, D. Jakubek, P. Lake, T. Lester, S. Lowther, C. Maloney, H. Matsuyama, K. Menzies, V. Mihai, J. Montier, G. Murawski, E. Muyllaert, G. Myers, P. Nelson, M. Nicholas, O. Nickel, S. O’Connor, J. O’Neill, W. Parentals, A. Plummer, G. Poyner, S. Richard, J. Ripero Osorio, J. Ritzel, D. Rodriguez Perez, R. Sabo, L. Shotter, P. Steffey, W. Strickland, D. Suessman, F. Teyssier, T. Vale, P. Vedrenne, A. Wargin, and A. Wilson. Many of these observers’ data were reproduced in Fig. 1. We apologize to any observers we neglected to acknowledge, with deep gratitude. We thank the staff of the AAVSO International Database and the ARAS Spectral Data Base, including E. Waagen and F. Teyssier.

We thank the Chandra team (PI: B. Wilkes), the Swift team (PI: N. Gehrels), and the Very Large Array team for the discretionary time, and for their technical help throughout the observations on which this paper is based. ABL thanks K. Mukai, T. Nelson, T. Iijima, the Chandra help desk, the Swift help desk, and the GAIA help desk for useful conversations.

ABL is supported by the NSF GRFP under grant DGE-1644869. ABL and JLS are supported by Chandra award DD6-17080X. JLS is supported by NSF AST-1616646. UM is partially supported by PRIN INAF 2017 (Towards the SKA and CTA era: discovery, localisation and physics of transient sources, P.I. M. Giroletti). NR acknowledges support from the Infosys Foundation through the Infosys Young Investigator grant. GJML is a member of the CIC-CONICET (Argentina) and acknowledges support from grant #D4598, ANPCYT-PICT 0478/2014 and 0901/2017. MJD acknowledges financial support from the UK Science and Technology Facilities Council. This work was supported in part by the UK Space Agency. ABL thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant #1829740, the Brinson Foundation, and the Moore Foundation; their participation in the program has benefited this work.

We gratefully acknowledge our use of pysynphot (STScI Development Team, 2013), PyAstronomy (https://github.com/sczesla/PyAstronomy), extinction (Barbary, 2016), the Mikulski Archive for Space Telescopes (MAST), the International Ultraviolet Explorer (IUE), the Sloan Digital Sky Survey (SDSS), NASA’s Astrophysics Data System (ADS), the SIMBAD database operated at CDS, and the National Institute of Standards and Technology (NIST).

Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts. The Liverpool Telescope is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias with financial support from the UK STFC. This work was based in part on data obtained with the Asiago 1.82m Copernico (INAF Padova) and Asiago 1.22m Galileo (University of Padova) telescopes.

Appendix A Distance, the burning question

In this section, we review and update the case for a distance of about 2.5 kpc (Schmid et al., 2001; Meier et al., 1996), which yields high accretion rate estimates that may be in some tension with some empirical and theoretical expectations for typical non-burning symbiotics, or suggest that a thermonuclear nova is imminent (§ 4.2). Each individual argument relies on debatable assumptions and the uncertainty is poorly constrained, but d$\approx$2.5 kpc is supported by several largely-independent lines of reasoning.

Geometric parallax: Adopting the Bailer-Jones et al. (2018) Bayesian analysis of GAIA DR2, the parallax distance is 2.4${}^{+1.3}_{-0.7}$ kpc. Here we assume that the parallax is not affected by binary motion, and that the Bayesian prior employed works for red giants in general and symbiotic stars in particular.

Reddening and dust maps: Three pieces of evidence point to an E(B-V) between 0.1 and 0.2: fits of field giant templates to the optical spectrum, the 2200Å absorption bump (Schmid et al., 2001), and the strength of interstellar absorption lines. For our part, we performed a fit of 5 Gaussians to the Na iD complex (one Gaussian for each of 4 resolved absorption lines, and one for what we suspect is remnant sky emission) in the publicly available pipeline reduction of the FEROS 07881 1990 November 14 spectrum, following the method for resolved lines described by Munari & Zwitter (1997). We obtained a total reddening of E(B-V)=0.1, which we regard as a lower limit due to uncertainties in the continuum placement and underlying emission. Using the Bayestar2017 3D map of Galactic extinction by Green et al. (2018) and incorporating their reported statistical uncertainties, E(B-V)$\geq$0.1 corresponds to a distance $\geq$2.1 kpc. The Stilism 3D dust map v. 4.1 is even more consistent with our conclusions, yielding E(B-V)=0.15 at 2.5 kpc (Lallement et al., 2018; Capitanio et al., 2017). However, the residuals between all types of 3D extinction maps and individual stars are often significant (e.g., see figure 9 in Green et al. 2015, and severe angular smoothing in Lallement et al. 2018).

Optical flux: For E(B-V)=0.15, the best-fitting RG in the optical (Appendix G) has a V band apparent magnitude of 11.8, corresponding to a 2.5 kpc distance. Here we assume that we have correctly disentangled the giant’s optical flux from the accretion disc flux, and that the giant has an absolute magnitude Mv = -0.4 – -0.5 consistent with both bulge and field giants (Allen, 1983; Whitelock & Munari, 1992).

Infared flux: MWC 560 has an average J band magnitude of 6.5, and Meier et al. (1996) used this to obtain a distance of 2.5 kpc. For a field giant with (V-J)0=5.4 (Koornneef, 1983) and E(B-V)=0.15, corresponding to Aj=0.87$\times$E(B-V)=0.13 (Savage & Mathis, 1979), we obtain a distance of 2.9 kpc. Similar results are obtained in the K band. If we assume that the accretion disc and the giant have the same reddening, the observed (V-J)0 of MWC 560 is more consistent with a field giant than a bulge giant. Still, if we did instead assume a bulge giant with higher reddening than the accretion disc, we would obtain a distance of about 1 kpc, illustrating the large model uncertainties in flux-based distance estimation.

Interstellar absorption line velocities: Schmid et al. (2001) argued for a large distance based on the radial velocities of apparently interstellar Na iD and Ca ii absorption. In the LSR frame, the resolved lines they observed actually correspond to -7, 0, +27, +35, and +52 km s*-1*. LSR radial velocities in the Galactic disc are typically around +10 km s*-1* above a flat rotation curve (V*⊙=220 km/s, d⊙=8.5 kpc) at this longitude, yielding +23 and +39 km s-1* at 1 and 2.5 kpc, respectively (Tchernyshyov & Peek, 2017, and private communication). This indeed supports a large distance $\gtrsim$2.5 kpc, although the affects of turbulence in the ISM are poorly constrained (e.g., Heiles & Troland, 2003).

Appendix B Optical and UV observation metadata

(appendices continue on following pages)

Appendix C Absorption correlated with flux: a photoionization effect?

During the 2016 outflow fast state, a clear correlation was observed between optical/NUV flux and the equivelant width of high-velocity ($|$v$|$>1500 km s*-1*) Balmer absorption. This correlation is demonstrated in Fig. 9. The NUV flux varied with the optical flux, sometimes with a larger amplitude. The strength of high-velocity absorption tracked both the optical and NUV flux.

The correlation between high-velocity Balmer absorption strength and optical/NUV flux is consistent with a photoionization effect. Balmer absorption from n=2 of neutral hydrogen occurs in the partially-ionized zone, just outside the fully-ionized H ii zone that surrounds the source of ionizing flux. As discussed in Hamann et al. (2019) and Williams et al. (2017), if the ionizing flux is not high enough to ionize the entire gas slab, then increasing the ionizing flux expands the partially-ionized zone and leads to more Balmer absorption (assuming that the column density is high enough and that the high density requirements to excite to n=2 are met). If the ionizing flux is further increased and the partially-ionized zone reaches the edge of the gas slab, the relationship reverses, and increasing the ionizing flux begins to reduce Balmer absorption and fully-ionize the whole slab. If the ionization parameter in the MWC 560 outflow is relatively low such that the former case applies, then we would expect a correlation between ionizing flux and Balmer absorption strength—which we observe. Assuming that the photoionizing flux was proportional to the optical flux, the amplitude of the optical depth variability is consistent with that expected in much of the parameter space modelled by Hamann et al. (2019) and Williams et al. (2017).

Alternatively, the correlation could be the result of real variability in the amount of high-velocity mass. If true, this would suggest that the 2016 outflow fast state was a period during which the accretion rate and the outflow power were correlated. But photoionization effects mean we must be cautious in treating the Balmer absorption line as if it were precisely representative of the outflowing mass itself in the absence of other evidence. Fortunately, as discussed in § 4.1.1 and partially thanks to radio and X-ray observations, this is not a significant problem for conclusions reached in the main text of this paper.

In order to detect the correlation between high-velocity Balmer absorption and optical/NUV flux, we smoothed the spectra to a common resolution of R=450. This was necessary due to the wide variety of spectral resolutions employed, the ambiguities of continuum placement, and the blending of emission features into the absorption line. The smoothing was performed using the instrBroadGaussFast function in the Python AstroLib (pyasl in PyAstronomy), after first re-binning with flux conserved to a uniform wavelength grid (using999http://www.astrobetter.com/blog/2013/08/12/python-tip-re-sampling-spectra-with-pysynphot/ (code by J. Lu) the spectrum and observation modules in pysynphot) where necessary. The smoothing decreases the absolute value of all the measured equivalent widths by an amount that depends only on the original resolution (i.e., by the same amount for spectra with the same original resolution), up to a few angstroms.

Finally, we present our high-resolution echelle velocity spectra during the 2016 outflow fast state for the main unblended lines for which we have high-SNR coverage: H$\alpha$ (Fig. 10), and H$\beta$ and Fe ii (Fig. 11). The Balmer lines show the same pattern of high-velocity absorption variability seen at low resolution. Fe ii exhibits a similar trend with more scatter.

Appendix D No detectable intra-epoch variability in X-rays

MWC 560 was marginally detected in Swift XRT with less than 3$\sigma$ significance, and we used this to check for variability on time-scales of months; none was detected. We obtained 3$\sigma$ upper limits with the conf task, after binning to at least 20 source-region counts and subtracting background, to the average Swift spectrum (from 2016 March 2 – June 1). Freezing all parameters to the Chandra best fit except the normalizations, we obtained upper limits $\approx$2 times the Chandra best fit values, meaning that the average X-ray flux from both components from March 2 through June 1 (sampled at roughly even intervals) was at most double the March 8-9 Chandra best fit. There is no evidence for variability; the Chandra fit was consistent with both the average Swift spectrum and the last month of Swift data (May 5 - June 1).

There was no evidence for X-ray variability within the Chandra observations either, although with only about 5 source counts per kilosecond, the system again was not bright enough for strong constraints. The model fit to the merged spectrum was consistent with both 25 ks March 8 and 9 spectra individually; following the same procedure used to constrain Swift variability, we found that the soft component normalization varied by less than $\pm$50% from the merged best fit and the hard component normalization varied by less than +100% and -75% from the merged best fit, between March 8 and 9 at 3$\sigma$. Likewise, visual inspection of full-energy-band Chandra light curves obtained with dmextract for a variety of time bins did not reveal any instances of variability in excess of the expected statistical noise. For that inspection, we estimated 1$\sigma$ errors using the Gehrels (1986) Equation 7 upper margin with S=1, a standard option in CIAO; the average statistical uncertainty for a single bin was 80%, 25%, and 15% of the source count rate for 1, 5, and 12.5 ks bins, respectively. Finally, for each binning, we calculated the observed standard deviation of the whole light curve and found it to be consistently 45-60% of these statistical uncertainties.

Appendix E X-ray fitting routine details

We used chi2gehrels statistics (low-count Poisson statistics from Gehrels, 1986) and the moncar method on the background-subtracted spectrum. The moncar method (a Monte Carlo routine that draws fit statistics for one variable at a time from parameter space, while allowing all other variables to float to their best fit at that point in parameter space) was motivated by a non-parabolic distribution of fit solutions in parameter space. To accommodate the extended, mildly bi-modal spread of allowable solutions, we increased the calculated confidence intervals to 95$\%$ ($\sim$2$\sigma$) confidence intervals, bracketing the solution spread by requiring a $\chi^{2}$ at least 3.84 above the minimum. The reduced $\chi^{2}_{\nu}$ is a little high (1.6 at the best fit), but visual inspection of Fig. 3 suggests that the fit suits the data and further model components are not motivated.

We checked our analysis against another common data reduction method, obtaining consistent results101010The best fits remained identical to within 1% for the soft component and 15% for the hard component, well within the uncertainties. The only more substantial change was a 100% and 25% stretching of the upper bound margins for the normalizations and for the hard component absorption, respectively, which does not impact our conclusions. when using a spectrum extracted from a larger 2.5 arcsec radius region and fitting from 0.3 to 10.0 keV.

The results for variability in observed flux between epochs, included in the main text, were obtained as follows. We froze all parameters to the 2016 best fit except for the normalization of the soft component, which was allowed to vary freely. When this frozen 2016 best fit was applied to the 2013 spectrum, the soft component normalization had a $1\sigma$ minimum of 0 and maximum of $2.5\times 10^{-6}$ (best fit $1.5\times 10^{-6}$). This range includes results from binning to both 20- and 40- source-region counts per bin (the latter avoids bins with negative counts in the background-subtracted spectrum, and gives our quoted best fit value). Trying to fit the background introduced too many free parameters to be useful, but we did try a version of low-count Cash statistics (wstat in Sherpa) on the un-binned background-subtracted spectrum and obtained a much smaller maximum. Similarly, when the frozen 2016 best fit was applied to the 2016 spectrum with Gehrels $\chi^{2}$ statistics, the soft component normalization was $1.9_{-0.2}^{+0.3}\times 10^{-5}$, almost independent of binning. In this way, we obtained that between the 2013 and 2016 X-ray epochs, the observed soft component flux increased by a factor of 13 (1$\sigma$ lower bound of 7).

For the inter-epoch variability in the hard component, an essentially identical method was employed: freezing the soft component to the 2013 fit obtained in the last section, freezing all hard component parameters except the normalization to the 2007 Stute & Sahai (2009) fit, and additionally fixing the weak Gaussian emission line component included in the Stute & Sahai (2009) model to have a normalization proportional to the hard component normalization.

Throughout the main text, observed X-ray flux and its derived quantities are quoted with rough 1$\sigma$ uncertainties from the fits with only 1 free parameter, while X-ray temperature and absorption columns are quoted with rough 2$\sigma$ uncertainties from the fits with 5 free parameters.

Appendix F Optical loading constraints

X-ray telescopes are susceptible to optical loading, in which optical photons reach the detector and masquerade as X-ray photons; red/infrared leaks can let optical light through the shielding. Symbiotic stars are particularly dangerous; not only do they contain bright cool giants with high red/infrared flux, they also contain colliding winds that produce soft X-rays with about the same energies as optical loading typically produces. Fortunately, the leaks and detectors on Swift and Chandra are well-enough characterized that we can check whether optical loading is an issue for any given object; here, we do so for MWC 560.

We ruled out optical loading in our Chandra observations. The Chandra optical leak is between the I and J bands. On the S3 chip, for sources dimmer than 5 magnitudes in those bands, no more than $\sim 1$ contaminant ADU appears on the brightest pixel per frame, much less than the $\approx 20$ ADU pixel*-1* required to be interpreted as an event count and much less than the 13 ADU pixel*-1* split threshold used for event grading; also, the 0.004 keV energy of an ADU is not enough to change the energy of real events significantly111111http://cxc.harvard.edu/cal/Hrma/UvIrPSF.html. The J band magnitude of MWC 560 is about 6.5 (Meier et al., 1996), and I band magnitudes around 8 were recorded by the AAVSO during the 2016 outflow fast state (Kafka, 2017).

It is interesting to note that optical loading on Chandra could be a very important issue for those symbiotic stars with I and J magnitudes brighter than -0.4, especially any observations taken on the S3 chip with a standard frame length early on in the life of the Chandra mission. As of 2006, that magnitude corresponded to the split threshold. Over time, Chandra’s filter has developed a thickening layer of contaminants that has helped plug the red leak; there is some evidence that it now takes a much brighter source to reach the split threshold, though these new tests have not yet been formally characterized (Chandra help desk 2017, private communication).

We estimated the optical loading in Swift XRT using the online tool provided by Swift121212http://www.swift.ac.uk/analysis/xrt/optical_tool.php, which required us to choose a Teff and Mbol-Mv for the source, and we adopted Teff=3330 K and Mbol-Mv=-2.15. A PC-mode XRT count threshold of 10*-4* counts s*-1* yield a V band threshold of 9.7 magnitudes, so the giant is dim enough not to cause optical loading on its own. The outflow fast state accretion disc is typically much bluer, and probably not an issue (the same threshold for a B0 star is V=8.6). However, the fast state accretion disc colour is not especially well characterized, particularly longer than $8000$Å. So we cannot totally rule out optical loading in the Swift XRT, especially during the brightest phases of the 2016 high state. Fortunately, our uncertainty on this point is not of great concern, because we only used the Swift XRT to check for variability (of which we found none).

Appendix G Estimating the accretion luminosity

To estimate the accretion disc luminosity at the 2016 peak and during quiescence, we took our spectra from 2016 April 3 (the brightest spectrum) and 2015 March 10, and dereddened them by E(B-V)$=0.15$ (Appendix A) with a Fitzpatrick & Massa (2007) reddening model using K. Barbary’s extinction python code131313https://extinction.readthedocs.io/en/latest/index.html. We then took the Case M5 RG template (Fluks et al., 1994) and scaled it by eye to match the red side of the quiescent spectrum (which leads to the smoothest residuals in both spectra), then subtracted this constant RG spectrum from the observed spectra before integrating the optical luminosity. The pre-subtraction dereddened spectra and the scaled RG template are shown in Fig. 12.

Our best guess accretion rate for the 2016 peak is 6 $\times 10^{-7}$ (d/2.5 kpc)2 (R / 0.01 R*⊙) (0.9 M⊙/ M) M⊙* yr*-1*, obtained as follows. The peak optical (3300–7930Å) luminosity in the 2016 outflow fast state was about 1100 L*⊙* (d/2.5 kpc)2. The peak Swift NUV flux (2100–2870Å) was about 70 L*⊙* (d/2.5 kpc)2, but judging by the shape of the NUV spectrum (see § 3.4), we suspect that the intrinsic unabsorbed 1200–3300Å range had a luminosity similar to the April 29 NUV+FUV IUE spectra, or about 730 (d/2.5 kpc)2 L*⊙*. Our total peak accretion disc luminosity is therefore about the same as that of the 1990 outburst (Tomov et al., 1992; Schmid et al., 2001).

We estimate a quiescent accretion rate minimum value of about 1 $\times$ 10*-7*(d/2.5 kpc)2 (R / 0.01 R*⊙) (0.9 M⊙/ M) M⊙* yr*-1*. This comes from about 205 (d/2.5 kpc)2 L*⊙* observed from 3300–7930Å in our RG-subtracted optical spectrum from 2015 March 8, 65 (d/2.5 kpc)2 L*⊙* observed from 1200–3300Å in the 1993 February 23 SWP and LWP spectra observed by the IUE during a quiescent period with similar optical flux, and an estimated 30 (d/2.5 kpc)2 L*⊙* absorbed by the Fe ii curtain. We ignore the possibility that up to half the luminosity could be emitted by the boundary layer in the extreme-ultraviolet and not re-emitted, which if true would double our estimated accretion rates.