Discovery of a New Molecular Bubble-Outflow Structure in the Taurus B18 Cloud

Star formation can produce bubbles and outflows, as a result of stellar feedback. Outflows and bubbles inject momentum and energy into the surrounding interstellar medium, and so are related to the overall energy balance of the molecular cloud. Molecular bubbles can be resolved by higher-resolution radio telescopes to quantify the effect of star formation on molecular clouds. We report here the identification of a new molecular bubble with an outflow, and an Herbig Haro object, HH319, located at the bubble center. Multi-wavelength data have been utilized to study its spatial structure, energy injection, and dynamical timescale. This bubble has a kinetic energy of $\rm 5.8 \times 10^{43}$ erg within the smallest radius of a bubble in Taurus, 0.077 pc. The bubble formed $\sim$70,000 years ago. According to the proper motion velocities of protostars from $Gaia$ EDR3, the T Tauri binary stars (FY Tau and FZ Tau) at the southwest edge of the bubble may have produced the outflow-bubble structure. This is an unusual new structure found in low- and intermediate-mass star formation regions. Only a bubble in Orion A, driven by V380 Ori, has a similar structure. The bubble-outflow structure provides additional observational evidence for the theory of stellar wind from T Tauri stars. It enhances our understanding of how stellar feedback acts on molecular clouds.


Introduction
Star formation is related to the structure and evolution of galaxies and giant molecular clouds, the synthesis of chemical elements, and the star formation rate and is intimately connected with a critical astrophysical parameter, the initial mass function (McKee & Ostriker 2007). Almost all young stellar objects (YSOs) have undergone periods of substantial mass loss. YSOs' fast and collimated stellar wind sweeps through the surrounding molecular gas. The molecular gas is expelled from the cavities and expands as irregular lobes and incomplete shells (Bachiller 1996). These become visible as molecular bubbles and molecular outflows. Herbig-Haro (HH) objects, consisting of knots of ionized gas, are the optical manifestations of this mass-loss process McKee & Ostriker 2007).
As a very common phenomenon in both high-mass and lowmass star formation regions, molecular outflows and molecular bubbles are two forms of stellar feedback (e.g., Lada 1985;Arce et al. 2007;Li et al. 2015). Both outflows and bubbles are important in the evolution of the surrounding interstellar medium, injecting momentum and energy into the cloud. Even very low-mass objects can inject energy into the molecular cloud (Phan-Bao et al. 2008). In the Taurus molecular cloud complex, the energy injection rate from the bubbles and outflows can match the turbulence dissipation rate (Li et al. 2015). Molecular bubbles have a greater impact on molecular clouds than outflows. The dynamical timescales of bubbles are longer than those of outflows in Taurus (Li et al. 2015). Thus, the energy injection from bubbles into the molecular clouds could last longer. The ensemble of bubbles in Taurus have about 110 times greater mass and 24 times higher energy than the total of the outflows (Li et al. 2015). Compared to supernova remnants, molecular bubbles could exist for a longer period of time (Arce et al. 2011). Due to their complex morphologies, bubbles are more difficult to identify than relatively clear outflow features (Arce et al. 2011). Therefore, the study of molecular bubbles is both important and challenging.
Bubbles arise from a largely spherical stellar wind. Observationally, bubbles are similar to complete or partial rings. Arce et al. (2011), Li et al. (2015), and Feddersen et al. (2018) have completed a systematic effort to identify bubbles in the Perseus molecular cloud complex, Taurus molecular cloud complex, and Orion A giant molecular cloud. From lowand intermediate-mass star formation regions to massive star formation regions, the minimum mass of a detectable bubble is 2 M e with an energy of 2 × 10 43 erg. The minimum radii of the bubbles are 0.28 pc and 0.050 pc for Taurus and Orion, respectively (Li et al. 2015;Feddersen et al. 2018). The previous bubble identifications in Taurus and Perseus are energetically complete (Liu et al. 2019). Current observations of the bubbles on small spatial scales are still limited by the telescopes' sensitivities and sky coverages. Bubbles from lowmass stars are capable of being resolved by higher-resolution telescopes. Studies of bubbles help us understand the variability of energy injection and spatial structure, the evolution of YSOs, and the evolution of stellar wind. Protostellar outflows are traced by radiative shock waves at near-UV, visible, and infrared wavelengths, or by high-velocity line wings on submillimeter to millimeter wavelength emission lines (Bally 2011). Extensive outflow surveys have been conducted in the Perseus, Ophiuchus, and Taurus molecular clouds (Arce et al. 2010;Nakamura et al. 2011aNakamura et al. , 2011bLi et al. 2015). In the Taurus molecular cloud complex, the minimum mass of outflows is 0.001 M e and the minimum energy is 1 × 10 41 erg (Li et al. 2015).
Increasing evidence suggests that the collimation and morphology of the outflow change with time (e.g., Arce & Sargent 2006;Seale & Looney 2008;Velusamy et al. 2014). At an early evolutionary stage, a young protostar could launch a collimated wind and jetlike outflow (e.g., Arce et al. 2007;Frank et al. 2014;Bally 2016). More evolved pre-mainsequence stars drive wide-angle or spherical winds (e.g., Castor et al. 1975;Bally 2011). Pre-main-sequence stars, such as T Tauri stars and Herbig Ae/Be stars, drive wide-angle winds that coexist with a collimated jetlike wind component (Arce et al. 2011).
At a distance of 140 pc (Torres et al. 2009), the Taurus molecular cloud complex is a widely studied low-mass starforming region. The B18 cloud (also known as TMC-2; Onishi et al. 1996) is located in the southern part of Taurus. The central part of B18 contains dozens of young stars and several HH objects. We have detected a new molecular bubble in the Taurus B18 cloud. A previously identified outflow is located at the center of the bubble (Li et al. 2015). The present observations allow us to derive the properties and potential interplay of the bubble and the outflow. With our analysis of the structure, estimation of the energy injection, and analysis of the formation history, we can understand the remarkable effects of low-mass star formation on a surrounding cloud.
In Section 2, we present the basic information about the IRAM 30 m observations. In Section 3, we discuss the spatial structure, velocity information, previous studies, and physical parameter estimation of the bubble, the outflow, and the HH objects. In Section 4.1, we discuss the possible driving sources and formation history of this structure. In Section 4.2, we compare our source with existing observations and models. In Section 5, we summarize the main conclusions of this paper.

Observations
We performed 12 CO, 13 CO, and C 18 O J = 2 → 1 observations toward the bubble by the IRAM 30 m telescope and 12 CO J = 3 → 2 observations toward the Taurus B18 cloud with the James Clerk Maxwell Telescope (JCMT). We collected multiwavelength data and summarize the images of the bubble-outflow structure in radio, infrared, and optical wavelengths and show the position of this structure in the Taurus molecular cloud in Figure 1. We present the IRAM 30 m 12 CO, 13 CO, and C 18 O J = 2 → 1 observations in Section 2.1 and the JCMT 12 CO J = 3 → 2 observations (Y. Duan et al. 2023, in preparation) in Section 2.2.
2.1. 12 CO J = 2 → 1 from IRAM 30 m Telescope The observations were performed with the IRAM 30 m telescope on January 10, 11, and 12 in 2019 (program ID: 161-18). We observed the 12 CO, 13 CO, and C 18 O J = 2 → 1 emission in a 6′ × 5′ region centered at 4 h 32 m 35 s +24°21′00″ (J2000). The rest-frame frequencies are 230.538000 GHz, 220.398684 GHz, and 219.560358 GHz for 12 CO, 13 CO, and C 18 O J = 2 → 1, respectively. The half-power beamwidths (HPBWs) at the three frequencies are 10°.7, 11°. 2, and 11°.2, respectively. The sensitivities of the antenna temperature for 12 CO, 13 CO, and C 18 O J = 2 → 1 are 0.35 K, 0.21 K, and 0.21 K for a velocity resolution of 0.25 km s −1 , 0.27 km s −1 , and 0.27 km s −1 , respectively. The spectra were treated and baseline-fitted using the GILDAS software. 9 The main-beam efficiency is 0.59 at 230 GHz. The selected OFF position for these observations was set incorrectly and shows some emission. As a consequence, the flux of the 12 CO and 13 CO observations is partly missing. We observed the OFF position independently to correct the data. We refer the reader to Appendix A for details. The shown data include the correction.
2.2. 12 CO J = 3 → 2 from JCMT We utilized a 1.3 deg 2 12 CO J = 3 → 2 mapping of the Taurus B18 cloud from a 14 hr JCMT-HARP observation (Y. Duan et al. 2023, in preparation). The observations were performed in band 3 on 2017 September 6, 11, and 13; 2017 November 14; and 2018 August 10 (program ID: M17BP027 and M18BP072). 12 CO J = 3 → 2 has a rest frequency of 345.79599 GHz and the selected spectral resolution is 61 kHz (0.05 km s −1 ). The telescope's HPBW is 14″ at this frequency. The sensitivity of the antenna temperature is 1.5 K for a velocity resolution of 0.05 km s −1 . The data were processed with the Starlink package (Currie et al. 2014).

Result and Analysis
In the 12 CO J = 3 → 2 map of Taurus B18 observed by JCMT (Y. Duan et al. 2023, in preparation), we discovered a new molecular bubble and acquired IRAM 12 CO, 13 CO, and C 18 O J = 2 → 1 mapping. Combining this mapping with the multiwavelength data, we conduct analysis of the bubble (in Section 3.1), the central outflow (in Section 3.2), the HH object (Haro 6-19 in Haro 1953, now known as HH 319; in Section 3.3), and the physical parameter estimation (in Section 3.4).
A three-color map of the identified bubble is shown in Figure 2, exhibiting the structure of the bubble, outflow, and YSOs at optical, infrared, and radio wavelengths. The three colors represent the integrated intensity maps of the 13 CO (red), 250 μm (green), and Hα (blue) images, which are in panels (j), (g), and (c) of Figure 1, respectively. The Hα image was obtained from the NOIRLab Astro Data Archive (formerly the NOAO Science Archive). 10 In Figure 1(c), there is Hα emission tracing HH 319, marked as a white cross. At the southwest of the bubble, there is strong Hα emission near the T Tauri binary stars (FY Tau and FZ Tau), marked as two white stars. In Figure 2, the distributions of the molecular gas traced by the 13 CO 5-6 km s −1 integrated intensity map and the dust traced by the Herschel 250 μm image are not coexistent. This may be due to the 13 CO and C 18 O images in velocities greater and less than 6 km s −1 showing an east-west symmetrical distribution (see Figure 4 for details). The 13 CO and C 18 O data integrated from 6 to 7 km s −1 should be associated with dust.
For the 5 to 6 km s −1 plots, the molecular gas is symmetrically distributed with respect to the dust emission.

Molecular Bubble
For the bubble survey of molecular clouds, we followed previous studies using similar steps to do the identification (Arce et al. 2010;Li et al. 2015;Feddersen et al. 2018). The steps in the commonly employed empirical procedure used to identify bubbles are summarized by Liu et al. (2019) as follows: (1) Search for bubble-like structures brighter than the surrounding molecular gas.
(2) Draw the position-velocity diagram (hereafter P-V diagram) and search for "V" or "∪" shapes in the figure. An expanding bubble should appear with redshifted or blueshifted gas shown in the P-V diagram. Here, we can see the bubble-like structure in the 12 CO, 13 CO, and C 18 O integrated intensity maps. Due to the influence of the outflow, the P-V diagram is more complicated. We draw the P-V diagram in Figure 3 avoiding the location of the outflow. We cannot distinguish clear "V" or "∪" shapes in the P-V diagram. The 13 CO gas in the bubble interior is redshifted and the bubble is blueshifted. This slight trend tends to form a "V" shape in the P-V diagram. We searched for nearby protostars and discuss them in Section 4.1.2. We checked the Herschel SPIRE map at 250 μm and found that the bubble is traced by the dust continuum (see Figures 1(d) and 2). The bubble formation pushed the dust away and formed a cavity. Using the above four steps, we identified this molecular bubble.
The channel maps of the 12 CO, 13 CO, and C 18 O J = 2 → 1 data are shown in Figure 4 for the bubble and in Figure 5 for  (Goldsmith et al. 2008;Narayanan et al. 2008). (b) Integrated intensity map of B18 cloud in Taurus, as traced by 12 CO J = 3 → 2 from JCMT (unpublished data; Y. Duan et al. 2023, in preparation). (c)-(m) Zoomed-in bubble-outflow images from the gray box region shown in panel (b). The sky coverages in panels (c)-(m) are the same. The images are an Hα image (c) from the NOIRLab Astro Data Archive; Spitzer 3.6 μm (d), 4.5 μm (e), and 8 μm (f) images (Fazio et al. 2004;Werner et al. 2004;IRSA 2022); a Herschel 250 μm image (g) (Eales et al. 2010;Valiante et al. 2016; H-ATLAS Team 2020); a 12 CO J = 3 → 2 image from JCMT (h) (Y. Duan et al. 2023, in preparation); 12 CO (i), 13 CO (j), and C 18 O J = 2 → 1 (k) images from the IRAM 30 m telescope; and 12 CO (l) and 13 CO J = 1 → 0 (m) images from the ∼100 deg 2 FCRAO CO survey in the Taurus molecular cloud (Goldsmith et al. 2008;Narayanan et al. 2008). All the CO maps are integrated in the velocity range of 5-6 km s −1 . The blue circles in the lower left corner of panels (g)-(m) are the beam sizes. HH 319 has a knotty structure. The white cross and squares are HH 319 and other HH objects presented by Magakian et al. (2002) (see Section 3.3 for details). The white stars are the T Tauri binary stars (FY Tau and FZ Tau). the outflow. The high-resolution mapping reveals a clear structure of the outflow located at the center of the bubble. In the 12 CO channel maps, the bubble presents a clear ring shape. To maximize the contrast, we drew the 13 CO and C 18 O maps in the range of 5-6 km s −1 , showing the asymmetrical distribution of the surrounding cloud.
We analyzed the Spitzer images and discovered that the bubble is brighter than its surroundings at 3.6 μm and 4.5 μm and can be seen in absorption at 8 μm in panels (d)-(f) of Figure 1. We show a three-color map of the three bands in Figure 6. This could be a case of "coreshine" emission (Pagani et al. 2010), caused by the scattering of infrared starlight from relatively large (∼0.5-1 μm radius) grains. The existence of the "coreshine" effect signifies an inner dense region lacking any UV or optical illumination.
The east side of this bubble, which is the original unperturbed material, is traced by the Herschel and Spitzer images. The cloud has been locally destroyed by stellar feedback, creating a rupture in the west, which is also seen as a clear deficit in the molecular gas. The stellar wind sweeps up the gas and creates a cavity at multiple wavelengths.

The Central Outflow
The outflow has been studied before. Lichten (1982) first reported the detection of a typically broad, low-intensity 12 CO J = 1 → 0 line wing there. The line wing at 4 h 32 m 45 s +24°25′ 18″ (J2000) has a full width of 30 km s −1 (Lichten 1982). The rms noise level is 0.09 K. The collimated bipolar outflow is considered an explanation for the broad molecular line wings. There is no specific anisotropic spatial distribution of red and blue lobes of an outflow. Lichten (1982) did not confirm the presence of a bipolar outflow. Using 12 CO J = 1 → 0 mapping with a higher resolution than that of Lichten (1982), Li et al. (2015) first identified the existence of a bipolar outflow here, named TMO_16 (SST 043231.7 + 242002). We detected this bipolar outflow with IRAM 12 CO J = 2 → 1 data.
We can see the high-speed gas flow in the P-V diagram of the outflow in Figure 7. In the P-V diagram, the blue lobe shows a high velocity of gas up to 1.9 km s −1 , at 24°21′-24°2 2′. The red lobe is weaker than the blue lobe. At 24°22′, the red lobe shows a faint gas of high velocity, greater than 6.5 km s −1 . The spatial distribution of the red lobe is more extended and weaker than that of the blue lobe. The average spectra in Figure 7 of the two lobes reveal the broad line wings of the outflow.

HH 319 in the Outflow
HH 319 is located at the center of the outflow, providing crucial information about its origin. In the Hα and [S II] images, HH 319 has an elongated, arrowhead-like shape and consists of several knots. Magakian et al. (2002) identified seven knots of HH 319 (named "HH 319A"-"HH 319G") and three adjacent independent HH objects (named "Obj.1," "Obj.2," and "Obj.3"; see their Figures 1 and 2). The coordinates of all the HH objects were given (see Magakian et al. 2002, Table 1). In the entire text, HH 319 is represented as a cross. The three other HH objects are represented as squares.
A weak bubble in the Hα image (see Figure 1 in Magakian et al. 2002) surrounds HH 319, at the same position of the molecular bubble, which is probably due to the shock. Magakian et al. (2002) considered the T Tauri stars (FY Tau and FZ Tau) could be the possible driving sources of HH 319. The position angle of the proper motion of HH 319 is 72°. The angle from the plane of the sky is 53° (Magakian et al. 2002). The proper-motion velocity of HH 319 is 0 05/yr (Luyten 1971), which corresponds to a tangential velocity of 30 km s −1 at the distance of 140 pc of Taurus (Magakian et al. 2002).
Knot A in HH 319 (refer to Figure 2 in Magakian et al. 2002) is the brightest part of HH 319, located at the apex of the shock. Throughout the text, the position of HH 319A represents HH 319. HH 319A is located in the southeast direction. To comply with its direction, the driving source of HH 319 must be located in the northwest, inconsistent with the position of the T Tauri stars. We explain that by giving a possible proper-motion trajectory of the T Tauri stars in Section 4.1.2. Figure 2. Three-color map of the bubble-outflow structure. The red represents the 13 CO integrated intensity. The velocity range is 5-6 km s −1 (range selected to provide the best contrast). The green represents the Herschel 250 μm image (Eales et al. 2010;Valiante et al. 2016; H-ATLAS Team 2020). The blue represents the Hα image. We mark HH 319 and other HH objects as a green cross and red squares, respectively. The green stars represent the pair of T Tauri stars. Figure 3. Left: gray-scale and contour map showing the integrated intensity map of 13 CO J = 2 → 1. The integrated velocity range is 5-6 km s −1 . The contour levels are 60%, 65%, ..., 95% of the peak value. The data are smoothed to 20″. The solid white lines represent the cut for the P-V diagram shown in the right panel. The black dashed line is the circle from the fitted radius (discussed in Section 3.4). The green cross and red squares are HH 319 and other HH objects. The green stars are the T Tauri binary stars (FY Tau and FZ Tau). Right: the 13 CO J = 2 → 1 P-V diagram along the decl. direction of a 36″ wide slice. This interval was chosen to avoid the area of the outflow. The contour levels are 20%, 30%, K, 90% of the peak value.

Physical Parameters
To compare the physical parameters of the bubble with the total Taurus bubble inventory, we follow the method of Li et al. (2015). The physical parameters of the outflow and the bubble are summarized in Table 1.
In Table 1, ΔV has different significations for the outflow and the bubble. The outflow velocity, , is the average observed velocity of the outflow (V obs ) relative to the cloud systemic velocity (V sys ) (Arce et al. 2010;Li et al. 2015). V sys is from the Gaussian fitting of 13 CO, which is 6.2 km s −1 . The observed average velocities of the blue and red lobes are 3.6 km s −1 and 8.4 km s −1 , respectively. For the outflow, ΔV(blue) = 2.6 km s −1 and ΔV(red) = 2.2 km s −1 . For the molecular bubble, ΔV is the expansion velocity of the bubble determined by visual inspection in channel maps (Li et al. 2015).
S in Table 1 is the area of the bubble and of the outflow projected onto the sky plane. For the outflow, the blue lobe is better resolved with the resolution of 11″ of the IRAM 30 m telescope. We estimate the area of the blue lobe as 0 8 × 0 8 rather than 2′ × 2′ from Li et al. (2015). For the red lobe of the outflow, the gas structure is extended. We confirm the 4′ × 4′ size found by Li et al. (2015).
The distance of Taurus is 140 pc (Torres et al. 2009). We measure the maximum length of the outflow, denoted as max l . R bubble is the radius of the bubble, which is fitted by leastsquares fitting. The 12 CO J = 2 → 1 image in the 5.69 km s −1 velocity channel is cut into a 5 × 5 array of equal-sized tiles. The brightest point of each tile is selected and points close to the bubble periphery are manually selected and used to fit the circle. The fitted center of the circle is 4 h 32 m 35 6 +24°21′46″ (J2000) with a radius of 0.077 pc.
We use the 13 CO and 12 CO data to calculate the column density of the bubble and the outflow, respectively. For the outflow, we use 12 CO J = 2 → 1 data integrated over 2.0-4.5 km s −1 for the blue lobe and 7.5-9.0 km s −1 for the red lobe. For the bubble, we use a 13 CO J = 2 → 1 map integrated over 4.7-8.0 km s −1 . The detailed derivations of column density (N tot ), mass (M gas ), momentum (P), energy (E), and energy injection rate (  E) for the outflow and the bubble are given in Appendix B. Referring to the Taurus outflow and bubble survey (Li et al. 2015), the excitation temperature, T ex , is assumed to be 25 K. We assume that 12 CO is optically thick and 13 CO is optically thin in the calculations of the column densities. We summarize the results in Table 1.
Among the 37 bubbles in Taurus studied by Li et al. (2015), the minimum and median radii are 0.28 pc and 0.70 pc, respectively. The IRAM 30 m telescope at CO J = 2 → 1 has one-quarter the beam size of the FCRAO 13.7 m telescope at CO J = 1 → 0. Our CO J = 2 → 1 molecular bubble has a radius of 0.077 pc, which is smaller than any in the Taurus bubble survey by FCRAO. The minimum and median masses of the 37 bubbles in Taurus studied by Li et al. (2015) are 2 M e and 12 M e , respectively. The minimum and median momenta are 2 M e km s −1 and 21 M e km s −1 , respectively. The minimum and median energy are 2 × 10 43 erg and 41 × 10 43 erg, respectively. The minimum and median energy injection rates are 0.20 × 10 31 erg s −1 and 3.1 × 10 31 erg s −1 , respectively. Our molecular bubble has a mass and momentum of 5.2 M e and 5.5 M e km s −1 , respectively. The bubble injects 5.8 × 10 43 erg and the energy injection rate is 2.6 × 10 31 erg s −1 . The mass, momentum, energy, and energy injection rate of this molecular bubble are all lower than the median value of the bubble survey and slightly higher than the lowest value of the survey.
The physical parameters of this outflow (named "TMO_16") have been published by Li et al. (2015). With a higherresolution mapping, we obtain a new estimate of the physical parameters for the outflow. The 56 outflows in Taurus identified by Li et al. (2015) have minimum and median lengths of 0.1 pc and 0.3 pc, respectively. The maximum length of our outflow is 0.046 pc and 0.23 pc for the blue lobe and red lobe, respectively. The red lobe is close to the median length, and the blue lobe is smaller than the minimum value. It is reasonable that a higher resolution could reveal the subtle structure of the same source. The minimum mass and median mass for the 55 outflows in Taurus are 0.001 M e and 0.034 M e , respectively. The minimum and median momenta are 0.003 M e km s −1 and 0.098 M e km s −1 , respectively. The minimum and median energies are 0.01 × 10 43 erg and 0.26 × 10 43 erg, respectively. The minimum and median energy injection rates are 0.004 × 10 31 erg s −1 and 0.053 × 10 31 erg s −1 , respectively (Li et al. 2015). The mass, momentum, and energy of the red lobe for this outflow are close to the median values of the outflow survey, and those of the blue lobe are slightly higher than the lowest values of the outflow survey. The energy injection rates of the blue lobe and red lobe are close to the median level.
To estimate the dynamical timescale, we expect that the protostellar outflow and bubble formed in a relatively quiescent cloud, and that they probably formed at different times. The . Channel maps of 12 CO, 13 CO, and C 18 O J = 2 → 1 from IRAM observations. The central velocity of each channel is indicated in the lower left corner. The white box in the leftmost 12 CO and 13 CO images delineates the zoomed-in images displayed in Figure 5. The green cross and red squares are HH 319 and other HH objects presented by Magakian et al. (2002). The green stars are the T Tauri binary stars. dynamical timescale (t dyn ) can be estimated by In Magakian et al. (2002), the angle of the HH 319 flow with respect to the plane of the sky is 53°. The angle between HH 319 and the line of sight should be 37°. We assume that the inclination angle of the outflow with respect to the line-of-sight direction coincides with the angle of HH 319. On the basis of the inclination angle, 37°, the P corr , E corr , t dyn,corr , and  E corr of the outflow are corrected to 1.3 times, 1.6 times, 1.3 times, and 1.2 times the uncorrected values, respectively. A molecular bubble should be a 3D sphere in space and its dynamical timescale requires no correction. We do not consider the effect of the sky projection for the bubble here. In Table 1, only the momentum, energy, dynamical timescale, and energy injection rate of the outflow (P corr , E corr , t dyn,corr , and  E corr ) are modified by the 37°inclination angle of the line-of-sight direction. The rest of the physical parameters in Table 1 are uncorrected by default, including max l and ΔV. The dynamical age of the bubble is 0.71 × 10 5 yr, lower than the minimum dynamical timescale of the 37 bubbles in Taurus (Li et al. 2015), ∼1 × 10 5 yr. The dynamical age of the blue lobe of the outflow is between 0.18 × 10 5 yr and 0.23 × 10 5 yr, which is lower than the minimum dynamical timescale of the Taurus outflow survey (Li et al. 2015), ∼0.3 × 10 5 yr. For the red lobe of the outflow, the dynamical timescale is between 1.0 × 10 5 yr and 1.4 × 10 5 yr, which is close to the median value in Taurus, ∼1.5 × 10 5 yr. This is reasonable as compared to the estimated dynamical age of the Taurus molecular cloud, ∼10 7 yr (Vázquez-Semadeni et al. 2018), and the dynamical ages of Taurus stellar feedback, ∼0.3 × 10 4 to 10 6 yr (Li et al. 2015).
The two lobes of the outflow have asymmetrical dimensions and inconsistent intensities. The morphology, spatial scale, and estimated dynamical timescale of the blue lobe are significantly different from those of the red lobe of the outflow. The blue lobe has a clearer spatial distribution in the integrated intensity map in Figure 7. The red lobe is more extended. In the P-V diagram, the high-velocity gas of the blue lobe deviates from the system velocity more than that of the red lobe. Here we only discuss the formation history between the clearer blue lobe and the bubble.

Formation History
The spatial distribution of the bubble and the outflow determines their past relationship. It requires us to find the driving source of the bubble-outflow structure. At the center of the bubble, there are no identified YSOs. 11 To find the driving source, we discuss two possibilities: (1) There may be a 13 CO Figure 5. Zoomed-in 12 CO and 13 CO images from the white box region shown in Figure 4. The markers are the same as those in Figure 4. core in the center of the bubble, which is not visible in the C 18 O integrated intensity map. A newborn protostar may be embedded in the core at the bubble's center. The previous YSO survey did not identify one. We analyze the mass, extinction, energy, gravitational energy, and kinetic energy of the core in Section 4.1.1.
(2) The known YSOs have moved away out of the bubble. In Section 4.1.2, we discuss the proper motions of the adjacent YSOs and study their moving trajectories in the estimated dynamical timescales of the bubble and outflow.

Physical Conditions of the Core
In the bubble center, IRAM 13 CO J = 2 → 1 data show a dense gas distribution. We check the C 18 O spectrum of this core. This core has C 18 O emission, which is shown in panel (B) in the last row of Figure 9 (bottom) in Appendix A. Point B in Figure 9 represents the position of the outflow and the core. Though the integrated intensity map (in the range 4-8 km s −1 ) of C 18 O does not trace the core revealed by 13 CO, the core can be faintly seen in the C 18 O channel map at 6.16 km s −1 (Figure 4). The large difference between the two isotopologues in revealing the core might be attributed to the outflow wings contributing to the 13 CO emission, which also explain the larger, non-Gaussian line width of 13 CO toward the core (spectrum B in Figure 9). Here we analyze the physical parameters of the 13 CO core and discuss the possibility of protostar formation. The blue and red contours show the integrated intensity map of 12 CO J = 2 → 1, over 2-4.5 km s −1 for the blue lobe and 7.5-9.0 km s −1 for the red lobe. The contour levels are 20%, 40%, 60%, and 80% of the peak value for the blue lobe, and 60%, 70%, 80%, and 90% of the peak value for the red lobe. The white solid line represents a cut for the P-V diagram. The data have been smoothed to 20″. The green stars, green cross, and red squares represent the T Tauri binary stars, HH 319, and nearby HH objects, respectively. Upper right panel: P-V diagram of 12 CO J = 2 → 1 along the line shown above at a position angle of 70°. The contour levels are 10%, 20%, K, 90% of the peak value. Lower left panel: the average lines of the blue lobe. The solid and dotted lines represent 12 CO and 13 CO J = 2 → 1. Lower right panel: the average 12 CO and 13 CO J = 2 → 1 spectra of the red lobe. Normally cloud cores are prolate ellipsoids (Myers et al. 1991;Li 2002;Li et al. 2013). We fit the cloud core with an ellipse. The fitting result is expressed as a core = 0.012 pc (2.6 × 10 3 au) for the semimajor axis and b core = 0.0046 pc (0.94 × 10 3 au) for the semiminor axis. It is smaller than the typical molecular core dimension (∼0.1 pc; Myers & Goodman 1988), and the Jeans length (λ J ∼ 0.83 pc) in the Taurus molecular cloud (Li et al. 2015).
We estimate the gas mass and extinction of the core. We use 13 CO J = 2 → 1 integrated intensity data within the velocity interval of 3.0 to 9.0 km s −1 to estimate the column density of the core. All calculations are based on the main-beam temperature scale. The derivation of the total column density for the core is shown in Equations (B12)-(B16) in Appendix B. The column density for the rotational upper level of the transition J = 2 → 1, N u ( 13 CO), is 2.7 × 10 15 cm −2 for the core. We assume that 12 CO is optically thick and 13 CO is optically thin, and the optical depth of 13 CO can be obtained from Li et al. 2015). The calculated optical depth of the 13 CO J = 2 → 1 emission is 0.44. The correction factor for the column density due to the optical depth, f τ , is 1.24. The excitation temperature, T ex , is assumed to be 25 K (Li et al. 2015). We estimate the fraction of 13 CO in the upper level, f up , to be 0.28. We assume the abundance ratio of 12 CO to 13 CO is 65 (Langer & Penzias 1993 Pineda et al. 2010). Based on our estimate, the maximum visual extinction of the core is 8.1 mag. This result is lower than the 16 mag extinction of the nearby core (core 11 in Pineda et al. 2010). Based on the ellipsoid shape of the core, we write the mass of the core as The mass of the core is 0.029 M e .
We estimate the kinetic energy and the gravitational energy, referring to previous discussion of the virial theorem (e.g., McKee & Ostriker 2007;Li et al. 2013). For an axisymmetric ellipsoid core with concentric density profiles, the gravitational binding energy can be expressed as where α = (1 − a/3)/(1 − 2a/5) for a power-law density profile ρ ∝ r − a . We assume a = 1.6, as in an appropriate isothermal cloud in equilibrium (Bonnor 1956). The geometry factor β = arcsine/e is determined by eccentricity e f 1 2 ( ) = -. The intrinsic axis ratio, f, should be less than the observed f obs under the projection of the sky plane. For our core, f obs is 0.37. We choose the prolate ellipsoid (Fall & Frenk 1983) and find f f f 2 0.5, 0.5, 0.5, 1.5, 1, 1 , 11  is the Appell hypergeometric function of the first kind (Li et al. 2013). We estimate that f = 0.26 and β = 1.4. The radius R in Equation (10) is the geometric average of the semimajor and semiminor axes of the ellipse, 0.0075 pc. The gravitational constant G is 6.67 × 10 −11 N m 2 kg −2 . We estimate the gravitational binding energy of this core to be −9.9 × 10 39 erg. The turbulent kinetic energy of the core is where the velocity dispersion σ is 0.57 × 10 3 m s −1 from the Gaussian fitting of the 13 CO line. We estimate a value of 2.8 × 10 41 erg for this core.
The ratio of the total kinetic energy to the gravitational binding energy is defined as the virial ratio For r vir = 1, the uniform, unmagnetized gas sphere is in virial balance ( 2 | | =  ). We estimate the virial ratio of this core to be r vir = 57. Although we do not have magnetic field data, it appears that this core is definitely not gravitationally bound.
The measured non-gravitationally bound condition of the 13 CO core could have two possible causes.
(1) The bubble and the outflow are from the core. A newborn protostar is buried deep in the core, which generated the bubble and the outflow. This process released enough kinetic energy into its surroundings to cause the core to be in a non-gravitationally bound state.
(2) The outflow and the bubble do not originate from the core, which may be in the process of being destroyed or dismantled. No protostar existed in the core to drive the bubble and the outflow. Both of the above scenarios are possible based on the estimation of the physical parameters of the core.

Nearby Known YSOs
Because of the lack of direct evidence, we can only speculate that a YSO drove the expansion of the bubble and then moved out of the bubble. To explore the origin of this structure and identify the driving source, we check the up-to-date Spitzer YSO catalog (Rebull et al. 2010). There are four possible YSOs: the common names are FY Tau, FZ Tau, CFHT-5, and MHO-8. FY Tau and FZ Tau are a pair of binary stars (Kenyon & Hartmann 1995;Kraus & Hillenbrand 2009;Akeson et al. 2019). We summarize the main information of the four YSOs in Table 2. The common names, coordinates, and YSO classifications are adopted from Rebull et al. (2010). The projected distances to the bubble, d bubble , are the distances between the YSOs and the fitted center of the bubble (4 h 32 m 35 6 +24°21′46″).
For these four nearby possible driving sources, we search for information from Gaia Early Data Release 3 (EDR3; Gaia Collaboration et al. 2016Collaboration et al. , 2020Collaboration et al. , 2021 archive data. The proper-motion velocity (μ α , μ δ ) and parallax of the four stars are taken from Gaia EDR3 (Gaia Collaboration et al. 2016, 2020. We calculate the distance to the Earth (d Earth in Table 2) from the parallax.
We examine the nearby YSOs that have a systematic bulk motion toward the southeast direction. Luhman (2018) used Gaia DR2 to study the systematic proper motions of stars for the Taurus molecular cloud and divided Taurus into nine regions. For one of them, dark clouds L1536, L1529 (i.e., the B18 cloud), and L1524 covered a submap of 3.5°× 4°.25 (see Figure 6 in Luhman 2018). In Figure 6 of Luhman (2018), there is a red population and a blue population distinguished from the parallax. Based on the spatial distribution, the four possible driving sources of the bubble in the B18 cloud here all belong to the red population. The red population in Figure 6 in Luhman (2018) has a coverage of about 2°× 4°. 25 and has a bulk motion of μ α = 7.2 mas yr −1 , μ δ = −24.1 mas yr −1 , Δμ α = −1.0 mas yr −1 , and Δμ δ = 1.5 mas yr −1 (see Table 2 in Luhman 2018). Considering that the error of the Gaia EDR3 data for the four YSOs is smaller than that of Gaia DR2, we choose to calculate the bulk-motion velocity of the YSOs using Gaia EDR3 (Gaia Collaboration et al. 2016, 2020. Acquiring the bulk motions for the four possible driving sources first requires defining the region in which to select the surrounding stars. We adopt the division of this region from Luhman (2018) (red population in Figure 6) as the scope. We set the angular radius of 4°.25 as a selection range for stars in the sky plane and line-of-sight direction, which ensures the full coverage of all Gaia sources in the B18 cloud. For each of the possible driving sources (R.A. 0 , decl. 0 , parallax 0 ), we use Gaia EDR3 to select the stars around it. In the sky projection plane, we select stars in the range In the line-of-sight dimension, star selection should convert 4°. 25 to parsecs and calculate the interval of the line-of-sight distance. This distance interval is converted into a parallax interval. Only stars with parallax in the range The Gaia data indicate that YSOs have systematic motions relative to other stars in the sky projection plane. The four possible driving sources of the bubble are YSOs. Thus, the four possible driving sources have systematic motions relative to stars other than YSOs. To deduce the relative motions of the four possible sources relative to the cloud, we take the median proper motion of stars excluding all YSOs as a reference. We use the following steps for the calculation.
(1) For the four possible driving sources, we select all stars around them within 4°. 25 in Gaia EDR3.
(2) All YSOs are eliminated from the selected stars. The catalog of YSOs is from the latest Spitzer YSO survey (Rebull et al. 2010).
(3) We calculate the median proper motion of all selected stars excluding YSOs to be μ α,bulk_motion , μ δ,bulk_motion ( Table 2). (4) We calculate the corrected proper motion of the driving sources relative to μ α,bulk_motion , μ δ,bulk_motion to be μ α,corr , μ δ,corr ( Table 2). (5) When we select stars in Gaia EDR3, we constrain their parallaxes using Equation (16). Among the four possible driving sources, CFHT-5 is 97 pc away from us, and the other three are 129 pc away from us. The selected stars around them and their bulk motions are different. We conduct the above process four times to obtain the corrected proper motions.
With the dynamical timescale of the bubble and the corrected bulk-motion velocities, we can obtain the motion trajectory, shown in Figure 8. We explore using a selection range that is larger or smaller than 4°.25. The proper-motion trajectories of FY Tau and FZ Tau are shown as movement from the north to the south. Compared to the other two YSOs, the T Tauri binary stars, FY Tau and FZ Tau, are the most probable ones to have originated from the bubble center to the present observed position. They are more likely to explain the bubble-outflow formation than the other two YSOs in many attempts. The distance of CFHT-5 in the line-of-sight direction is 97 pc. There are relatively fewer stars around CFHT-5. Regardless of whether the selection range is greater or less than 4°. 25, CFHT-5 has a different bulk motion compared with the other three possible sources. In the case of MHO-8, many attempts have shown that it is difficult to get a trajectory from the interior of the bubble to where it is now.
T Tauri stars drive wide-angle winds, which in some cases coexist with a collimated jetlike wind component (Arce et al. 2011). They cannot move too far from their birthplace during their lifetime (about 10 6 yr; Gomez et al. 1993). Li et al. (2015) and Magakian et al. (2002) identified the pair of T Tauri stars as the driving source of this bipolar outflow and HH 319. In conclusion, there is a possibility the T Tauri binary stars, FY Tau and FZ Tau, are the driving source of the bubble-outflow structure. We have developed the following possible explanation. The T Tauri binaries (FY Tau and FZ Tau) formed this bubble and outflow. The formation of the bubble preceded the formation of the outflow. Otherwise, the process of the bubble expanding would destroy the structure of the outflow. We reconstruct the trajectories of their movements over the dynamical ages, shown in Figure 8.
Observations have demonstrated that a subset of HH outflows and jets exhibit C-symmetric bending when the source is moving through the medium Bally 2016). In our concept, FY Tau and FZ Tau moved from the bubble interior to their present positions. C-symmetric bending is not present in any optical images (e.g., the Hα image in Figure 1). We consider two possible reasons for its invisibility.
(1) The intensity of the outflow is changing. It may be visible initially and become weak or invisible with the passage of time. The HH object only traces gas that passed through a shock recently. It would become invisible after the post-shock gas cools Bally et al. 2012).
(2) HH flows and objects are formed by the interaction of the YSOs with the surrounding medium. It is possible that the surrounding gas is not dense enough to produce clearly visible structures. If future data can reveal this C-symmetric bending structure, the associated source could be demonstrated to be the driving source of the bubble-outflow structure. At present, FY Tau and FZ Tau are the most probable sources to generate this bubble and outflow.

The Outflow
Both bubble and outflow originate from a period of mass loss during star formation. Generally, the collimation of the stellar wind decreases with the star's evolution (e.g., Arce & Sargent 2006;Seale & Looney 2008;Velusamy et al. 2014). Recently, several models of stellar winds have been developed to explain the different types of morphology, velocity, collimation, and momentum injection (Cabrit et al. 1997). Arce et al. (2007) summarized four types of outflow models that are widely accepted: (1) wind-driven shells (e.g., Li & Shu 1996;Lee et al. 2001), (2) jet-driven bow shocks (e.g., Chernin & Masson 1995;Cliffe et al. 1996;Hatchell et al. 1998;Lee et al. 2001), (3) jet-driven turbulent flows (e.g., Canto & Raga 1991;Chernin & Masson 1995;Bence et al. 1996), and (4) circulation flows (e.g., Fiege & Henriksen 1996;Lery et al. 1999). These models are used to explain the observed spatial structures and velocity fields of outflows. Some outflows are more complex and cannot be explained by a single model. For example, a wide-angle wind and a collimating episodic wind could cooperate to produce the observed red lobe of the HH 46/47 outflow (Arce et al. 2013). Lee et al. (2002) observed some CO outflows with multipolar structures and proposed that a model combining both jet-driven and wind-driven components may best match the observation. There are some other observations. T Tauri stars traced by forbidden emission line profiles show two-velocity components: a high-velocity component that is argued to arise from a jet and a low-velocity component that might result from a disk wind (Kwan & Tademaru 1995;Li & Shu 1996).
The present outflow was identified as a bipolar outflow in a previous study (Li et al. 2015). The two lobes of the outflow traced by 12 CO have different morphologies. This is reasonable and common, similar to the case of the HH 46/47 outflow (Arce et al. 2013) and the GK/GI Tau outflow (Arce & Sargent 2006). The blue lobe of the present outflow is similar to the morphology described by the jet-driven bow shock model (e.g., Chernin & Masson 1995;Cliffe et al. 1996;Hatchell et al. 1998;Lee et al. 2001;Arce et al. 2007). Limited by the low velocity resolution, we are unable to make any more determinations for the blue lobe and the red lobe of this outflow.

The Bubble
Molecular bubbles originate from spherical or very wide angle winds and are powered by young stars in molecular  . The color background is the integrated intensity map of 12 CO J = 3 → 2 from 5 km s −1 to 6 km s −1 (Y. Duan et al. 2023, in preparation). The contours are the Herschel SPIRE 250 μm data. The contour levels are 50%, 55%, K, 95% of the peak value. The white and cyan stars represent the positions of the YSOs at present and 7.1 × 10 4 yr ago, respectively. The arrows show the motion paths of the YSOs. The red crosses represent the upper limit and lower limit of the corrected proper motions for the YSOs. The white cross and squares are HH 319 and other HH objects. The gray square represents the zoomed-in region shown in the lower panel. Lower panel: a zoomed-in view of the above panel, showing the proper-motion paths of the T Tauri binary (FY Tau and FZ Tau). clouds. Norman & Silk (1980) built a dynamical model of molecular clouds that is dominated by collisions between bubbles from the interaction of T Tauri stars with the surrounding quiescent gas. In subsequent systematic surveys of molecular bubbles in Perseus, Taurus, and Orion A giant molecular clouds (Arce et al. 2011;Li et al. 2015;Feddersen et al. 2018), powerful spherical winds from bubbles have been proved as one of the main energy origins of molecular clouds for maintaining turbulence.
In the recent bubble surveys for the three molecular clouds (Arce et al. 2011;Li et al. 2015;Feddersen et al. 2018), we only find one example of a bubble coexisting with an outflow and HH flows. In the Orion A giant molecular cloud, Feddersen et al. (2018) identified two nested expanding shells, Shell 39 and Shell 40. Shell 40 is nested inside Shell 39. A multiple stellar system, V380 Ori, was proven to drive the two molecular shells and an identified molecule outflow (Morgan et al. 1991;Feddersen et al. 2018). V380 Ori (IRAS 05339-0644) is a Herbig Ae/Be intermediate-mass emission line star, one of the youngest Herbig Ae/Be stars (Herbig 1960;Allen & Davis 2008). In contrast to our molecular bubble in Taurus, there is also a far-infrared dark cavity excavated by protostellar jets from the V380 Ori system and located 0.1 pc away from molecular Shell 40 (Stanke et al. 2010). The group associated with V380 Ori drives multiple HH flows. It can drive highvelocity large-scale HH flows, like the 5.3 pc long HH 222/ 1041 flow (Reipurth et al. 2013). Moreover, it can be the driving source for small-scale HH flows, like HH 1031/130 and HH 35, which may represent more recent dynamical interactions. Both types of flows play a role in shaping the shells (Feddersen et al. 2018). An accretion-driven outburst of V380 Ori could explain the observed large-scale HH flows, increased mass-loss rates, and spherical winds (Feddersen et al. 2018). In our case, there is no sign of such an outburst for the T Tauri binary stars, FY Tau and FZ Tau.
Herbig Ae/Be stars and T Tauri stars have different masses and internal structures. Unlike the low mass (1.5 M e ) of T Tauri stars, the masses of Herbig Ae/Be stars are in the range of 1.5-10 M e (Pinzón et al. 2021). Herbig Ae/Be stars are mainly radiative and T Tauri stars are mostly fully convective (Iben 1965). However, similar observational features are found around both Herbig Ae/Be stars and T Tauri stars, such as coexisting outflows, bubbles, HH flows, and far-infrared dark cavities. FY Tau, FZ Tau, and V380 Ori demonstrate that lowand intermediate-mass stars can significantly impact the interstellar medium's physical and chemical properties.
As the possible driving source, FY Tau and FZ Tau are likely to be the first T Tauri binary found to drive a molecular outflow and a molecular bubble together. Our observation of this structure again demonstrates the potential effect of bubbles and outflows for dark clouds. Higher-resolution observations will be required to ascertain the origin of this bubble-outflow structure.

Conclusions
In this paper, we have studied a bubble-outflow structure using observations of multiple transitions of CO together with existing infrared data. We analyze the energy injection, dynamical timescale, possible driving source, and formation history of the bubble and the outflow, and existing observations and models. We draw the following conclusions: 1. We have discovered a new molecular bubble. A known bipolar outflow (Li et al. 2015) is located at the interior of the bubble. The blue lobe of the outflow is coincident with the bubble center. HH 319 is in the blue lobe. 2. The molecular bubble has a mass, momentum, and kinetic energy of 5.2 M e , 5.5 M e km s −1 , and 5.8 × 10 43 erg, respectively, within a radius of 0.077 pc. The energy injection rate is 2.6 × 10 31 erg s −1 . We estimate that the bubble formed 7.1 × 10 4 yr ago. This bubble has the smallest radius and the shortest kinematic timescale among the known 37 bubbles in Taurus (Li et al. 2015). 3. The bipolar outflow has a total mass of 0.039 M e .
Considering a projection angle of 37°with respect to the line-of-sight direction, the outflow has a total angular momentum between 0.087 M e km s −1 and 0.11 M e km s −1 . The energy is between 1.9 × 10 42 erg and 3.0 × 10 42 erg. The energy injection rate is between 0.98 × 10 30 erg s −1 and 1.2 × 10 30 erg s −1 . 4. According to our estimate, the core at the center of the bubble is non-self-gravitationally bound. Its kinetic energy to gravitational binding energy ratio, r 2 vir =  | |  , is 57. 5. The T Tauri binary stars located southwest of the bubble may be the driving sources. After subtracting their bulk motion, we find the T Tauri binary stars (FY Tau and FZ Tau) are the most likely ones to have moved from the bubble center to their present position. Their corrected proper-motion velocities are μ α,corr = −0.71 ± 0.07 mas yr −1 , μ δ,corr = −1.11 ± 0.04 mas yr −1 for FY Tau and μ α,corr = −0.24 ± 0.08 mas yr −1 , μ δ,corr = −1.15 ± 0.03 mas yr −1 for FZ Tau.
We compare the existing systematic bubble surveys (Arce et al. 2011;Li et al. 2015;Feddersen et al. 2018). This is the first example of bubble-outflow coexistence in a low-mass star formation region. In Orion A, there is an example of a colocated bubble and outflow. Our discovery demonstrates the potential of bubbles and outflows in energy injection and modification of structures within dark clouds. It demonstrates the ability of the stellar wind from T Tauri binary stars to make a dramatic impact on the surrounding environment. Further observations are needed to reveal the origin of the bubble and the outflow. The LTE partition function (for hB/kT ? 1) is Q(T) = kT/hB, where B is the rotation constant of CO, B = 57,635.968 MHz 12 (Li et al. 2015;Xie et al. 2021 ] is assumed to be 10 4 , the mean molecular weight μ g = 2.72, and m(H) = 1.67 × 10 −24 g. The momentum and energy of the outflow can be calculated from where T( 12 CO) and T( 13 CO) are the brightness temperatures of 12 CO and 13 CO, respectively. Assuming that 12 CO is optically thick and 13 CO is optically thin, the optical depth of 13 CO can be obtained from The bubble's gas mass (M gas ), momentum (P), energy (E), and energy injection rate (  E) are calculated by the same equations for the outflow.