© The Authors 2024

1 Introduction

Red novae are a class of stellar eruptions believed to be caused by the merger of two non-compact stars (Soker & Tylenda 2003, 2006; Tylenda & Soker 2006; Soker et al. 2007). Red novae are characterised by low temperatures (~3000 K), multi-peaked light curves with an intermediate peak luminosity between those typical of classical novae and supernovae¹ (Kasliwal 2012), and high molecular abundances and dust. Sources such as V838 Mon (Munari et al. 2002; Bond et al. 2003; Tylenda 2005; Kamiński et al. 2015; Mobeen et al. 2021), V4332 Sgr (Martini et al. 1999; Banerjee et al. 2003; Kamiński et al. 2010; Kamiński & Tylenda 2011, 2013; Tylenda et al. 2015), BLG-360 (Tylenda et al. 2013), and CK Vul (Banerjee et al. 2020; Kamiński et al. 2020, 2021) have all been previously identified as red novae, the latter of which is the oldest known red nova and had previously been misidentified as a classical nova. Multiple red novae have also been discovered in other galaxies (Pastorello et al. 2019, and references therein).

V1309 Sco erupted in 2008 (Nakano 2008) and remains the only example of a Galactic red nova whose progenitor has been regularly observed. The OGLE survey (Udalski 2003) photometrically observed the progenitor binary system over multiple years, revealing that the binary system was made up of two K-type sub-giants with masses of (0.5–0.9) and (1.1–1.3) M_⊙. Its orbital period of P ≈ 1.4 days was exponentially decreasing when approaching the common envelope (CE) phase and the merger event (Stepień 2011; Tylenda et al. 2011; Pejcha 2014). This is the only example to date of a binary system of non-compact stars observed during the spiral-in phase of a merger event. V1309 Sco therefore provides a useful testing ground for establishing connections between the properties of the progenitor binary and the merger remnant.

The nature of the stellar remnant is currently unknown due to steady formation of dust obscuring the star, starting from 2yr after the merger event (Nicholls et al. 2013; McCollum et al. 2013), but the origins and nature of the circumstellar environment (CSE) surrounding the remnant are of interest. In particular, we wish to understand the evolution of post-merger remnants and find connections between the kinematic structure and chemical markers associated with such events. The CSE has been predicted to come from a variety of sources: L₂ mass loss before the merger, circumbinary disk formation, common envelope (CE) phase, merger ejecta and stellar winds (e.g. Zhu et al. 2013; Nandez et al. 2014; Pejcha 2014; Pejcha et al. 2017; Matsumoto & Metzger 2022; MacLeod & Loeb 2020; MacLeod et al. 2022).

Several observational studies attempted to infer the structure of the remnant. Kamiński et al. (2018; hereafter K18) revealed a two-component structure in the sub-millimetre/millimetre spectrum of V1309 Sco obtained with the Atacama Large sub-Millimeter Array (ALMA), indicating a velocity gradient in the moment-1 maps (intensity-weighted velocity field)² of SiO and showing a similar structure to another Galactic red nova, V4332 Sgr, which was shown to exhibit bipolar outflows. Mason & Shore (2022) independently detected two absorption components in spectra taken soon after the 2008 eruption, which they associate with material ejected both before and during the coalescence. One of these components switches from absorption and emission 25–35 days after the eruption, showing that the circumstellar material evolved rapidly at this time.

The chemical composition is also worthy of study. V1309 Sco is known to be oxygen-rich, and so is abundant in oxygen-bearing molecules (Mason et al. 2010). Kamiński et al. (2015) reported the detection of CrO in V1309 Sco for the first time, a molecule only detected in V4332 Sgr and in no other classes of objects. The presence of such a rare molecule in red nova remnants indicates that non-standard chemistry occurs in such objects.

This study aims to gain a better understanding of the kinematic and spectroscopic structure of the remnant of V1309 Sco after 2012. In Sects. 2 and 3, we present observations, results, and an analysis of ALMA and X-shooter observations, respectively. Section 4 presents our discussions of the results and Sect. 5 presents our conclusions. Also provided are several appendices showing tables of detected lines and bands in both ALMA and X-shooter, as well as details on the modelling of optical molecular bands and statistics of H₂ shock models used to examine the shock properties in the CSE.

2 ALMA

We start with presenting the most recent ALMA observations of V1309 Sco.

2.1 Observations

V1309 Sco was observed with ALMA in band 7 on 17 January, 26 March, and 9 April 2016. The details of these observations are described in detail by K18. Band 7 observations were repeated on 26–28 August 2019 (PI: T. Kamiński). The setup used in 2019 had longer baselines of up to 3.6 km and science integration time of ~2.5 h, compared to maximum baselines of 460 m and a science integration time of ~48 min in 2016. The longer baselines meant a spatial resolution of almost one order of magnitude better than in 2016. A summary of the comparison between the observational epochs is shown in Table 1. The 2019 baselines gave a clean primary beam size (FWHM) of 80×60 mas at uniform weighting, whilst the spectral coverage was 342.5–358.1 GHz with a gap between ~346–355 GHz due to the heterodyne setup. The spectral windows were shifted from the setup used for the 2016 observations by ~2 GHz (see Table 1) in order to cover the H¹³CN (4–3) line (v_rest = 345.34 GHz), which was subsequently not detected. The ALMA data was reduced with CASA pipeline version 4.5.3 using the default calibration script. The source was too faint for self-calibration. The data was imaged using the CASA routine tclean.

As noted in K18, the ALMA field of view covers an unidentified sub-millimetre source at $\text{RA}={17}^{\text{h}}{57}^{\text{m}}{32}^{\text{s}}.6768$, $\text{ Dec }=-{30}^{°}{43}^{\prime }{14.}^{\prime \prime }157$ (J2000). It is located ~5″ south-west from V1309 Sco and is bright in the 2019 observations. This source is likely an uncategorised background galaxy.

2.2 Line identification

To identify the sub-millimetre/millimetre spectral features, we used the CASSIS spectroscopic analysis tool³ (Vastel et al. 2015) to extract spectroscopic information from the Cologne Database for Molecular Spectroscopy (CDMS; Endres et al. 2016) and the Jet Propulsion Laboratory database (JPL; Pickett et al. 1998). All identified lines are listed in Table A.1. As the ALMA spectrum is mostly unchanged in 2019 compared to 2016 (see Fig. 1), the majority of lines seen in the 2019 spectrum could be identified without major effort. However, due to the aforementioned shift in the spectral windows for the 2019 observations relative to 2016, the SO₂ 20(0,20)–19(1,19), 25(3,23)–25(2,24), SO 9(8)– 8(7), and ²⁸SiO 8–7 transitions are not observed in 2019. Instead, ²⁹SiO 8–7 and SO₂ 12(4,8)–12(3,9) are the newly observed lines, not covered in 2016 (Table A.1). From the moment-0 maps of CO (3–2), the flux has increased from 181.5 Jy/(beamkm s⁻¹) in 2016 to 226.0 Jy/(beamkm s⁻¹) in 2019 (26% increase).

Perhaps the most significant is the identification of HCO⁺ for the first time in V1309 Sco (Table A.1). As seen in Fig. 1, the HCO⁺ (4–3) line is blended with the SO₂ 10(4,6)–10(3,7) line in 2016 and could not be identified. In 2019 the strength of HCO⁺ relative to SO₂ increased. This observation, coupled with the better spectral resolution of the 2019 ALMA spectrum, allowed HCO⁺ to be (partially) resolved from SO₂.

The ALMA spectrum is dominated by SO₂ emission, with 20 individual lines identified. The average Local Standard of Rest (LSR) peak velocity across all identified lines in Table A.1 is −81 km s⁻¹.

2.3 Source size and distance

Using moment-0 (integrated intensity) maps combined from the full measurement sets, we used the CASA routine imfit to fit elliptical Gaussians to the spatial distribution of the source in each epoch. The images used represent all emission within all spectral windows. The resulting source sizes (deconvolved from the main beam) measured were 178 (±21) × 125 (±36) and 87.9 (±1.8) × 59.8 (±1.7) mas respectively for epochs 1 and 2. As a sanity check, the images were cleaned using natural and uniform weightings to compare the results, with the different weightings providing the same results per epoch within 1σ errors. As each epoch had different beam sizes, we smoothed the epoch 2 data to the same beam size $\left(0{.}^{″}51\times 0{.}^{″}45\right)$ as epoch 1. The imfit results provided the same results as those previously found for the unsmoothed epoch 2 data. At face value, comparing the source size in each epoch would indicate that the source has contracted between 2016 and 2019. However, we find this highly unlikely and rather consider the uncertainties in epoch 1 measurements to be largely underestimated.

Using the CO (3–2) emission from the 2019 observations, we estimate the kinematical distance of V1309 Sco. The CO (3–2) distribution has a semi-major axis of 45.5 mas, assumed to be the semi-major FWHM/2 of the 2D Gaussian fitted to the CO (3–2) spatial distribution using imfit. We assume that the radial velocity is equal to the velocities observed in the far wings of the CO (3–2) line, and so take half the baseline width as the tangential velocity in one direction, which turned out to be v_e = 232.5 km s⁻¹. The light curve peak of V1309 Sco occurred on 2.5 September 2008 (Julian date 2 454 712; Tylenda et al. 2011) and the first date of the epoch 2 observations was 26 August 2019. We assume that the CO gas was ejected at the point of the merger, and so t=4012 days (=3.464×10⁸ s) had elapsed. This gives a physical radius of the CO emission region. This also provides an upper limit on the distance, as the cold CO gas may have been ejected prior to the merger, meaning the time elapsed is potentially much longer and the CO emission radius is therefore larger. Converting the $FWHM/2=0{.}^{″}0455$ to a value⁴ of $5\sigma (=0{.}^{″}097)$ and adopting this as the source size, θ, we derive the distance d = (v_e × t)/ tan(θ). The kinematic distance we determine is therefore 5.6 kpc, which is consistent with the distance estimate presented in Mason & Shore (2022). This corresponds to a distance of the most extended CO gas from the central star of ~540AU.

Fig. 1 ALMA band 7 spectra of V1309 Sco. The reference frame is LSR. Black shows the 2016 spectrum, and blue shows the 2019 spectrum. The key lines are indicated by the labels. Red vertical lines indicate the position of the multiple identified SO2 lines shifted to the observed positions. The 2019 spectrum was extracted from data smoothed to the same beam size as the 2016 data. Both spectra represent the entire source. Both spectra are smoothed.

2.4 Kinematics and structure

2.4.1 Moment maps

As the dominant emission originates from the CO (3–2) and ²⁹SiO (8–7) lines, both were used to probe the spatio-kinematic structure in 2019. Moment-0 and 1 maps were constructed using the CASA routine immoment. We extracted the emission across the same velocity range of −310 to 160 km s⁻¹. Any pixels with values <5σ than the root-mean-square (rms) noise were neglected. The CO moment-1 map from 2019 is shown in Fig. 2 (top panel). Even though the overall source structure is still largely unresolved in 2019, the better spatial resolution reveals a velocity gradient spanning from northeast to south-west, with the north-east region moving away from the observer. The gradient axis approximately matches the direction of a slight elongation seen in the CO (3–2) moment-0 maps.

We also constructed moment-0 and 1 maps of the ²⁹SiO (8–7) line (Fig. 2, bottom panel). The moment-0 map also shows a nearly point-like structure, as seen in CO (3–2). The ²⁹SiO (8–7) moment-1 map shows an identical structure to the observed velocity gradient in the²⁸ SiO (8–7) line from 2016 (Kamiński et al. 2018, their Fig. 9). At Solar composition, ²⁹Si constitutes only 5% of the main isotope ²⁸Si, and thus ²⁹SiO emission should be optically thinner than²⁸ SiO (cf. Sect. 2.4.3). The velocity structure of the ²⁹SiO (8–7) emission shows a red-shifted lobe in the north-east, with less extreme velocities at extended distances. This is inconsistent with a rotating sphere, which would show the most extreme velocities in the moment map at the most extended distances.

2.4.2 Radiative transfer modelling

To characterise the molecular gas observed with ALMA, we used the local thermodynamic equilibrium (LTE) radiative transfer tool available in CASSIS. We modelled all detected molecules (shown in Table A.1) as well as ³⁴SO₂.

The initial model fitting was performed only on SO₂ using a Monte-Carlo Multi-Chain (MCMC) χ¹ fitting algorithm, simulating 20 detected SO₂ lines to get the best estimates of the SO₂ kinetic temperature and column density. The upper energy level range of the SO₂ transitions cover 48-521 K, slightly wider than the range covered by the 2016 observations. This should result in better constraints of the excitation temperature in 2019. We then assume that emission of all detected molecules is produced in the same LTE conditions, meaning we assume the same temperature for all molecules for each gas component. We model two gas kinematic components for all molecules except for²⁹ SiO, HCO+, and CO, which are better represented by a single Gaussian. The key change made to fit the model was the gas column density for each molecule. The results for the SO₂ fitting and the overall radiative transfer modelling are presented in Tables 2 and 3, respectively, and the best model is shown in Fig. 3.

We therefore abbreviate these two components to the broad component (BC) and narrow component (NC). BC is found at the less negative velocity of −43 km s⁻¹ and has a larger line width of 233 km s⁻¹. NC is at a velocity of −123 km s⁻¹ and has a FWHM of 112km s⁻¹. As the epoch 1 dataset was re-analysed independently, the column densities for epoch 1 presented in Table 3 are different to those found in K18. As the values presented in K18 and in this paper are of the same order of magnitude, any discrepancy does not change the overall conclusions drawn from this analysis.

The detection of HCO⁺ is important for recognising the late role of outflows in the CSE. However, the HCO⁺ (4–3) line is heavily blended with SO₂. We use the model of SO₂, shown in Table 3 and Fig. 3 to recover the intrinsic line profile of HCO⁺ by subtracting the SO₂ model from the ALMA 2019 spectrum. The profile is further discussed in Sect. 4.2.

Fig. 2 Moment-1 maps of molecules detected in V1309 Sco with ALMA, extracted across -310 to 160 km s-1. Top: Moment-1 map of CO (3-2). Bottom: Moment-1 map of29 SiO (8-7). Black contours represent 5, 10, 20, 40, 80 and 95% of the peak flux for both molecules. Pixels with values below 5σ noise level were blanked on all maps.

2.4.3 Optical depth

The CASSIS software applies a limited correction on line saturation effects. It is important to know the opacity of the strongest emission features, since we use these to trace the source structure. Therefore, the optical depth (τ) for CO, ²⁸SiO, and ²⁹SiO was calculated as $$\tau =\frac{c^3}{8\pi v^3}\frac{1}{\text{Δ}v}\frac{N}{Z}{A}_{ul}{g}_u{e}^{\frac{-F_{\Psi }}{T}}\left[e^{\frac{hv}{kT}}-1\right],$$(1)

where N is the column density, T is the gas temperature, Δv is the FWHM of the modelled molecular lines, v is the rest wavelength and all other constants are associated quantum values of each transition (Goldsmith & Langer 1999). All measurable values (N, T, Δv) were obtained from the LTE radiative transfer modelling described in Sect. 2.4.2. The results are presented in Table 4.

The determined values of τ indicate optically thick emission. However, in the wings of the lines, which probe the outflowing gas shown in Fig. 2, the optical depth is considerably lower compared to the values given in Table 4. Therefore, the optically thin wings of the lines represent well the source structure.

Fig. 3 ALMA band 7 spectrum (black) with the best fitting LTE model overlaid (red).

3 X-shooter

Here, we describe optical observations of V1309 Sco taken in 2016 and 2022.

3.1 Observations

V1309 Sco was observed in 2016 using X-shooter (Vernet et al. 2011), a mid-resolution spectrometer on UT 2 of the Very Large Telescope (VLT) that allows spectra to be taken simultaneously across three arms (UVB, VIS, NIR) with a total spectral coverage of 3000–24 800 Å. The spectra were acquired in stare mode. The slit width was $1{.}^{″}3$ for the UVB arm and $0{.}^{″}9$ for the VIS and NIR arms. The pixel binning for all arms was 1×1. The data was observed in four separate runs. Run A was executed on 20 May 2016, runs B and C on 23rd June 2016, and run D on 6 July 2016. During the first run the telescope was incorrectly centred on a background star, and so run A was repeated on 6 July 2016.

In total, four exposures were acquired in the UVB and VIS arms and 20 exposures for NIR. The exposure times used were 2935, 2840, and 600 s for the UVB, VIS and NIR arms respectively. The data was reduced using X-shooter pipeline version 3.6.1 (Modigliani et al. 2010) in the Esoreflex environment, and further processed using standard IRAF routines. The data was calibrated using several standard stars: Hip084982, HD190285, Hip076069, Hip067973, Hip094378, Hip093049, Hip11900, Hip08254, and Hip017734. Dereddening was performed assuming E(B−V)=0.8 mag (Kamiński et al. 2015). Table 5 presents a summary of the X-shooter observations.

X-shooter observed V1309 Sco again on 26th October 2022, using a similar spectral setup to that used in 2016. However, in 2022, a narrower slit of $0{.}^{″}6$ was used for the NIR spectral arm, as the NIR emission was expected to have increased from the previous observation epoch. The exposure times were reduced, and more exposures in total taken. The total exposure time in each arm in 2022 is therefore significantly less than in 2016. The standard star used to calibrate the spectra was Hip110573.

3.2 Line identification and analysis

The X-shooter spectrum from 2016 is dominated by atomic and molecular emission. Atomic lines were identified by comparison to spectra in other Galactic red novae (Kamiński et al. 2009, 2015; Tylenda et al. 2015). In order to confirm uncertain lines, the NIST database⁵ was used to look for lines of the same multiplet, as well as compare relative intensities. Identified atomic lines are listed in Table B.1.

Examining the presence of molecules in the spectra, several oxygen-bearing molecules are identified: AlO, CrO, ScO, TiO, and VO. Multiple electronic systems have been identified across the three spectral arms.

Flux and width measurements of the 2016 observations were done using splot, as part of the onedspec package in IRAF. Errors are estimated at 5σ, where 1σ is equal to the noise rms, and multiplying by the square root of the width of the line in pixels. Sections 3.2.1–3.2.3 describe in detail the detected atomic and molecular emission from each X-shooter arm in the 2016 and 2022 spectra.

3.2.1 UVB range

The signal-to-noise ratio (S/N) of the UVB spectral range is lower than for the VIS and NIR ranges. Despite this, we still identified 26 individual atomic emission lines of different species, including S, Fe, Ca, and Cr. Amongst the strongest lines identified were those of Ca I λ4226, Cr I λ4254, and Fe I λ4375. The Hβ emission was weaker than expected relative to other lines, but since Hβ is partially blended with the AlO B²σ⁺−X²σ⁺ (1,1) band, the flux measurement on Hβ is not reliable. Other significant features are the [SII] and [Fe II] lines, which are valuable due to their sensitivity to electron density and temperature. The Mn I triplet at ~4028 Å and the Cr I doublet at ~5295 Å could not be deblended and so flux measurements in Table B.1 are presented as integrated flux measurements across the blended features. These measurements are shown for completeness and were not used in any detailed analysis.

The UVB range features molecular emission from all observed molecules except ScO, although less dominant than in the VIS and NIR spectra. For AlO and CrO, we see the B–X electronic systems; for VO we see C–X; TiO is traced through the α band emission. All such features at λ > 5000 Å could be modelled via the ExoMol-ExoCross⁶ tool (Yurchenko et al. 2018), whereas very few molecular bands were reliably identified at shorter wavelengths (see Fig. D.1). As seen in Table C.1, very few molecular bands were identified in the UVB spectrum. Although there is a weak continuum seen in the UVB range, it is not strong enough to derive the spectral type and so is not analysed.

3.2.2 VIS range

The VIS range has a similar number of atomic lines identified to UVB with 25 lines, although 6 are uncertain. One such line, found at a peak wavelength of 8662.33 Å, was originally identified as Ca II λ8662. but it was removed due to the inconsistent peak velocity with nearby atomic lines and other identified Ca II lines. The reasons for ambiguity in identification were mostly due to the observed central wavelengths not reliably matching the velocity measured for the majority of the lines, as well as due to possible contamination by molecular emission. For example, the unidentified line at 8601Å is heavily contaminated by the broad VO B⁴Π−X⁴σ⁻ (0,1) band.

We see multiple forbidden lines of [O I], [S II], [Ca II] and [Fe II] as well as lines of Na I, K I, and Rb I. The strongest lines within the VIS range is the K I λλ7664, 7698 doublet. We also see the semi-forbidden [Ca I line at 6572 Å and the Li I 6708 Å line (for a summary of lithium in red nova remnants, see Kamiński et al. 2023). The forbidden lines are important tools for constraints on the properties of the CSE and so are analysed in detail in Sect. 3.4.

The [S II] λ6730 line is contaminated by the TiO γ (1,0) F₁− F₁ and TiO γ (2,1) F₃−F₃ bands. As [S II] λ6730 is important for the constraints of the electron temperature and density, an accurate flux measurement for this line is important. We used the Pgopher⁷ tool (Western 2017) to model the TiO emission across a short wavelength range spanning the [S II] λ6730 line as well as the TiO γ (0,0) bands in order to remove the molecular contamination. TiO was modelled to best fit the TiO γ (0,0) F₂−F₂ band, as it had sufficient S/N and was devoid of saturation effects, unlike the TiO γ (0,0) F₃−F₃ band. The best-fitting temperature was 240 K, with a Gaussian smoothing kernel of 3.5 Å. Pgopher was used rather than ExoMol–ExoCross for convenience, but is less accurate as it is parameterised by a single (rotational) temperature. The resulting flux measurements of the recovered [S II] λ6730 line is shown in Table B.1.

The molecular emission is far richer in the VIS range than the UVB. We see emission from all molecules listed in Table C.1, with multiple bands detected for several molecules. The VIS spectrum is the only band where we observe ScO (A–X band between 6000 and 6100 Å). At λ < 7300 Å the molecular spectrum is dominated by TiO and CrO whereas at longer wavelengths, VO is dominant. The key detections in the VIS spectrum are the various bands of the CrO B–X system, which has only been observed in one other astrophysical source (V4332 Sgr, Tylenda et al. 2015). V1309 Sco and V4332 Sgr remain the only sources in which we detect CrO emission. Additionally, as in Kamiński et al. (2015), we detect emission from the CrO A′–X band observed between 8350 and 9500 Å (Table C.1). Unlike in early spectroscopic observations of the remnant (Kamiński et al. 2015), there is no visible continuum in the VIS range.

3.2.3 NIR range

As in 2012 (Kamiński et al. 2015), there are very few atomic features seen in the NIR range. Ten atomic emission lines are seen, made up of [S II] and 6 unidentified lines. Three lines between 10 290 and 10 310 Å cannot be identified, despite being relatively strong. They do not appear to be sky lines due to the broad widths comparable to the widths of the detected lines.

The molecular emission in the NIR range is just as rich as that seen for the VIS range. The NIR spectrum is dominated by the AlO A–X electronic bands, whilst CrO A–X, VO A–X, and ro-vibrational emission of AlO (X–X) are also seen. The strongest observed molecular feature across all spectral arms in absolute flux units, observed between 12 100 and 12 800 Å, is the blended feature of the AlO v=9−0 ro-vibrational emission, AlO A²Π_1/2−X²σ⁺ (3,0), AlO A²Π_3/2−X²σ⁺ (5,1), and CrO A⁵σ⁺− X⁵Π (0,0) bands. The ro-vibrational H₂ 1–0 S(1) and 1–0 Q(1) lines are detected near 2.1 and 2.4 µm, respectively.

3.3 2022 spectrum

The 2022 X-shooter spectrum shows several changes in the spectroscopic appearance of the V1309 Sco remnant. The 2016 and 2022 spectra are compared in Figs. E.1 and E.2. A significant difference is the absence of many atomic lines, including several forbidden iron and sulphur lines. Many of the atomic lines have also decreased significantly in intensity. Some of the more notable absences are the Rb I λ7947 line, the [S II] multiplet expected between 10280 and 10 400 Å, and the [Fe II] λλ5159,5262 doublet. Additionally, the VI line at 3980.5 Å was observed in the UVB spectrum that had not been previously seen in 2016. The weakening of the optical spectrum can be attributed to increased extinction due to growing dust column densities, indicating continuing dust formation since 2016. Therefore, constraining the physical properties of the optical emitting region of the remnant is becoming increasingly difficult for V1309 Sco, and is not attempted using the 2022 observations.

3.4 Atomic line diagnostics

V1309 Sco shows no photospheric signatures from the merger remnant due to heavy circumstellar extinction caused by dust. The only visible source of emission we can study is the surrounding CSE, the origins of which are uncertain. The CSE could be made up of material ejected during the binary interaction phase preceding the merger, or of material carried out by winds or outflows associated with the coalescence. By using diagnostic line ratios, we aim to examine the physical structure of the remnant.

In order to examine whether LTE conditions apply to the atomic emission region, we compared selected line ratios of Fe I to simulated intensity ratios from the NIST service. The NIST-simulated ratios are calculated in LTE conditions, using the Boltzmann distribution and the Saha equation.

In Fig. 4, we compare observed and simulated intensities of the detected Fe I lines relative to Fe I λ4216. The first set of simulations were extracted between 1000-10 000 K at intervals of 1000 K. We estimated the best fitting electron temperature using a χ² test on each model for the three line ratios, which gave a value of 3000 K. We then reduced the temperature range to 3000-4000 K with an interval of 250 K to refine the constraints. Figure 4 shows the results for the refined grid. Theχ² test yields a temperature of 3500 ± 250 K, with χ² = 0.979. The poor consistency between the 3500 K model and the observed line ratios for the 4376 Å and 5166 Å lines may indicate a non-LTE environment.

Using PyNeb v1.1.16 (Luridiana et al. 2013, 2015), we attempted to get a more robust constraint on the electron temperature and density. The PyNeb tool makes no underlying assumptions about LTE, allowing the physical conditions to be constrained in non-LTE conditions. We constructed diagnostic maps similar to Figs. 2 and 3 in Keenan et al. (1996). For our diagnostics, we adopted the [S II] ratios from Keenan et al. (1996) and the [Fe II] ratios from Bautista et al. (2015). For [SII], we used $\frac{4068}{6716+6731}$ vs. $\frac{6716}{6731}$, and replaced the 4068 Å line with the other detected [SII] lines (except the [SII] λλ6716,6731 doublet) to produce five different grids. For [Fe II], we plotted the $\frac{5159}{7155}$, $\frac{5262}{7155}$, $\frac{5159}{8617}$ and $\frac{5262}{8617}$ ratios against the $\frac{5159}{5262}$ and $\frac{7155}{8262}$ ratios. All line diagnostic grids are presented in logarithmic scale.

Figure 5 shows the [SII] $\frac{4068}{6716+6731}$ vs. $\frac{6716}{6731}$ diagram. The location of the observed ratios on the grid suggests a lower limit on the electron density of log(N_e) ≤ 3.5. No constraint can be derived for the electron temperature T_e. We find that the grids involving the NIR [S II] lines are degenerate and therefore do not contribute to the constraints. Thus, little can be revealed about V1309 Sco by the [S II] diagnostics alone.

The top panel of Fig. 6 shows the diagnostic grid of [Fe II] $\frac{5262}{8617}$ vs$\frac{7155}{8617}$. The position of the observed ratios and associated errors are indicative of an N_e range of 3 ≤ log(N_e) ≤ 5.5, as well as a lower T_e limit of T_e ≥ 5000 K. The N_e constraint inferred from Fig. 6 is consistent with the lower limit inferred by Fig. 5. However, the constraints are loose as the parameter space is densely filled, and so constraints with lower uncertainty cannot be derived from [Fe II] $\frac{5262}{8617}$ vs. $\frac{7155}{8617}$.

When we compare these constraints with the diagnostic grid of both [Fe II] $\frac{5159}{8617}$ vs. $\frac{5159}{5262}$ (Fig. 6, middle panel) and [Fe II] $\frac{5262}{8617}$ vs $\frac{5159}{5262}$ (Fig. 6, bottom panel), the constraints can be tightened for both N_e and T_e. The position of the observed line ratios and associated errors infer 4 ≤ log(N_e) ≤ 5 and 5000 ≤ T_e ≤ 15 000 K, consistent with the observed quantities and diagnostic grids of [S II] $\frac{4068}{6716+6731}$vs. $\frac{6716}{6731}$ and [Fe II] $\frac{5262}{8617}$vs. $\frac{7155}{8617}$. Howeveras the 1σ errorbars extend out between certain grid lines, we provide a more conservative constraint for N_e to include the possibility that the errors are underestimated and due to the fact that the resolution of the diagnostic grids are quite low. The conservative constraint therefore is extended to 3 ≤ log(N_e) ≤ 5.

As described in Sect. 3.1, the data were corrected for interstellar reddening assuming E(B–V) = 0.8 mag (Tylenda et al. 2011; Kamiński et al. 2015). In order to test this assumption, we used the python-implemented tool extinction⁸. We removed the previously applied extinction from the data and applied new corrections with different values of E(B–V) to see if the constraints are affected. We re-reddened the data with E(B–V)=0.6 and 1.0 mag and plotted these values on the diagnostic grids. There was no significant effect of different E(B–V) values on the diagnostics we derive.

To constrain the properties of shocks present in the CSE, we used the MAPPINGS V code (Sutherland et al. 2018). We compared the measured integrated fluxes of atomic emission lines available in the MAPPINGS V model database to the models themselves. The models were primarily dependent on the gas density and shock velocity. Lines modelled in MAPPINGS include all the forbidden [S II] and [Fe II] lines, as well as [O I] and Fe I. However, the errors on the line fluxes proved too large to constrain either the gas density or shock velocity. Another issue was that for the densities we explored, which were dictated by the constraints on the gas density established from the PyNeb forbidden line diagnostics (Figs. 5–6), the shock models were difficult to resolve at different densities for many of the modelled emission lines.

Fig. 4 NIST-simulated intensities of the Fe I lines at 4375, 5110 and 5166 Å relative to the Fe I λ4216. The simulations were calculated at temperatures between 3000 and 4000 K, at intervals of 250 K. The χ2 test showed that the best fitting temperature was 3500 K (green points). The black points indicate the observed line ratios.

Fig. 5 Diagnostic grid of [S II] $\frac{4068}{6716+6731}$ vs. $\frac{6716}{6731}$. Solid lines represent lines of constant temperature and dotted lines represent lines of constant density. Larger black dots show the intersecting points between solidand dotted lines, and the red point and errorbars indicate the location of the observed quantity for both line ratios.

3.5 Molecular modelling

Using the ExoMol-ExoCross modelling described in Appendix C, we attempted to constrain the rotational (T_rot) and vibrational temperatures (T_vib) of all detected molecules in our X-shooter spectra. The results are T_rot = 200 and 300 K and T_vib = 1700 and 2400 K for AlO and TiO, respectively. The fitting routine could not converge for CrO, ScO, and VO, and so we instead take the average ro-vibrational temperatures fitted for AlO and TiO, T_rot = 200 K and T_vib = 2100 K, and apply these to the other three detected molecules to obtain their simulated spectra. The overall model, as well as individual contributions from each molecule, is presented in figures of Appendix D.

Fig. 6 Same as Fig. 5, but for [Fe II] $\frac{5262}{8617}$ vs. $\frac{7155}{8617}$ (top panel), [Fe II $\frac{5159}{8617}$vs. $\frac{5159}{5262}$(middle panel) and $\frac{5262}{8617}$vs. $\frac{5159}{5262}$(bottom panel).

4 Discussion

4.1 CSE structure

Very little theoretical work has been done on the post-merger evolution of red novae, with studies concentrating primarily on the merging process (e.g., Zhu et al. 2013; Pejcha 2014; Pejcha et al. 2017; Nandez et al. 2014; Metzger & Pejcha 2017; Iaconi et al. 2017, 2020; Soker 2023a). Therefore, any connections between theoretical and observational studies of red novae cannot be fully verified. Nevertheless, the angular momentum of the binary is expected to impose some form of bipolar structure in the remnant (Soker 2018, 2023b), but whether it could be observed at any stage and how much mass it carries are open issues.

Our ALMA maps only marginally resolve the cool component of the V1309 Sco remnant. The spatial structure is not easy to recognise at the achieved resolution, but a slight elongation can be seen in mapped molecular lines. Additionally, moment-1 maps show a velocity gradient along the same axis. These features are consistent with a bipolar configuration. A large part of the molecular gas must also be located near the systemic (average) velocity and near the central source to produce the centrally-peaked moment-0 maps shown with contours in Fig. 2. Such observational characteristics are consistent with what ALMA observed for V4332 Sgr, the older and larger clone of V1309 Sco (K18). Using ALMA data of V4332 Sgr, which better resolved the source, K18 constructed a full 3D (plus kine-matical) model of the remnant, which prominently shows a pair of bipolar lobes formed by opposite wide-angle outflows. We postulate that the remnant of V1309 Sco displays very similar architecture.

This bipolar structure can be considered to be reflected in our radiative transfer modelling of the sub-millimetre lines, in which the emission features are represented by two kinematical components, NC and BC. Since these components are not readily separated in the spectra and overlap for all observed species, deciding on modelling two components is merely a simplification which was thought to produce more robust constraints on the gas properties than any single-component model. Fitting more than two components would require even more parameters, and is less informative. However, the observations may suggest that more than two kinematic components form the observed spectral profiles. Indeed, the line profiles of CO and ²⁹SiO representing the entire source show an excess within the red wings of the lines that is not replicated in our CASSIS model (Fig. 3). This excess is also clearly seen in Fig. 7. There is some asymmetry in the gas distribution, which we are not able to resolve with the current sub-millimetre data.

In our simplifying model, the two postulated lobes are represented by the narrow and broad Gaussian components, whose properties for both ALMA epochs are displayed in Table 2. The automatic fit is not consistent in assigning the positions and widths of the components when the two epochs are compared.

In the earlier epoch, it is the blue component which is wider, whereas in the 2019 fit the blue component is narrower. It is unlikely to be a real change in the physical parameters of the flow and rather reflects the large uncertainties in modelling the multiple overlapping components. As discussed in Sect. 2.4.3, the moderate to high optical depths of the mapped emission may be another source of uncertainty. A FWHM of 150 km s⁻¹, representing an average of all simulated components, is probably the best educated guess on the actual (projected) velocity dispersion in both lobes.

The bulk of material in the CSE is likely represented by the molecular emission observed with ALMA. This material has cooled to temperatures of 35–113 K, and recombined to form molecules. Slightly warmer molecular gas at a temperature of 200 K is traced by the optical spectra presented in Appendix A. Here, we assume that the rotational temperature is close to the excitation temperature, whilst the vibrational temperature of 2100 K represents the colour temperature of the radiation field responsible for the fluorescent emission in the electronic molecular bands. Atomic gas traced by the same optical spectra exhibits much higher temperatures of 5–15 kK, revealing hot gas with partial ionisation.

In Fig. 7, we compare the velocity profiles of Fe I and [Fe II], tracing the hot gas, with those of the sub-millimetre lines of CO and ²⁹SiO which represents the coolest gas observed. The profiles overlap closely, although the sub-millimetre lines show stronger red wings relative to the optical lines at positive velocities. It appears that the kinematics of the hot and cold gas regions are similar within the remnant, but this does not mean that they are chemically mixed. In classical cool circumstellar envelopes (such as those surrounding AGB stars), warmer gas is usually found closer to the central star which heats the gas via radiation or shocks (Bell 1993; Olguin et al. 2020; Massalkhi et al. 2020). It is therefore likely that the atomic gas is much more compact than the sub-millimetre molecular gas.

Therefore attempts to resolve the remnant are likely to be more successful using sub-millimetre observations rather than optical/IR. Better angular resolution is also required to resolve the overall structure of the V1309 Sco remnant.

Fig. 7 Velocity profiles of HCO+ (4–3), H2 1–0 S(1) and H2 1–0 Q(1), normalised or otherwise scaled and shifted (both panels). The top panel also shows the CO (3–2) and Fe I λ5110 profiles, whilst the bottom panel shows [Fe II] λ8616 and 29SiO (8–7). The grey dashed line represents the average peak velocity of the atomic emission detected in X-shooter at−134 km s−1.

4.2 Role of shocks

In analogy with other Galactic red nova remnants (Kamiński et al. 2009, 2010; Tylenda et al. 2015), the remnant star of V1309 Sco is expected to be a cool (2000–3000 K) giant or super-giant. Section 4.1 describes the detection of atomic gas with temperatures exceeding 5000 K. The presence of hot gas in the vicinity of a cool star would therefore be significant. The gas may be material that was heated during the 2008 eruption, cooling on a timescale of decades and recombining into molecular gas. However, the detection of NIR H₂ and sub-millimetre HCO⁺ emission suggests that the material may be heated by active shocks propagating through the CSE, subsequently changing its molecular composition.

4.2.1 HCO⁺

In Sect. 2.2, we postulated the identification of the HCO⁺ J=4−3 line in the 2019 spectrum, which must have been much weaker or absent in the 2016 spectrum (see Fig. 1). Although the HCO⁺ line is heavily blended with SO₂, the radiative transfer model of SO₂ is good enough to extract the intrinsic line profile of HCO⁺. This profile is shown in Fig. 7.

The recovered HCO⁺ profile has a peak velocity of −53 km s⁻¹, and its profile partially overlaps in velocity with the profiles of other sub-millimetre and optical lines. The HCO⁺ emission is the most redshifted emission detected across both the optical and sub-millimetre. The shift between HCO⁺ and the average velocity of the atomic gas is 81 km s⁻¹.

HCO⁺ is often considered a shock tracer, especially in cooler environments. Shock dissociation of stable molecules such as H₂O and CO, and subsequent recombination can form HCO⁺ in circumstellar media (Pulliam et al. 2011; Sanchez Contreras et al. 1997, 2000, 2015). Alternatively, ultraviolet and X-ray fluxes can influence HCO⁺ abundances through photo-dissociation (Kimura et al. 2012; Cleeves et al. 2017). Due to a lack of UV or X-ray sources in V1309 Sco, shocks are the most likely explanation for the formation of HCO⁺.

4.2.2 H₂

Ro-vibrational 1–0 S(1) and Q(1) lines of H₂ are detected in our X-shooter NIR spectrum at 2.10 and 2.41 µm, respectively. Ro-vibrational H₂ lines have long been considered a shock tracer, and have often been used in conjunction with NIR [Fe II] to probe different regions of the same shock (Sternberg 1989; Davis et al. 2003; Kokusho et al. 2020; Mohan et al. 2023). We detect no such NIR [Fe II] in our X-shooter spectrum. In Fig. 7, we plot both detected H₂ line profiles compared to sub-millimetre and optical lines. The line profiles of each line are found at the same peak LSR velocity of −216km s⁻¹, and partly overlap with the velocity profiles of CO, Fe I and [Fe II]. The blue wings of both H₂ lines do not overlap at all with other observed species. The shift between the average atomic gas velocity and the velocity of the H₂ lines is 82 km s⁻¹.

In order to derive shock parameters from H₂, we examine the H₂ shock model grid presented in Kristensen et al. (2023). The grid is calculated using the Paris-Dunham code (Godard et al. 2019), covering six parameters; pre-shock density, shock velocity, transverse magnetic field strength, external UV radiation, H₂ cosmic-ray ionisation rate and fractional abundance of polycyclic aromatic hydrocarbons (PAHs). The authors identify the twenty-five dominant cooling lines, eight of which are covered in our X-shooter spectra. Of those eight, three lines lie within the telluric absorption bands (see Fig. D.2) and so are unreliable in estimating upper limits of fluxes. Three more lines, 1–0 Q(3), 1–0 S(2), and 1–0 Q(5) are not detected. The 1–0 S(2) line is heavily obscured by AlO, and so we measure the flux of the AlO feature and apply the integrated flux as a conservative upper limit. The other two lines are not detected above rms noise levels. We therefore simulated the fluxes of a Gaussian with peak flux equal to the rms noise in the region where we expect the line to be detected, with the same FWHM as the average of the two detected lines. Using the two detected and three undetected lines, we calculate the observed line ratios of the 1–0 S(1) line and three upper limits against 1–0 Q(1), the stronger of the two detected lines. We then search through the model grid to find models which are in agreement with the observed ratios to 5σ accuracy for the detected line ratio (where σ is the line ratio error calculated using standard quadrature equations), and are below the ratio upper limits calculated for the simulated lines.

Nicholls et al. (2013) find no evidence of PAHs in their analysis of the mid-IR SED of V1309 Sco. Therefore, we only examine models with the lowest fractional abundance of PAHs modelled (=10⁻⁸). However, we find it has no real effect on the statistics presented in Appendix F. Using the conditions described above, we find that 139 out of 14364 models are in agreement. We present statistics of the consistent models in Appendix F.

We find some initial constraints on the shocks related to the H₂ emission (Fig. F.1), although the external UV radiation and H₂ cosmic-ray ionisation rate cannot be constrained. The models indicate a pre-shock density lower limit of 10⁷ cm⁻³. This is a lower limit as larger densities are not covered by the models, but may still be the true pre-shock density. The scaled value of the magnetic field covers 0.1 to 1.0 (while the grid spans to values as high as 10). For the shock velocity, we see two separate subsets of models at lower and higher shock velocities (Fig. F.1, third panel). Models with low (4v_s < 10 km s⁻¹) shock velocities are typically found at lower magnetic fields, whereas higher shock velocities (20–30 km s⁻¹) are found across a wide range of magnetic field strengths (Fig. F.2). Finally, almost all models fulfilling our observational constraints exhibit J-type shocks.

4.2.3 The properties of the shocks

We have proposed that both NIR H₂ and sub-millimetre HCO⁺ emission probes shock-excited regions in V1309 Sco, but – as shown in Fig. 7 – these shock tracers have very different kinematics with H₂ observed at extreme blueshifted flow and HCO⁺ in the extreme redshifted flow. Figure 7 shows an almost equal magnitude in the velocity shift of the shock tracers (~80 km s⁻¹) relative to the average velocity of the atomic emission. This symmetric location can be explained if the emission arises in two opposite flows, as schematically illustrated in Fig. 8. Since the fastest ejected gas is expected farthest from the central object, we place the shock regions at the apexes of the bipolar flows, which are seen in cool molecular gas observed with ALMA. However, it is unclear why the blueshifted lobe produces only ro-vibrational H₂ emission, whereas we only see redshifted HCO⁺ emission. This suggests that the shocks propagating in opposite directions interact with the ambient circumstellar gas in very different ways. This could depend on either the properties of the ambient gas into which the shocks are propagating, or even the shocks themselves. It could also be that the NIR H₂ emission excited by shocks via the redshifted lobe is compact and could be obscured by circumstellar dust. However, this does not explain why we do not see blueshifted HCO⁺ emission. The sub-millimetre emission is not affected by dust extinction.

The results presented in Appendix F show that the shocks that excite H₂ could be classified as J-type shocks. The results also suggest that the shock velocities lie in one of two regions, in the range of 2–10 km s⁻¹ or are faster than 15 km s⁻¹. High-velocity J-type shocks with shock velocities greater than 20 km s−¹ would be expected to dissociate H₂ (Flower et al. 2003, and references therein), but this may not be the case for lower velocity shocks. If the shocks associated with HCO⁺ do have higher shock velocities, the absence of ro-vibrational H₂ in the HCO⁺ excitation region may be due to H₂ being efficiently dissociated. On the other hand, if the shocks propagating via the blueshifted lobe had lower shock velocities, it is possible that the shocks have sufficient energy to excite H₂, but not enough to dissociate stable molecules such as CO that would go on to form HCO⁺.

An alternative explanation would be that there is an increased column density of CO surrounding the redshifted lobe that could lead to increased HCO⁺ formation. Indeed, as mentioned, the CO 3–2 line profile does show excess emission in the red wing.

Fig. 8 Schematic of the bipolar structure hypothesised in V1309 Sco. The velocities for the blueshiſted and redshiſted lobes are measured from the 2019 ALMA data.

4.3 Effects of dust obscuration

In the X-shooter spectra, we detect a plethora of neutral and singly ionised emission, as well as forbidden line emission (see Sect. 3.4). It has been shown that the line profiles of atomic emission in core-collapse supernovae is affected by dust (Lucy et al. 1989; Pozzo et al. 2004; Bevan et al. 2017), with a bluewards shift in the peak velocity and asymmetries between the red and blue wings. Shore et al. (2018) applied a similar method to the ejecta of classical novae, which have simpler geometries and lower ejecta velocities and masses. Similar effects are seen, as well as that in cases where the inclination is low so that the line of sight barely passes through the ejected material, asymmetries will still be seen. To examine the effect of dust on the neutral and ionised atomic emission, we looked at the profiles of neutral and singly ionised Fe and Ca. The profiles are shown in Fig. 9. The average profiles of [Ca II] and [Fe II] in 2016 show a slight asymmetry skewed towards shorter wavelengths, with the peak velocity also blueshifted. This would indicate that ionised gas is more obscured by dust and so is likely to be found at smaller radial distances from the central star, assuming the winds or outflows are biconical and sufficiently inclined to the observer that both cones are visible. Such a geometry is consistent with the ALMA velocity maps shown in Fig. 2, which in 2019 are similar to that seen for V4332 Sgr (Kamiński et al. 2018). In K18, it is noted that the inclination angle may be even lower than that found for V4332 Sgr (13°). If this is the case, the assumed geometry in the simulations of Shore et al. (2018) can be applied to V1309 Sco.

4.4 V1309 Sco is not a blue straggler

The VISTA Variables in Vía Lactea (VVV Minniti et al. 2010) survey observed V1309 Sco across multiple epochs between 2010 and 2015. Analysis by Ferreira et al. (2019) showed that the near-infrared J – K_s colour decreased from 1.40 to 0.42 mag between 2010 and 2015. In addition, the I – K_s colour (where I is the OGLE I-band, see Tylenda et al. 2011) changed from 3.54 mag in 2010 to 2.75 mag in 2015. Ferreira et al. (2019) use the apparent shift towards bluer colours between 2010 and 2015, as well as the asymptotic decline in the K_s band, to conclude that V1309 Sco is a blue straggler (Sandage 1953). The J and covers a wavelength range of 1–1.4 µm, and the K_s band covers 1.8–2.6 µm. Our 2016 X-shooter spectra presented in Fig. 10 in those spectral ranges show that the / filter is dominated by a prominent molecular band featuring emission from AlO and CrO. Kamiński et al. (2015) note that in 2009, V1309 Sco spectra showed molecular absorption bands which diminished towards the end of the year. By 2012, much of the molecular absorption had transitioned to pure emission in their spectra. It seems as though the 2010 VVV observations were taken during this transition period, with the 2015 observations taken after the molecular emission had been established. As the J band is dominated by molecular emission rather than continuum, the apparent increase in blue luminosity in the near-infrared is attributed to the emergence of molecular emission rather than the evolution of V1309 Sco to a blue straggler. Regarding the I-K_s colour evolution, the K_s band is less dominated by molecular emission, although some AlO A – X emission is seen at the edge of the spectral coverage. As seen in Fig. 10, no continuum is detected within the K_s band. Our X-shooter observations are much more sensitive and yet no continuum is detected that can be attributed to V1309 Sco. V1309 Sco is located in a crowded field, and so it is possible for mistakenly associating continuum emission with V1309 Sco. Kamiński et al. (2015) note that the continuum source disappears between 2009 and 2012, meaning that the evolution of the I-K_s colour can be attributed to a disappearance of continuum entirely. We therefore believe that V1309 Sco is not evolving towards a blue straggler.

Fig. 9 Averaged profiles of neutral and singly ionised Ca and Fe. Black and red lines show Fe and Ca respectively, whilst solid and dashed lines show neutral and singly ionised line profiles. The profiles of ionised species show asymmetry and a shift towards bluer velocities, possibly indicating greater obscuration by dust.

Fig. 10 X-shooter spectra between 11300 and 14100 Å from 2016 (black) with VVV transmission curves overplotted. The top panel shows the transmission curve of the J band (blue), and the bottom panel shows the transmission curve of the Ks band (red).

5 Summary

Using optical spectroscopy from 2016 and sub-millimetre interferometry from 2019, we examined the circumstellar environment of the stellar merger remnant V1309 Sco in order to understand the kinematical structure of the environment and any recent changes. We considered two components present in the line profiles of many sub-millimetre molecular lines, except for CO and SO which exhibit more complex line profiles. Radiative transfer modelling reveals that these two components have physically distinct properties, including column density, temperature, and line width. We associate the broad and narrow components with the redshifted and blueshifted lobes of outflowing gas, respectively, forming a bipolar structure. These lobes interact with pre-existing ambient material surrounding the central star through shocks, inducing non-uniform chemistry within the circumstellar environment. Via shocks, the redshifted lobe dissociates molecules such as CO to form HCO⁺, whereas the blueshifted lobe excites H₂, which is detected via ro-vibrational lines. Shock models indicate that the blueshifted lobe shocks, where H₂ is excited, can be classified as J-type, and have a pre-shock density n_H ≥ 10⁷ cm⁻³. The detection of only HCO⁺ in the redshifted lobe, and ro-vibrationally excited H₂ only in the blueshifted lobe, indicates different properties of the shocks, and therefore likely different kinematics in either the outflows themselves or in the properties of the ambient medium.

The diagnostics of the atomic emission suggests gas densities of 3 ≤ log(N_e) ≤ 5, and an electron temperature of 5000 ≤ T_e ≤ 15 000 K, which is unexpectedly high around a star with an expected effective temperature of around 3000 K. Modelling of the molecular bands present in the optical spectra reveal low rotational temperatures of 200–300 K. The moment-1 maps of CO and ²⁹SiO support the presence of molecular outflows. Using sub-millimetre SO₂ emission as a thermometer, we derived the excitation temperature of the molecular gas to be 35–113 K. It likely represents the coolest emission regions of the remnant.

We therefore present V1309 Sco as a kinematically and chemically complex object. The presence of a bipolar outflow presents an analogy with other Galactic red novae such as V4332 Sgr and CK Vul (Kamiński et al. 2020, 2021; Zain Mobeen et al. 2023). The inconsistencies in distribution of the shock tracers and the kinematics of the bipolar lobes may indicate additional components to the circumstellar environment that we did not resolve, either spatially or spectroscopically. Further observations at higher angular resolution in sub-millimetre may reveal the puzzling nature of the outflows. Due to decreasing S/N in optical spectra, James Webb Space Telescope observations may serve as a better method to fully constrain the outflow and shock properties.

Acknowledgements

Based on observations collected at the European Southern Observatory under ESO programme 097.D-009 (PI T. Kamiński). This paper makes use of the following ALMA data: ADS/JAO.ALMA2018.1.00336.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. T.S. and T.K. acknowledge funding from grant SONATA BIS no 2018/30/E/ST9/00398 from the Polish National Science Center (PI T. Kamiński). We thank K. Menten, M.Z. Mobeen, R. Tylenda and the rest of the Mergestars team for their invaluable discussions.

Appendix A ALMA line table

Appendix B X-shooter atomic line table

Appendix C Molecule identification and modelling

Identification of molecular bands relied on comparison to previous studies (Kaminski et al. 2015) as well as using the interactive molecular modelling tool Pgopher¹⁰ (Western 2017) to model the detected molecules. Although the identification of molecular lines and bands are more difficult than for atomic lines, the specific shape of some molecular bands help to identify them.

We generated models from the molecular line lists available in 2022 on the ExoMol online database¹¹ (Patrascu et al. 2015; McKemmish et al. 2016, 2019; Chubb et al. 2021). The models are calculated using both the vibrational (T_vib) and rotational temperature (T_rot) as well as the quantum properties for each ro-vibrational state. ExoMol files are available for AlO, TiO, and VO. We also modelled ScO and CrO using custom-made files generated from published spectroscopic constants of the electronic bands covered (see Kamiński et al. (2009) for details). For each molecule, a model grid was generated between 100 and 3500 K for both T_rot and T_vib at intervals of 100 K, with the condition that T_rot < T_vib. As there are regions of the spectrum where certain molecules are not observed, the models were not generated across the entire spectrum for all molecules to save time. AlO was modelled from 5000 to 25000 Å, TiO was modelled from 4880 to 8700 Å and VO is modelled between 3400 and 16000 Å. The simulations were then fitted via minimised residuals to find the best fitting rotational and vibrational temperatures.

The fitting routine was carried out across prominent features of each molecule. The fitting region for AlO from 14200 Å to 17400 Å covered multiple AlO A²π−X²Σ⁺ bands. For TiO, the fitting region we used was from 6110 Å to 7260 Å, covering multiple TiO γ features including the prominent TiO γ (0,0) band between 7040 and 7160 Å. For VO, we used the region between 7090 and 9010 Å to fit the VO models to the spectrum. This covered multiple VO B⁴π−X⁴Σ⁻ bands in the VIS spectrum. The CrO B⁵π −X⁵π band was modelled across 6380 and 6500 Å. Although this range is narrow, it is the only CrO B⁵π −X⁵π emission that is uncontaminated from atoms or other molecules. CrO A⁵σ⁺−X⁵π bands were modelled across the main CrO NIR feature between 12185 and 12655 Å. Finally, ScO was modelled across the only system detected, A²π−X²σ⁺, prominent between 5610 and 6340 Å. The full ExoMol model, together with the individual components of each molecule, is presented in Figs. D.1–D.2.

Due to issues with the fitting of VO, CrO, and ScO, we are unable to present reliable constraints on the rotational and vibra-tional temperatures of these molecules. The output of the fitting routine returned the maximum possible T_vib of 3500 K that we simulated in our model grid, no matter what input parameters were used. We extended the model grids to higher temperatures, but only achieved convergence at unrealistic temperatures of ~5000 K. Furthermore, upon visual examination, the fit residuals were larger than those at T_vib ≤ 3500 K for many prominent features. Therefore, we instead present in Figs. D.1–D.2 the VO model simulated with the average ro-vibrational temperatures that we found between TiO and AlO.

The identification of features of the AlO spectrum presented in this paper is improved compared to that presented in Kamiński et al. (2015) owing to the more complete list of lines in the Exo-Mol database. In particular, a feature near 12200 Å marked in Kamiński et al. (2015) as an AlO perturbation turned out to be a ro-vibrational v=9−0 band of AlO X²σ⁺ (within the ground electronic state of AlO).

Appendix D Molecular models

Fig. D.1 Plot of the X-shooter UVB (top) and VIS (bottom) spectra and best-fitting molecular model. The red line shows the full model with all components added, and the green, light blue, gold, magenta, and dark blue lines underneath represent the best-fitting VO, AlO, TiO, ScO and CrO models, respectively. The atmospheric transmission spectrum is represented by the grey curve and the grey regions represent the strongest telluric absorption bands.

Fig. D.2 Same as Fig. D.1 but for the X-shooter NIR spectrum.

Appendix E Epoch comparison

Fig. E.1 X-shooter UVB spectra of V1309 Sco taken in 2016 (black), and 2022 (red). Prominent spectral lines are indicated by vertical grey lines. The Cr I triplet at ~427O Å is highlighted by the yellow area in the top panel. The 2022 spectrum is shifted along the y-axis.

Fig. E.2 Same as Fig. E.1, but for the VIS and NIR spectra in the top and bottom panels, respectively. The yellow area in the bottom panel highlights the four [S II] lines found between 10280-10370Å.

Appendix F H₂ models

Fig. F.1 Histograms of models from Kristensen et al. (2023) in agreement with the observed H2 line ratios and upper limits. From top to bottom, the models show the following: the pre-shock density, scaled transverse magnetic field, shock velocity, and shock type.

Fig. F.2 Histogram of models consistent with low (vs <10 km s−1, blue) and high (vs ≥10 km s−1, magenta) shock velocities relative to the input magnetic field strength.

References

All Tables

All Figures

Fig. 1 ALMA band 7 spectra of V1309 Sco. The reference frame is LSR. Black shows the 2016 spectrum, and blue shows the 2019 spectrum. The key lines are indicated by the labels. Red vertical lines indicate the position of the multiple identified SO2 lines shifted to the observed positions. The 2019 spectrum was extracted from data smoothed to the same beam size as the 2016 data. Both spectra represent the entire source. Both spectra are smoothed.

In the text

	Fig. 2 Moment-1 maps of molecules detected in V1309 Sco with ALMA, extracted across -310 to 160 km s-1. Top: Moment-1 map of CO (3-2). Bottom: Moment-1 map of29 SiO (8-7). Black contours represent 5, 10, 20, 40, 80 and 95% of the peak flux for both molecules. Pixels with values below 5σ noise level were blanked on all maps.
In the text

	Fig. 3 ALMA band 7 spectrum (black) with the best fitting LTE model overlaid (red).
In the text

	Fig. 4 NIST-simulated intensities of the Fe I lines at 4375, 5110 and 5166 Å relative to the Fe I λ4216. The simulations were calculated at temperatures between 3000 and 4000 K, at intervals of 250 K. The χ2 test showed that the best fitting temperature was 3500 K (green points). The black points indicate the observed line ratios.
In the text

	Fig. 5 Diagnostic grid of [S II] $\frac{4068}{6716+6731}$ vs. $\frac{6716}{6731}$. Solid lines represent lines of constant temperature and dotted lines represent lines of constant density. Larger black dots show the intersecting points between solidand dotted lines, and the red point and errorbars indicate the location of the observed quantity for both line ratios.
In the text

	Fig. 6 Same as Fig. 5, but for [Fe II] $\frac{5262}{8617}$ vs. $\frac{7155}{8617}$ (top panel), [Fe II $\frac{5159}{8617}$vs. $\frac{5159}{5262}$(middle panel) and $\frac{5262}{8617}$vs. $\frac{5159}{5262}$(bottom panel).
In the text

	Fig. 7 Velocity profiles of HCO+ (4–3), H2 1–0 S(1) and H2 1–0 Q(1), normalised or otherwise scaled and shifted (both panels). The top panel also shows the CO (3–2) and Fe I λ5110 profiles, whilst the bottom panel shows [Fe II] λ8616 and 29SiO (8–7). The grey dashed line represents the average peak velocity of the atomic emission detected in X-shooter at−134 km s−1.
In the text

	Fig. 8 Schematic of the bipolar structure hypothesised in V1309 Sco. The velocities for the blueshiſted and redshiſted lobes are measured from the 2019 ALMA data.
In the text

	Fig. 9 Averaged profiles of neutral and singly ionised Ca and Fe. Black and red lines show Fe and Ca respectively, whilst solid and dashed lines show neutral and singly ionised line profiles. The profiles of ionised species show asymmetry and a shift towards bluer velocities, possibly indicating greater obscuration by dust.
In the text

	Fig. 10 X-shooter spectra between 11300 and 14100 Å from 2016 (black) with VVV transmission curves overplotted. The top panel shows the transmission curve of the J band (blue), and the bottom panel shows the transmission curve of the Ks band (red).
In the text

Fig. D.1 Plot of the X-shooter UVB (top) and VIS (bottom) spectra and best-fitting molecular model. The red line shows the full model with all components added, and the green, light blue, gold, magenta, and dark blue lines underneath represent the best-fitting VO, AlO, TiO, ScO and CrO models, respectively. The atmospheric transmission spectrum is represented by the grey curve and the grey regions represent the strongest telluric absorption bands.

In the text

	Fig. D.2 Same as Fig. D.1 but for the X-shooter NIR spectrum.
In the text

	Fig. E.1 X-shooter UVB spectra of V1309 Sco taken in 2016 (black), and 2022 (red). Prominent spectral lines are indicated by vertical grey lines. The Cr I triplet at ~427O Å is highlighted by the yellow area in the top panel. The 2022 spectrum is shifted along the y-axis.
In the text

	Fig. E.2 Same as Fig. E.1, but for the VIS and NIR spectra in the top and bottom panels, respectively. The yellow area in the bottom panel highlights the four [S II] lines found between 10280-10370Å.
In the text

	Fig. F.1 Histograms of models from Kristensen et al. (2023) in agreement with the observed H2 line ratios and upper limits. From top to bottom, the models show the following: the pre-shock density, scaled transverse magnetic field, shock velocity, and shock type.
In the text

	Fig. F.2 Histogram of models consistent with low (vs <10 km s−1, blue) and high (vs ≥10 km s−1, magenta) shock velocities relative to the input magnetic field strength.
In the text