Resurrection of Type IIL Supernova 2018ivc: Implications for a Binary Evolution Sequence Connecting Hydrogen-rich and Hydrogen-poor Progenitors

The evolution of massive stars toward core collapse and supernova (SN) explosions has not been fully understood, especially in their post-main-sequence (MS) evolution (e.g., Langer 2012; Maeda 2022). Studying the pre-SN mass-loss history through observational data after the SN explosion sheds light on this issue; it is essentially the one and only method to probe massive star evolution in the final millennium to even less than 1 yr (e.g., Smith 2017). The radio synchrotron emission, as created at the shock wave formed by the interaction between the SN ejecta and the circumstellar matter (CSM), is a direct probe to the CSM distribution (e.g., Chevalier 1998; Chevalier & Fransson 2006; Maeda 2012; Matsuoka et al. 2019). For the typical SN shock velocity of ∼10,000 km s⁻¹ and the CSM velocity of ∼20 km s⁻¹ (for a case of a supergiant progenitor; Smith 2017), the nature of the CSM probed by the synchrotron signal a few years after the explosion translates to the mass-loss rate ∼1000 yr before the explosion.

The SNe IIb are a particularly interesting type of SN in this context. They are defined by the initial appearance of hydrogen lines in their optical spectra that eventually disappear (Filippenko 1997), indicating that they do not have an H-rich envelope as massive as canonical SNe II (red supergiants) while still keeping ∼0.01–1 M_⊙ of the envelope (e.g., Shigeyama et al. 1994; Woosley et al. 1994; Bersten et al. 2012; Hiramatsu et al. 2021), the features that distinguish SNe IIb from the more stripped envelope analogs, i.e., SNe Ib (He stars) or SNe Ic (C+O stars). As such, SNe IIb bridge SNe II and SNe Ib/c in terms of the pre-SN mass-loss mechanism, especially that responsible for ejecting the H-rich envelope (e.g., Fang et al. 2019). The binary evolution channel is a leading scenario for the envelope stripping from SNe II toward SNe IIb/Ib/Ic (e.g., Ouchi & Maeda 2017; Yoon 2017), whose progenitors probably share a similar initial mass range (e.g., Fang et al. 2019). Details are, however, yet to be clarified.

Discovered by the D < 40 Mpc SN survey (Tartaglia et al. 2018) on 2018 November 24 (UT) within a day of the explosion (Valenti et al. 2018; Bostroem et al. 2020), SN 2018ivc is an outlier of the SN IIL subclass with its faint and rapidly evolving optical light curve. Maeda et al. (2023) presented and analyzed the data taken by the Atacama Large Millimeter/submillimeter Array (ALMA) covering 4–198 days since the (well-determined) explosion date, together with the optical and X-ray light curves. They suggested that SN 2018ivc is indeed an analog of an SN IIb explosion (potentially with a slightly more massive H-rich envelope), with its different observational features from canonical SNe IIb originating in its powering mechanism, the SN–CSM interaction. Such a scenario may also be in line with its optical spectral features showing broad and boxy emission lines (Dessart & Hillier 2022). Maeda et al. (2023) further proposed that SN 2018ivc and a fraction of SNe IIL may be a direct link between SNe IIb and SNe IIP in the binary evolution where the sequence of SNe Ic/Ib–IIb–(some) IIL–IIP is mainly controlled by the initial orbital separation. This is along the same line as previously suggested by Nomoto et al. (1995, 1996) but with a substantial difference in the mode of binary interaction (see Section 4 for more details).

The proximity (87 east and 161 north) of SN 2018ivc to the core of the well-studied Seyfert galaxy NGC 1068 (M77), for which we assume a distance of ${10}_{-1.5}^{+1.8}$ Mpc (Tully et al. 2009), requires a high angular resolution to resolve this SN, but at the same time it offers a unique opportunity for long-term monitoring; the SN location may be covered by archival observations targeting the core of NGC 1068. Recently, its detection at ∼1000 days after the explosion in the centimeter emission, at 6.5 GHz by the Karl G. Jansky Very Large Array (VLA), and at 6.3 GHz by the enhanced Multi-Element Radio Linked Interferometer Network was reported (Mutie et al. 2022). In this paper, we report the result of our long-term monitoring of SN 2018ivc with ALMA. The paper is structured as follows. In Section 2, we present the observation and data reduction. Results are presented in Section 3, with the analysis of the nature of the CSM at a large scale around SN 2018ivc. Its implications for the progenitor evolution of SN 2018ivc and its relation to other classes of SNe are discussed in Section 4. The paper is summarized in Section 5 with concluding remarks.

In this paper, we report late-time observations of SN 2018ivc (≳1000 days) based on data taken through ALMA Director's Discretionary Time (DDT) program 2021.A.00026.S (PI: K. Maeda), complemented by the data at earlier epochs published in Maeda et al. (2023) and archival data newly presented in this work. The log of the data used in the present work, together with the measured fluxes, is shown in Table 1. The DDT program was conducted on 2022 August 18 (1364 days since the explosion), essentially in the same spectral setup with our early-phase observations up to ∼200 days. On-target exposure time is 5.0 minutes for band 3 and 9.6 minutes for band 6. The central frequencies in the continuum bands are 100 and 250 GHz. The angular resolution is 081 × 071 in band 3 and 033 × 028 in band 6 (Figure 1). Since this observation targeted the SN, i.e., the point source, the imaging processes are simple. We use the data calibrated through the standard ALMA pipeline with Common Astronomy Software Applications (CASA; The CASA Team et al. 2022). The continuum flux densities associated with the SN in each image are measured based on the imfit task in CASA.

Zoom In Zoom Out Reset image size

Table 1. ALMA Measurements of SN 2018ivc

MJD	Phase a	Frequency b	Fν c	Array	Resolution d	ID	PI
	(days)	(GHz)	(mJy)
Band 3
58,449.1 e	4.1	100.0	4.25 ± 0.22	C43-4	070	2018.1.01193.T	K. Maeda
58,452.1 e	7.1	100.0	7.42 ± 0.38	C43-4	094	2018.1.01193.T	K. Maeda
58,462.1 e	17.1	100.0	9.05 ± 0.46	C43-4	117	2018.1.01193.T	K. Maeda
58,466.0	21.0	93.5	7.41 ± 0.99	C43-4	114	2018.1.01506.S	S. Viti
58,643.6 e	198.6	100.0	0.336 ± 0.026	C43-9	006	2018.A.00038.S	K. Maeda
58,658.5	213.5	92.1	0.322 ± 0.036	C43-9	005	2018.1.01135.S	J. Wang
58,751.2	306.2	103.0	0.264 ± 0.044	C43-6	033	2018.1.01684.S	T. Tosaki
59,469.4	1024.4	92.1	1.128 ± 0.069	C43-9	005	2019.1.00026.S	M. Imanishi
59,809.4	1364.4	100.0	0.955 ± 0.093	C-5	076	2021.A.00026.S	K. Maeda
Band 6
58,449.1 e	4.1	250.0	4.21 ± 0.43	C43-4	028	2018.1.01193.T	K. Maeda
58,451.2 e	6.2	250.0	4.32 ± 0.44	C43-4	042	2018.1.01193.T	K. Maeda
58,462.1 e	17.1	250.0	2.49 ± 0.28	C43-4	050	2018.1.01193.T	K. Maeda
58,643.6 e	198.6	250.0	0.120 ± 0.022	C43-9	003	2018.A.00038.S	K. Maeda
59,809.4	1364.4	250.0	0.286 ± 0.067	C-5	030	2021.A.00026.S	K. Maeda

Notes.

^a The phase is measured from the putative explosion date (MJD 58,445.0). ^b Central frequency. ^c With 1σ error. ^d Average of the major and minor axes. ^e From Maeda et al. (2023).

Download table as:  ASCIITypeset image

In addition, we use ALMA band 3 archival data that cover the position of SN 2018ivc. We selected the data in which the SN is detected at a signal-to-noise ratio (S/N) of >5. We downloaded the processed continuum images by the Japanese Virtual Observatory. In the cases of 2018.1.01135.S and 2019.1.00026.S, we manually performed the imaging processes, i.e., the clean task in CASA with primary beam correction, because the archival images do not cover the SN position. For 2018.1.01135.S, we conducted standard self-calibration to minimize the effects by side lobes from the bright nuclear emission. The S/N was improved by self-calibration, but the flux value associated with the SN is consistent regardless of the self-calibration processes. For 2019.1.00026.S, the data were taken during several visits over the course of ∼2 weeks, and we select those that provide the best S/N at the SN position.

Figure 1 shows the reconstructed images of SN 2018ivc on days 199 (Maeda et al. 2023) and 1364 (this work). After the flux decay in the earlier phase up to ∼200 days, the SN has begun brightening again. There is a hint of a slow decay between days 1024 and 1364, but they are also consistent with the same flux level within 1σ (after correcting for the central frequency difference).

Figure 2 compares the spectral energy distributions (SEDs) on days 199 and 1364, with the addition of the centimeter flux on day 1217 taken from Mutie et al. (2022). It is clear that SN 2018ivc remains optically thin in the high-frequency emission in the ALMA bands. Assuming f_ν ∝ ν^α , α = −1.32 ± 0.15 (1σ) is measured on day 1364. Given the SED slope of α = −1.12 ± 0.22 on day 199 (Maeda et al. 2023), there is a hint of the spectral steepening, but it is also consistent with no change within 1σ (i.e., α = −1.22 ± 0.27 for the combined analysis of the data on days 199 and 1364). The SED as found in our ALMA observations also indicates that the centimeter emission must be in the optically thick regime. This highlights the power of the higher-frequency observation despite the generally lower flux level; catching the optically thin portion of the synchrotron emission is essential to constrain the CSM properties (Maeda et al. 2021, 2023).

Zoom In Zoom Out Reset image size

Figure 3 shows the light curves. Rebrightening at 100 and 250 GHz is clearly seen, which probably started at ∼300 days. Since it has been optically thin in the ALMA bands both before and after the rebrightening (Figure 2), the flux increase means that the CSM is not smoothly distributed, and that the CSM density in the outer region must be higher than the extrapolation from the inner region.

Zoom In Zoom Out Reset image size

We have computed a series of synchrotron emission models (see Maeda et al. 2021, 2023) by considering a CSM extended outward from ∼4 × 10¹⁶ cm (i.e., with the main interaction taking place at ∼500 days for the SN ejecta velocity of ∼10,000 km s⁻¹). In practice, we simply assume that the CSM is distributed with a density structure of ρ_CSM ∝ Dr^−s without a hole in the computation (where $D=\dot{M}/4\pi {v}_{{\rm{w}}}$, $\dot{M}$ is the mass-loss rate, and v_w is the mass-loss wind velocity). Once the interaction has fully developed at ∼500 days (at ∼4 × 10¹⁶ cm), the mass of the swept-up CSM is expected to soon be dominated by the outer component; therefore, the evolution of the SN–CSM interaction will quickly lose the memory of the earlier interaction and follow the asymptotic behavior determined by the outer CSM component. We thus apply this model to the epoch after ∼500 days.

The ejecta structure and microphysics parameters are the same as those adopted by Maeda et al. (2023) for the early-phase modeling: M_ej = 3 M_⊙, E_K = 1.2 × 10⁵¹ erg, and the outer ejecta density slope of n = −10. The centimeter emission is in the optically thick regime, and we find that the synchrotron self-absorption alone would not explain the cutoff at the low frequency. We thus include free–free absorption, which is, however, poorly constrained by theory. We simply vary the optical depth by changing the constant electron temperature in the unshocked CSM. It is just for demonstration, and in most of the subsequent analyses, we do not use constraints from the centimeter observations due to the large uncertainty in the model framework.

Figure 3 shows models with four different CSM distributions: s = −2 (i.e., increasing density), 0 (constant), 1, and 2 (steady-state mass loss). These CSM distributions are plotted in Figure 4. By changing the density scale (D), all of the models can reproduce the flux levels at the latest epochs (1024 and 1364 days). The flat (s = 0) and increasing (s = −2) CSM distributions are not favored, given the nearly flat or slowly decaying evolution observed during this time window. The model with s = 1 provides the best representation of the evolution, while our focus in the present work is the CSM density scale without aiming at constraining the value of s. It is possible to explain the flux levels at 6.5/6.3 GHz in the optically thick emission; within a simple treatment of the free–free absorption, we are, however, not able to explain the flat or slowly decaying evolution also seen in the centimeter emission. We emphasize that the centimeter emission, including both the model and the calibration of the data, is not a topic of this work.

Zoom In Zoom Out Reset image size

While there is a degeneracy in the CSM slope given the limited temporal coverage, the CSM density at ∼10¹⁷ cm is strongly constrained to be ∼10⁻¹⁸ g cm⁻³; this is where the data points exist. It is substantially above the extrapolation from the inner structure, and it is about an order of magnitude larger in terms of the mass-loss rate (assuming that v_w did not change).

In the SN IIb scenario with an extended progenitor (Maeda et al. 2023), the mass-loss velocity would be ∼20 km s⁻¹ (e.g., Groh 2014; Smith 2017). Then, the location of the "outer" dense CSM (∼10¹⁷ cm) corresponds to the mass-loss history ∼1500 yr before the explosion. Based on the binary evolution model toward SNe IIb (Ouchi & Maeda 2017), Maeda et al. (2023) suggested that the progenitor of SN 2018ivc might have experienced an extreme case of the case C mass transfer, where the binary mass transfer is initiated in the advanced evolution stage after the core He ignition (e.g., Langer 2012), to explain the high mass-loss rate just before the explosion as inferred from the early-phase observation (see also Maeda et al. 2015). As compared to other SNe IIb, it is suggested that the strong binary interaction has taken place later in its evolution, closer to the time of the core-collapse explosion, which is attributed to a larger initial orbital separation.

Ouchi & Maeda (2017) indeed predicted that the extreme case C mass transfer may bridge the evolution toward SNe IIb and SNe IIP, where the SNe as surrounded by a dense CSM will be identified as either SNe IIL or IIn. In Figure 4, we show two such models showing late-time binary interaction (Ouchi & Maeda 2017) for the initial masses of the primary and secondary of 16 and 9.6 M_⊙, with the only difference in the initial orbital period (1800 days for model 25 and 1950 days for model 28). In both models, the binary first experiences rapid mass transfer until the mass ratio is inverted, during which most of the H-rich envelope (≳5 M_⊙) is ejected. This is then followed by a relatively high mass-transfer rate of ≳10⁻⁵ M_⊙ yr⁻¹ as the primary keeps filling the Roche lobe (RL). The final leftover envelope mass is ∼2 M_⊙ in these models, which is somewhat large for an SN IIb but qualitatively consistent with the extended SN IIb progenitor scenario with a slightly more massive envelope than canonical SNe IIb; the final envelope mass should be dependent on the binary parameters, especially the mass ratio that roughly controls the mass of the leftover envelope after the rapid-transfer phase (see above). The initial rapid transfer takes place ∼5800 and ∼300 yr before the explosion for models 25 and 28, respectively.

While the timing is different, the mass-loss properties during the initial rapid-transfer phase are similar (Figure 4). It is thus possible, in principle, to have the strong interaction phase at ∼1500 yr where most of the envelope, ≳5 M_⊙, is ejected, which forms the dense CSM at ∼10¹⁷ cm, by a different combination of the component masses and the initial separation. This is demonstrated in Figure 4 by comparing the mass-loss history derived for SN 2018ivc and those predicted by models 25 and 28 of Ouchi & Maeda (2017). Indeed, if we shift the mass-loss history of model 25 toward the explosion date, assuming slightly different binary parameters (e.g., the initial separation and/or initial mass ratio), it can explain that derived for SN 2018ivc reasonably well.

In optical wavelengths, SNe IIb from "extended" progenitors tend to show signs of the SN–CSM interaction in the late phase (≳1 yr) as the emergence of strong H_α emission and light-curve flattening (Patat et al. 1995; Matheson et al. 2000; Maeda et al. 2015; Fremling et al. 2019). However, it is more likely that the late-time emergence of the SN–CSM interaction in these cases is simply due to the decreasing importance of other powering mechanisms (e.g., ⁵⁶Ni/Co decay), and it would not require an increase of the SN–CSM interaction power by a distinct high-density CSM component (Maeda et al. 2015). It is naturally expected if the timing of the strong binary interaction occurred earlier for these systems than for SN 2018ivc; in this case, the CSM created by the first rapid mass-transfer phase must be located far outside.

To our knowledge, this is the second SN IIb for which a clear radio rebrightening has been detected; the first example is the somewhat peculiar SN IIb 2003bg, suggested to be a broad-lined SN IIb (Hamuy et al. 2009) that showed a jump in the radio flux by a factor of 2 or 3 at ∼100 days (Soderberg et al. 2006). We also note that SN IIb 2001ig showed a modulation superimposed on a long-term decline in its radio light curves (Ryder et al. 2004). These SNe IIb may share some common evolutionary path to SN 2018ivc. Possible modulations were also seen in SN IIb 2004C (DeMarchi et al. 2022), though the sparse sampling at such late times makes it hard to ascertain their magnitude and timescale. In any event, SN 2004C did not show a sustained transition from decay to rebrightening. It is nevertheless one of the only two SNe IIb included in the sample of late-time radio monitoring by Stroh et al. (2021), with the other example lacking information on its temporal evolution.

Radio rebrightening is rare for SNe, and only a handful of events are known so far (e.g., Bauer et al. 2008). Of particular interest are a few SNe Ib/c showing a transition to SNe IIn and the late-time radio rebrightening with the timescale covering ∼1 yr to decades (Anderson et al. 2017; Chandra et al. 2020; Thomas et al. 2022); Margutti et al. (2017) showed that about 10% of SNe Ib/c exhibit radio (centimeter) rebrightening ∼1 yr after the explosion with typically an order of magnitude enhancement in the flux. It should be interesting to compare the characteristics of SN 2018ivc with those events, especially SN 2014C as the most well-studied example. Margutti et al. (2017) argued that the inner region is more like a cavity surrounded by an extremely dense CSM; an extremely high mass-loss rate (∼1 M_⊙ yr⁻¹) ∼10–1000 yr before the explosion was followed by a mass loss consistent with a wind from a Wolf–Rayet star (≲10⁻⁵ M_⊙ yr⁻¹ Smith 2017). Margutti et al. (2017) suggested as one possibility the scenario where the common envelope is responsible for the creation of the detached dense CSM, leaving a bare He star. The case for SN 2018ivc has two remarkable differences: (1) the initial, main mass loss (∼1500 yr before the explosion) is less dramatic, and then (2) it keeps a relatively high mass-loss rate toward the explosion. Therefore, it fits better with the scenario proposed here (see also Ouchi & Maeda 2017; Maeda et al. 2023); the initial mass transfer was rapid but did not enter into the common envelope, and the rapid phase was over when the mass ratio was inverted. At that moment, ∼1 M_⊙ of the H-rich envelope was left, which kept filling the RL. An interesting alternative scenario may be a merger of the cores following the common envelope, in which the excited merger product may keep a high mass-loss rate up until the explosion (Nomoto et al. 1995, 1996).

Millimeter observations of SNe are extremely rare; therefore, virtually no samples are available to make any data-based estimate of the frequency of SNe showing rebrightening at millimeter wavelengths. Indeed, a well-defined control sample is still missing even at centimeter wavelengths (Bietenholz et al. 2021), which makes even a qualitative estimate of the frequency of rebrightening events somewhat fraught, with the possible exception of SNe Ib/c (see above; Margutti et al. 2017). Alternatively, we could estimate an expected rate of such events based on the scenario proposed here (i.e., the boundary between the single and binary evolution channels). Maeda et al. (2023) gave such an estimate for the "SN 2018ivc–like" events based on the model sequence of Ouchi & Maeda (2017) as being ∼10% of SNe IIb and ∼3%–6% of stripped-envelope supernovae (SESNe). This estimate indeed indicates that a nonnegligible fraction of SNe Ib/c showing rebrightening (∼10% of SNe Ib/c; see above) may be linked to this binary evolution channel. On the other hand, the same fraction was estimated to be ∼10%–20% of SNe IIL, indicating that the evolutionary channel considered here is not a major channel toward SNe IIL, which might reflect potential "mixed populations" within SNe IIL (Arcavi et al. 2012; but see also Anderson et al. 2014). As compared to SNe IIP, SN 2018ivc–like events are very rare, with an expected rate of only ∼3% of SNe IIP (with large uncertainties), where we assume that the relative proportion of SNe IIP to SNe IIL is ∼7:1 (Li et al. 2011). These crude estimates suggest that the SN 2018ivc–like events, or the events linking the single and binary evolution channels, are intrinsically rare in the population of H-rich SNe. The SNe IIL would make especially interesting targets, as the above estimate suggests that the bulk of SNe IIL are unlikely to have originated in the specific binary evolutionary channel proposed here. A future sample of SNe IIL followed up in the radio, particularly at millimeter wavelengths, may reveal whether such rebrightening is exclusively associated with the SN 2018ivc–like objects.

Another possibility for the detached CSM is an external effect. For the detached CSM shell associated with SN 2014C, Milisavljevic et al. (2015) considered photoionization confinement of the CSM by an external radiation field imposed by neighboring stars in a stellar cluster based on the scenario proposed by Mackey et al. (2014). This scenario, however, is unlikely to work for SN 2018ivc; Bostroem et al. (2020) placed an upper limit of the MS mass of a progenitor star being ≲11 M_⊙ that existed at the SN location in pre-SN Hubble Space Telescope images, even if it was a single star that exsited at the SN location.

We have presented long-term monitoring of the peculiar SN IIL 2018ivc (∼4–1364 days) with ALMA. It started showing rebrightening in its synchrotron flux ∼300–500 days after the explosion. The radial CSM distribution as reconstructed from the long-term millimeter light curves shows a very high mass-loss rate exceeding 10⁻³ M_⊙ yr⁻¹ ∼1500 yr before the explosion, followed by a moderately high mass-loss rate of ≳10⁻⁴ M_⊙ yr⁻¹ in the remaining evolution toward the SN explosion. This behavior is in line with the expectation from the scenario in which SN 2018ivc is indeed more like an "extended SN IIb" in its ejecta properties (i.e., a He star surrounded by ∼1 M_⊙ of the H-rich envelope) produced through the binary with a relatively large initial separation at a boundary between SNe IIb and SNe IIP (Ouchi & Maeda 2017; Maeda et al. 2023).

To our knowledge, it is the first example where rebrightening is detected in the millimeter synchrotron emission from SNe. Further, in the abovementioned scenario, it is the second example where the radio synchrotron rebrightening has been observed for an SN IIb. Indeed, SNe showing radio rebrightening are very rare. The comparison to the properties of the well-observed SN 2014C showing a transition from SN Ib to SN IIn shows that they likely share a similar (binary) evolution channel, albeit with some different details. It is likely that the strong binary interaction started at similar epochs (∼100–1000 yr before the explosion), while the outcome is different; SN 2018ivc is more likely explained by the rapid RL overflow (nonconservative) mass transfer leaving an RL-filling primary star, while SN 2014C entered into the common envelope and stripped away essentially all of the H-rich envelope. Identifying the cause of the difference is beyond the scope of the present work, but it may be due to, for example, the masses of the binary components (e.g., the mass ratio).

The models presented here suggest that ∼1 M_⊙ of the CSM has been swept up by day ∼1500. As this is comparable to the expected envelope mass, the shock will start experiencing rapid deceleration afterward. Therefore, we predict that the synchrotron flux will shortly turn into the decay phase again. Further, if the mass budget for the CSM is ∼5–10 M_⊙, as expected in the binary evolution scenario for SNe IIb (and SESNe), we may assume that the characteristic shock velocity would be ∼5000 km s⁻¹ after a substantial amount of this CSM is swept up. Given that there is additional CSM for another (0.5–1) × 10¹⁷ cm, the shock will eventually break out of this CSM in 3–6 yr, after which the interaction will be essentially over. As such, we plan to continuously monitor SN 2018ivc, which could be done in part of the monitoring observations of NGC 1068. In addition, it could be an interesting target for continuous monitoring with very long baseline interferometry; we expect to start spatially resolving the shock expansion at ≳5 yr for an angular resolution of ∼0001. It could be an interesting target not only for existing facilities but also for future facilities such as the next-generation VLA (Murphy et al. 2018).

The authors thank the referee for constructive comments. This paper makes use of the following ALMA data: ADS/JAO.ALMA #2018.1.01135.S, #2018.1.01193.T, #2018.1.01506.S, #2018.1.01684.S, #2018.A.00038.S., #2019.1.00026.S, and #2021.A.00026.S. Some of the ALMA data were retrieved from the JVO portal (http://jvo.nao.ac.jp/portal/) operated by ADC/NAOJ. 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. K.M. acknowledges support from Japan Society for the Promotion of Science (JSPS) KAKENHI grants JP18H05223, JP20H00174, and JP20H04737. K.M. was supported by the ALMA Japan Research Grant of the NAOJ ALMA Project, NAOJ-ALMA-269. T.M. appreciates the support from NAOJ ALMA Scientific Research grant No. 2021-17A. T.M. is supported by JSPS KAKENHI grant No. JP22K14073. H.K. was funded by Academy of Finland projects 324504 and 328898. M.I. acknowledges support from JSPS KAKENHI grant JP21K03632. The work is partly supported by the JSPS Open Partnership Bilateral Joint Research Projects between Japan and Finland (K.M. and H.K.; JPJSBP120229923).