A bipolar structure and shocks surrounding the stellar-merger remnant V1309 Scorpii

Context. V1309 Sco is an example of a red nova, a product of the merger between non-compact stars. V1309 Sco is particularly important within the class of red novae due to the abundance of photometric data of the progenitor binary before the merger. Aims. We aim to investigate the spatio-kinematic and chemical properties of the circumstellar environment, including deriving the physical conditions and establishing the origins of the di ff erent circumstellar components. Methods. We use radiative transfer modelling of molecular emission in sub-mm spectra to examine the properties of the molecular gas, and use forbidden line diagnostics from optical spectra to constrain electron density and temperature using forbidden line diagnostics. We compare line intensities from shock models to observations to look for and constrain shocks. Results. We derive a new kinematical distance of 5.6 kpc to the source. The detection of ro-vibrational H 2 and sub-mm HCO + emission in 2016 and 2019, respectively, indicate active shock interactions within the circumstellar environment. The velocity profiles of both H 2 and HCO + , as well as the moment-1 maps of sub-mm CO and 29 SiO, indicate a bipolar structure that may be asymmetric. The sub–mm and optical molecular emission exhibits temperatures of 35–113 and 200 K, respectively, whilst the atomic gas is much hotter, with temperatures of 5–15 kK, which may be due to shock heating. Conclusions


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 (∼2000 K), multipeaked light curves with an intermediate peak luminosity between those typical of classical and supernovae, and high molecular abundances and dust.Sources such as V838 Mon (Bond et al. 2003;Kamiński et al. 2015), V4332 Sgr (Kamiński et al. 2010;Kamiński & Tylenda 2011, 2013;Tylenda et al. 2015), BLG-360 (Tylenda et al. 2013) and CK Vul (Kamiński et al. 2020(Kamiński et al. , 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 Ktype 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 ⋆ email: thomas@ncac.torun.pl 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 two years 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 2 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-millimeter/millimeter spectrum of V1309 Sco obtained with the Atacama Large sub-Millimeter Array (ALMA), indicating a velocity gradient in the moment-1 maps 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 detect 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 d 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 oxygenbearing molecules.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 nonstandard 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 analysis of ALMA and XSHOOTER 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 XSHOOTER, as well as details on the modelling of optical molecular bands and statistics of H 2 shock models used to examine the shock properties in the CSE.

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

Observations
V1309 Sco was observed with ALMA in band 7 on 17 th January, 26 th March and 9 th April 2016.The details of these observations are described in detail by K18.Band 7 observations were repeated on 26 th -28 th 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 hours, 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 13 CN (4-3) line (ν rest =345.34GHz), 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-mm source at RA=17 h 57 m 32 s .6768,Dec=−30 • 43 ′ 14 ′′ .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.

Line identification
To identify the sub-mm/mm spectral features, we used the CAS-SIS spectroscopic analysis tool1 (Vastel et al. 2015) to extract spectroscopic information from the Cologne Database for Molecular Spectroscopy (CDMS; Endres et al. 2016) and the Jet  A.1). From the moment-0 maps of CO (3-2), the flux has increased from 181.5 Jy/(beam km s −1 ) in 2016 to 226.0 Jy/(beam km s −1 ) 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 2 10(4,6)-10(3,7) line in 2016 and could not be identified.In 2019 the strength of HCO + relative to SO 2 increased.This observation, coupled with the better spectral resolution of the 2019 ALMA spectrum, allowed HCO + to be (partially) resolved from SO 2 .
The ALMA spectrum is dominated by SO 2 emission, with 20 individual lines identified.The average LSR peak velocity across all identified lines in Table A.1 is -81 km s −1 .

Source size and distance
Using moment-0 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 (0 ′′ .51×0 ′′ .45) 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 e =232.5 km s −1 .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 d (=3.464×10 8 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 2 of 5σ (=0 ′′ .097) and adopting this as the source size, θ, we derive the distance d = ( 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 ∼540 AU.

Moment maps
As the dominant emission originates from the CO (3-2) and 29 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 −1 .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 north-east 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 29 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 29 SiO (8-7) moment-1 map shows an identical structure to the observed velocity gradient in the 28 SiO (8-7) line from 2016 (Kamiński et al. 2018, their Fig.9).At Solar composition, 29 Si constitutes only 5% of the main isotope 28 Si, and thus 29 SiO emission should be optically thinner than 28 SiO (cf. Sect. 2.4.3).The velocity structure of the 29 SiO (8-7) emission shows a redshifted 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 Here, σ refers to the standard deviation of the Gaussian profile fitted to the source, calculated using σ=FWHM/2 √ 2 ln(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 34 SO 2 .The initial model fitting was performed only on SO 2 using a Monte-Carlo Multi-Chain (MCMC) χ 2 fitting algorithm, simulating 20 detected SO 2 lines to get the best estimates of the SO 2 kinetic temperature and column density.The upper energy level range of the SO 2 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 29 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 2 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 −1 and has a larger line width of 233 km s −1 .NC is at a velocity of -123 km s −1 and has a FWHM of 112 km s −1 .
The detection of HCO + is important for recognizing the late role of outflows in the CSE.However, the HCO + (4-3) line is heavily blended with SO 2 .We use the model of SO 2 , shown in Table 3 and Fig. 3 to recover the intrinsic line profile of HCO + by subtracting the SO 2 model from the ALMA 2019 spectrum.The profile is further discussed in Sect.4.2.

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, 28 SiO, and 29 SiO was calculated as where N is the column density, T is the gas temperature, ∆v is the FWHM of the modelled molecular lines, ν 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.

XSHOOTER
Here, we describe optical observations of V1309 Sco taken in 2016 and 2022.Notes.Excitation temperatures were measured as the peak brightness temperature of the spectral feature.The values used for both isotopologues of SiO, such as temperature and partition function, are the same, with only column density varying between the calculations.

2016
V1309 Sco was observed in 2016 using XSHOOTER (Vernet et al. 2011), a mid-resolution spectrometer on UT 2 of the Very Large Telescope (VLT) which allows spectra to be taken simultaneously across three arms (UVB, VIS, NIR) with a total spectral coverage of 3000-24800 Å.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 th May 2016, runs B and C on 23 rd June 2016, and run D on 6 th July 2016.During the first run the telescope was incorrectly centered on a background star, and so run A was repeated on 6 th July 2016.In total, 4 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 XSHOOTER 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 per- formed assuming E(B-V)=0.8mag (Kamiński et al. 2015).Table 5 presents a summary of the XSHOOTER observations.

2022
XSHOOTER observed V1309 Sco again on 26 th 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 spectra were reduced and combined using the same pipeline version and method described in Sect.3.1.1.The standard star used to calibrate the spectra was Hip110573.

Line identification and analysis
The XSHOOTER 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. 2009a(Kamiński et al. , 2015;;Tylenda et al. 2015).In order to confirm uncertain lines, the NIST database 3 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 XSHOOTER arm in the 2016 and 2022 spectra.

UVB range
The signal-to-noise ratio (SNR) 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 2 Σ + -X 2 Σ + (1,1) band, the flux measurement on Hβ is not reliable.Other significant features are the [S II] 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 3 https://physics.nist.gov/PhysRefData/ASD/lines-form.html 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-ExoCross4 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.

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 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 4 Π-X 4 Σ − (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 1 -F 1 and TiO γ (2,1) F 3 -F 3 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 Pgopher5 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 2 -F 2 band, as it had sufficient SNR and was devoid of saturation effects, unlike the TiO γ (0,0) F 3 -F 3 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 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-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.

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 10290 and 10310 Å 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.

2022 spectrum
The 2022 XSHOOTER spectrum shows several changes in the spectroscopic appearance of the V1309 Sco remnant.The 2019 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-10400 Å, and the [Fe II] λλ5159,5262 doublet.Additionally, the V I 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 2019.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.

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.Could the CSE 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-10000 K at intervals of 1000 K.We estimated the best fitting electron temperature using a χ 2 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 χ 2 test yields a temperature of 3500 ± 250 K, with χ 2 =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(Luridiana et al. , 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 diagnos-tics, we adopted the [S II] ratios from Keenan et al. (1996) and the [Fe II] ratios from Bautista et al. (2015).Figure 5 shows the [S II] 4068 6716+6731 vs 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] 5262 8617 vs 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] 5262 8617 vs 7155 8617 .When we compare these constraints with the diagnostic grid of both [Fe II] 5159 8617 vs 5159 5262 (Fig. 6, middle panel) and [Fe II] 5262 8617 vs 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 ≤ 15000 K, consistent with the observed quantities and diagnostic grids of [S II] 4068 6716+6731 vs 6716 6731 and [Fe II] 5262 8617 vs 7155 8617 .However, as 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 was corrected for interstellar reddening assuming E(B-V)=0.8mag (Tylenda et al. 2011;Kamiński et al. 2015).In order to test this assumption, we used the python-implemented tool extinction6 .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.

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 XSHOOTER 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

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. 2017Iaconi et al. , 2020;;Soker 2023).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, 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 recognize 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 fea-tures 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 centrallypeaked 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 3-dimensional (plus kinematical) 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 modeling of the sub-mm 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 29 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-mm 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 Sec.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 −1 , 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-mm lines of CO and 29 SiO which represents the coolest gas observed.The profiles overlap closely, although the sub-mm 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-mm molecular gas.
Therefore attempts to resolve the remnant are likely to be more successful using sub-mm observations rather than optical/IR.Better angular resolution is also required to resolve the overall structure of the V1309 Sco remnant.

Role of shocks
In analogy with other Galactic red nova remnants (Kamiński et al. 2009b(Kamiński et al. , 2010;;Tylenda et al. 2015), the remnant star of V1309 Sco is expected to be a cool (2000-3000 K) giant or supergiant.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 2 and sub-mm HCO + emission suggests that the material may be heated by active shocks propagating through the CSE, subsequently changing its molecular composition.

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 2 , the radiative transfer model of SO 2 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 −1 , and its profile partially overlaps in velocity with the profiles of other sub-mm and optical lines.The HCO + emission is the most redshifted emission detected across both the optical and sub-mm.The shift between HCO + and the average velocity of the atomic gas is 81 km s −1 .
HCO + is often considered a shock tracer, especially in cooler environments.Shock dissociation of stable molecules such as H 2 O and CO, and subsequent recombination can form HCO + in circumstellar media (Pulliam et al. 2011;Sánchez Contreras et al. 1997, 2000, 2015).Alternatively, ultraviolet and X-ray fluxes can influence HCO + abundances through photodissociation (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 candidate for the formation of HCO + .

H 2
Ro-vibrational 1-0 S(1) and Q(1) lines of H 2 are detected in our XSHOOTER NIR spectrum at 2.10 and 2.41 µm, respectively.Ro-vibrational H 2 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 XSHOOTER spectrum.In Fig. 7, we plot both detected H 2 line profiles compared to sub-mm and optical lines.The line profiles of each line are found at the same peak LSR velocity of -216 km s −1 , and partly overlap with the velocity profiles of CO, Fe I and [Fe II].The blue wings of both H 2 lines do not overlap at all with other observed species.The shift between the average atomic gas velocity and the velocity of the H 2 lines is 82 km s −1 .
In order to derive shock parameters from H 2 , we examine the H 2 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 2 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 XSHOOTER 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 −8 ).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 2 emission (Fig. F.1), although the external UV radiation and H 2 cosmic-ray ionisation rate cannot be constrained.The models indicate a pre-shock density lower limit of 10 7 cm −3 .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 (v s <10 km s −1 ) shock velocities are typically found at lower magnetic fields, whereas higher shock velocities (20-30 km s −1 ) 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.

The properties of the shocks
We have proposed that both NIR H 2 and sub-mm HCO + emission probes shock-excited regions in V1309 Sco, but -as shown in Fig. 7 -these shock tracers have very different kinematics with H 2 observed at extreme blueshifted flow and HCO + in the extreme redshifted flow.Fig. 7 shows an almost equal magnitude in the velocity shift of the shock tracers (∼80 km s −1 ) 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 2 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 2 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-mm emission is not affected by dust extinction.
The results presented in Appendix F show that the shocks that excite H 2 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 −1 or are faster than 15 km s −1 .Highvelocity J-type shocks with shock velocities greater than 20 km s −1 would be expected to dissociate H 2 (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 2 in the HCO + excitation region may be due to H 2 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 2 , 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.

Effects of dust obscuration
In the XSHOOTER 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 (CC-SNe 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, asymmetrries 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. can be applied to V1309 Sco.

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 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 XSHOOTER 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.

Summary
Using optical spectroscopy from 2016 and sub-mm interferometry from 2019, we examine 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 consider two components present in the line profiles of many sub-mm molecular lines, except for CO and SO which exhibit more complex line profiles.Radiative transfer modelling reveal 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 2 , which is detected via ro-vibrational lines.Shock models indicate that the blueshifted lobe shocks, where H 2 is excited, can be classified as J-type, and have a pre-shock density n H ≥ 10 7 cm −3 .The detection of only HCO + in the redshifted lobe, and ro-vibrationally excited H 2 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 suggest gas densities of 3 ≤ log(N e ) ≤ 5, and an electron temperature of 5 000 ≤ 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 29 SiO support the presence of molecular outflows.Using sub-mm SO 2 emission as a thermometer, we derive 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(Kamiński et al. , 2021;;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 do not resolve, either spatially or spectroscopically.Further observations at higher angular resolution in sub-mm may reveal the puzzling nature of the outflows.Due to decreasing signal-to-noise in optical spectra, James Webb Space Telescope observations may serve as a better method to fully constrain the outflow and shock properties.

Fig. 1 :
Fig. 1: ALMA band 7 spectra of V1309 Sco.The reference frame is the local standard of rest (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 SO 2 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.

Fig. 2 :
Fig. 2: Top: Moment-1 map of the 2019 CO (3-2) line, extracted across -310 to 160 km s −1 .Black contours represent 5, 10, 20, 40, 80 and 95% of the peak CO (3-2) flux.Bottom: Moment-1 map of the 29 SiO (8-7) line, extracted across -260 to 160 km s −1 .Black contours represent 5, 10, 20, 40, 80 and 95% of the peak 29 SiO (8-7) flux.The colourbars cover the same velocity range for both maps.Pixels with values below 5σ noise level were blanked in all maps.(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 TableC.1, with multiple bands detected for several molecules.The VIS spectrum is the only band where we observe ScO (A-X band Fig. 3: ALMA band 7 spectrum (black) with the best fitting LTE model overlaid (red).

Fig. 4 :
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-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.
For [S II], we used 4068 6716+6731 vs 6716 6731 , and replaced the 4068 Å line with the other detected [S II] lines (except the [S II] λλ6716,6731 doublet) to produce 5 different grids.For [Fe II], we plotted the 5159 7155 All line diagnostic grids are presented in logarithmic scale.

Fig. 5 :
Fig. 5: Diagnostic grid of [S II] 4068 6716+6731 vs 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 solid and dotted lines, and the red point and errorbars indicate the location of the observed quantity for both line ratios.
Fig. 8: Schematic of the bipolar structure hypothesised in V1309 Sco.The velocities for the blueshifted and redshifted lobes are measured from the 2019 ALMA data.

Fig. 9 :
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. 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. 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 band covers a wavelength range of 1-1.4 µm, and the K s band covers 1.8-2.6 µm.Our 2016 XSHOOTER spectra presented in Fig. 10 in those spectral ranges show that the J 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

Table 1 :
Summary of ALMA observations in 2016 and 2019.

Table 2 :
Parameters of SO 2 gas for both ALMA epochs.

Table 3 :
Column densities (in cm −2 ) of best fitting model for each species constrained by the LTE modelling in CASSIS.

Table 5 :
Summary of XSHOOTER observations of V1309 Sco.

Table A .
1: Molecular transitions identified in the spectrum of V1309 Sco over 2 epochs.Notes.a The HCO + line is covered by epoch 1 observations, but is not resolved from the SO 2 10(4,6)→10(3,7) line.Rest frequency, transition quantum numbers and upper level kinetic temperature (E u ) are specified, as taken from the JPL and CDMS catalogues accessible via CASSIS.
Table B.1: Identified atomic and H 2 emission lines in 2016 XSHOOTER spectrum.