11institutetext: Departamento de Física, Escuela Politécnica Superior de Jaén,
Universidad de Jaén, Campus Las Lagunillas s/n, A3-420, 23071 Jaén, Spain
11email<EMAIL_ADDRESS>22institutetext: Departamento de Ingeniería Mecánica
y Minera, Escuela Politécnica Superior de Jaén, Universidad de Jaén, Campus
Las Lagunillas s/n, A3-008, 23071 Jaén, Spain
22email<EMAIL_ADDRESS>33institutetext: Departamento de Ciencias Integradas,
Centro de Estudios Avanzados en Física, Matemática y Computación, Universidad
de Huelva, 21071 Huelva, Spain
33email<EMAIL_ADDRESS>44institutetext: Grupo de Investigación
FQM-322, Universidad de Jaén, Campus Las Lagunillas s/n, A3-065, 23071 Jaén,
Spain

# A runaway T-Tauri star leaving an extended trail

Josep Martí 1144 Pedro L. Luque-Escamilla 2244 Estrella Sánchez-Ayaso 3344

(Received September XX, 2022; accepted XXXXXX XX, XXXX)

###### Abstract

Aims. We address the problem of young stellar objects that are found too far
away from possible star formation sites. Different mechanisms have been
proposed before to explain this unexpected circumstance. The idea of high-
velocity protostars is one of these mechanisms, although observational support
is not always easy to obtain. We aim to shed light on this issue after the
serendipitous discovery of a related stellar system.

Methods. Following the inspection of archival infrared data, a peculiar
anonymous star was found that apparently heads a long tail that resembles a
wake-like feature. We conducted a multiwavelength analysis including
photometry, astrometry, and spectroscopy. Together with theoretical physical
considerations, this approach provided a reasonable knowledge of the stellar
age and kinematic properties, together with compelling indications that the
extended feature is indeed the signature of a high-velocity, or runaway,
newborn star.

Results. Our main result is the discovery of a low-mass young stellar object
that fits the concept of a runaway T-Tauri star that was hypothesized several
decades ago. In this peculiar star, nicknamed UJT-1, the interaction of the
stellar wind with the surrounding medium becomes extreme. Under reasonable
assumptions, this unusual degree of interaction has the potential to encode
the mass-loss history of the star on timescales of several $\sim 10^{5}$
years.

###### Key Words.:

Stars: formation – Stars: winds, outflows – ISM: jets and outflows – Shock
waves – Stars: variables: T Tauri, Herbig Ae/Be

##  1 Introduction

Star formation usually occurs in the deep cores of molecular clouds (McKee &
Ostriker, 2007). However, sometimes infant stars appear to be isolated, and
their space velocity is too low (few km s-1) for them to have reached their
current position from any plausible cradle site within their young age (a few
million years; Neuhäuser (1997)). While dispersion of the parental molecular
cloud might eventually overcome this issue (Hoff et al., 1998), the concept of
runaway T-Tauri stars (TTS), also known as RATTS, was offered as a competing
alternative scenario (Sterzik & Durisen, 1995), although only a few candidates
are reported decades later (Neuhauser et al., 1997; Neuhäuser et al., 1996).
Moreover, none of these candidates exhibits unambiguous signatures of a fast-
moving object. In this paper, we present the discovery of a high-velocity TTS
escaping from its natal molecular cloud and leaving behind an extremely long
wake that we estimate to be at least $\sim 10$ pc. This elongated feature
rivals the length of similar features in the Galaxy that are associated with a
single star (Martin et al., 2007).

This finding not only revives the RATTS paradigm, but might also enable us to
recover the mass-loss history of a protostellar object back to several hundred
thousand years. In addition, the long tail behind this TTS also provides a
unique test bench for studying the interplay between turbulence and
instabilities at very large Galactic scales. The present work, dealing with an
obscured low-mass star, complements other recent studies on runaway stars
using Gaia data that were mostly focused on early-type luminous stars (Hattori
et al., 2019; Neuhäuser et al., 2020) or giant unobscured stars (Li et al.,
2021). In these evolved contexts, dynamical ejection in multiple systems and
supernova explosions in close binaries currently appear as the most likely
scenarios.

The paper is organized as follows. After describing the discovery
circumstances and early observational work, we present the evidence supporting
the RATTS nature of the target star. The spectral energy distribution (SED) of
the star itself and of its bow-shock tail is addressed in detail. Next, in the
discussion section, we devote most of our attention to the mechanical
scenario. We scale the properties of the stellar wake, determine the type of
instabilities that eventually develop in its flow, and describe stellar wind
evolution. Anticipating our main conclusion, we appear to have found a RATTS
that looses mass in time due to variable stellar winds and strong interaction
with its surrounding interstellar medium (ISM). Finally, four appendices with
equation formalisms of cooling timescales, hydrodynamical phenomena, and
speculation about the origin of the star are also included.

Figure 1:  View of UJT-1 and its environments. This image uses the WISE All-
Sky data release at 3.4, 4.6, and 12 $\mu$m, coded as blue, green, and red
layers, respectively. The dashed circle shows the position of the star, and
the arrow vector represents the proper motion direction corrected for Galactic
rotation. Axes are labeled in equatorial coordinates; north is rotated to the
right and east is up. The stellar bow shock wake is clearly visible extending
northward. The inset shows a zoomed view of the stand-off region of the bow
shock.  Figure 2:  Least-squares fit to astrometric observations. Data from
Table 1 are plotted both in right ascension (top) and declination (bottom) as
a function of time. The two celestial coordinates are expressed in the ICRS
system. The dashed lines represent linear least-squares fits yielding the
estimated proper motions. Error bars correspond to one standard deviation.
Table 1: Log of astrometric observations

Source of image | Epoch | Right Ascensiona | Declinationa
---|---|---|---
| (year) |  $308^{\circ}\leavevmode\nobreak\ 07^{\prime}$ $+$ $(^{\prime\prime})$ |  $43^{\circ}\leavevmode\nobreak\ 57^{\prime}$ $+$ $(^{\prime\prime})$
1st Digitized Sky Survey | 1953.449 | $36.443\pm 0.100$ | $42.317\pm 0.105$
2nd Digitized Sky Survey | 1992.468 | $36.453\pm 0.071$ | $42.077\pm 0.067$
Isaac Newton Telescope | 2003.613 | $36.478\pm 0.013$ | $42.012\pm 0.014$
Pan-STARRS1 | 2011.772 | $36.476\pm 0.015$ | $41.962\pm 0.018$
Gaia DR3 | 2016.000 | $36.459\pm 0.001$ | $41.937\pm 0.001$
UJA Telescope | 2019.433 | $36.443\pm 0.042$ | $41.884\pm 0.041$
Nordic Optical Telescope | 2019.441 | $36.422\pm 0.061$ | $41.936\pm 0.059$
UJA Telescope | 2021.763 | $36.431\pm 0.030$ | $41.915\pm 0.032$
111(a) In the ICRS. Errors given are one standard deviation.

##  2 Serendipitous discovery and initial observational follow-up

While inspecting the bow shock of BD+43∘3654, a Cygnus massive runaway star
(Benaglia et al., 2010), in the All-Sky Data Release of the Wide-field
Infrared Survey Explorer (WISE), we detected an unusual filamentary structure,
nearly one degree long and with a measured position angle along
$178^{\circ}\pm 2^{\circ}$ , with a conspicuous star-like source at a distance
of about $0.35$ arcmin from the frontal sharpest edge (Fig. 1). This anonymous
object was later realized to have been present, but gone unnoticed in
previously published images of the field (Toalá et al., 2016). It also
corresponds to entry 2070567522539111168 in the third Gaia data release (DR3),
where it is flagged as variable. No proper motion or parallax information is
available, however. Intrigued by these facts, we first conducted intensive
astrometric and photometric observations with the 0.4 m University of Jaén
Telescope (UJT; MPC code L83, Martí et al. (2017)), during which the nickname
UJT-1 was assigned. To overcome the shortage of Gaia parameters, astrometry
was performed on different archival images from nearly seven decades ago. This
was also complemented with astrometry based on modern CCD images including the
most recent image in 2021 with the UJT. The different image sources are
collected in the first column of Table 1.

Accurate astrometric solutions for all images, including the historical ones,
were established based on numerous Gaia reference stars in the field corrected
for proper motions. The IRAF package was used for bias, dark current, and
flat-field corrections when needed and to determine the plate solutions. In
particular, its tasks daofind, ccxymatch, ccmap, and cctran were key for this
last purpose. The resulting sky coordinates in the International Celestial
Reference System (ICRS) are presented in Table 1. They were computed using the
full astrometric solutions (usually up to third-order terms) to account for
plate distortions in the focal plane. This translates into astrometric
residuals of reference stars ranging from about 0.1 arcsecond for historical
photographic survey images to $\sim 0.01$ arcsecond for more recent electronic
frames selected among the best seeing conditions. The statistical errors given
in Table 1 typically correspond to about 1/10 to 1/20 pixel since the target
star was always detected with a good signal-to-noise ratio. The final outcome
are the proper motions ($\mu_{\alpha}\cos{\delta}=-1\pm 1$ mas yr-1 and
$\mu_{\delta}=-6\pm 1$ mas yr-1) that resulted from a weighted linear least-
squares fit to all Table 1 positions (see Fig. 2).

The UJT astrometric run in 2019 also included absolute photometry, yielding
$R=16.48\pm 0.06$ and $I=14.70\pm 0.03$ using Landolt standard stars (Landolt,
1992). In a similar way, imaging with the 1.23 m telescope at the Centro
Astronómico Hispano-Alemán in Calar Alto (CAHA, Spain) on 2019 December 10
provided $V=19.2\pm 0.1,$ but not a simultaneous detection in the blue ($B\geq
20.5$).

##  3 Spectral type and distance determination

Figure 3:  Optical spectrum of UJT-1. Data obtained with the NOT telescope.
The most prominent stellar spectral lines, H$\alpha$ in emission, and the CaII
triplet absorption are marked together with $H_{2}O$ and $O_{2}$ telluric
absorption features.

We also obtained optical spectroscopy of UJT-1 using the ALFOSC spectrograph
of the Nordic Optical Telescope (NOT) at the Observatorio del Roque de los
Muchachos (ORM) (Fig. 3). The spectra were taken on 2019 June 10 with a total
of 1800s exposure time, and using the ALFOSC grism 4 that covers the 3200-9600
Å region with medium resolution. Data reduction, also with IRAF, included bias
and flat-field correction followed by spectrum extraction with removal of sky
background. Wavelength calibration was achieved using Th-Ar lamps. Flux
calibration was tied to a nearly simultaneous observation of the standard star
BD+17 4708 (a subdwarf star of spectral type F8), at roughly the same air
mass. Unfortunately, the target lines in the blue region of the spectrum were
not accessible due to high interstellar extinction. From the continuum flux
level, we can still estimate the approximate but simultaneous magnitudes
$B=21.4\pm 0.2$ and $V=18.8\pm 0.1,$ indicating a color $B-V=2.6\pm 0.2$.
Beyond the highly reddened and absorbed continuum, the most distinctive target
feature was a noticeable, blueshifted H$\alpha$ emission component. The
H$\alpha$ energy flux and equivalent width amounted to $(4.1\pm 0.1)\times
10^{-15}$ erg s-1 cm-2 and $-10.7\pm 0.2$ Å. Other easily recognizable
spectral features were the near-infrared Ca II triplet in absorption. This
immediately indicates a late spectral-type classification for our target. A
second ALFOSC observation on 2019 October 18 using grism 19 confirmed the
H$\alpha$ emission and revealed no trace of lithium I 6707 Å.

Figure 4:  Comparison of the spectrum of UJT-1 with different spectral
templates. Templates are taken from a CaII triplet spectral library (Cenarro
et al., 2001). All spectra are rest-frame corrected.

In order to improve the spectral classification, the first NOT spectrum was
continuum rectified to be compared to a library of ionized calcium stellar
spectra (Cenarro et al., 2001) convolved to the same resolution. Under a
least-squares criterion, the closest match corresponded to a G5 III type (Fig.
4). We point out here that this refers only to a location in Hertzsprung-
Russell (HR) diagram and does not necessarily reflect the young evolutionary
status of the star. Based on the absolute magnitude and intrinsic $B-V$ color
at this HR position, the most plausible distance range to the source is
roughly estimated to be $d=4.5\pm 1.0$ kpc, with an interstellar absorption as
high as $A_{V}=4.8\pm 0.5$ magnitudes. This is equivalent to a color excess in
the range $E(B-V)=1.6\pm 0.2$ magnitudes. A G-dwarf location in the HR diagram
is ruled out because a very nearby distance (less than 1 kpc) would be
required, which is clearly inconsistent with the highly reddened spectrum in
Fig. 3. Based on the quality of the NOT medium-resolution spectrum, an F-dwarf
spectral classification could still be conceivable with similar $A_{V}$ and
$E(B-V)$ values. This would imply a closer location, in the range 1.5-2.0 kpc,
placing UJT-1 in the outskirts of the Cygnus X region. However, according to
the tridimensional reddening maps by Capitanio et al.
(2017)222https://stilism.obspm.fr, a distance of 2 kpc or higher is needed to
reach the estimated color excess for the UJT-1 line of sight. Therefore, the
farthest distance of 4.5 kpc appears to be slightly more favored and is
preferentially adopted in this work unless stated otherwise. Finally, a late-
type supergiant classification does not apply because it would imply a
distance far beyond the Milky Way limits.

A subsequent observation of UJT-1 on 2019 September 23 with the Isaac Newton
Telescope (INT), also at the ORM, and its intermediate dispersion
spectrograph, enabled a better measurement of the Ca triplet wavelengths. The
R1200R grism with 1200 s exposure time and similar data processing was used,
yielding a heliocentric radial velocity estimate of $-33\pm 4$ km s-1 under
the assumption that this is a single object.

##  4 Stellar run-away velocity

.

Figure 5:  Range of peculiar velocities and position angles of peculiar
motion. The gray and brown shaded areas cover the 90% and 50% confidence
regions that are consistent with the estimated proper motions. The vertical
lines illustrate the narrow interval of position angles that is consistent
with the sky orientation of the UJT-1 trail.

Based on the astrometric and distance parameters found above, it is possible
to estimate the peculiar velocity of UJT-1 with respect to its regional
standard of rest (RSR). For this purpose, the velocity components due to the
rotation of the Milky Way and the motion of the Sun in the local standard of
rest (LSR) have to be subtracted from the observed heliocentric velocity
resulting from the previous astrometric and spectroscopic measurements. The
equation formalism is well documented in the literature (Moffat et al., 1998)
and can be applied to both the tangential and radial velocity components. When
the peculiar velocities and associated proper motions are available, it is
also possible to compute the position angle of the true stellar motion
relative to its RSR, its line-of-sight inclination angle, and the peculiar
velocity modulus. In our case, the lack of accurate Gaia proper motions is a
handicap that renders our analysis more difficult than usual. To overcome this
limitation, we decided to explore the 90% confidence region around the
($\mu_{\alpha}\cos{\delta},\leavevmode\nobreak\ \mu_{\delta}$) values found
above, assuming that they follow uncorrelated Gaussian distributions. The low
($\sim 0.1$) absolute values of the corresponding covariance matrix elements
for Gaia stars in the field justifies this last assumption. The result of this
exploration is presented in Fig. 5, where the sampled parameter space overlaps
well with the observed position angle of the UJT-1 wake feature
($178^{\circ}\pm 2^{\circ}$). A 4.5 kpc distance was used here. The
corresponding modulus of the peculiar velocity relative to the RSR reaches
values in the range 15 to 80 km s-1 at the 90% confidence level and agrees
with the position angle of the wake. About 75% of this interval is well above
the usual threshold ($\sim 30$ km s-1) for the definition of runaway stars
(Eldridge et al., 2011). Therefore, UJT-1 is a strong candidate of this
category. An intermediate plausible value of 45 km s-1 is adopted for
discussion purposes, and the corresponding line-of-sight angle of $82^{\circ}$
is practically in the plane of the sky. Conversion from angular distance to
travel time is achieved using the factor $\mathscr{v}_{*}\sin{i}/d=2.3$ mas
yr-1.

When a shorter distance (2 kpc) is assumed, a plot similar to Fig. 5 is
obtained, but with a lower intermediate plausible velocity of 30 km s-1 with a
line-of-sight angle of $45^{\circ}$. This is still compatible with the runaway
classification. Therefore, UJT-1 remains consistent with being a high-velocity
object despite the distance uncertainties.

##  5 Spectral energy distribution of the star

Figure 6:  SED of the star UJT-1. It is based on CAHA + UJT photometry, NOT
spectroscopy, and different multiwavelengths surveys listed in the figure
inset. The thick line represents a tentative SED fit based on libraries of
theoretical spectra for YSO (dashed line) and dust (dotted line) components
(Draine & Li, 2007; Robitaille et al., 2006). Error bars correspond to one
standard deviation.

The SED of UJT-1 can be assembled from observations of our own work (UJT, NOT,
and CAHA), together with cataloged data obtained from different infrared
telescopes such as the Two Micro All Sky Survey (2MASS), WISE, AKARI, Spitzer,
Midcourse Space Experiment (MSX), and Planck. The resulting SED (Fig. 6) is
dominated by two clear maxima that are usually interpreted in terms of a
central protostar and disk or ambient dust emission, respectively. This figure
also includes a tentative attempt to model fit the data points using libraries
of young stellar objects (YSOs) and thermal dust spectra (Draine & Li, 2007;
Robitaille et al., 2006). Although this fitting exercise was hampered by high
optical extinction, strong variability (Heinze et al., 2018), and poor angular
resolution at long wavelengths, a plausible agreement with a YSO
interpretation is supported. In this context, it was difficult to distinguish
among spectral templates, but acceptable fits with the central protostar
temperature in agreement with a late spectral type (below 7500 K) were
obtained.

At X-ray energies, a marginal detection of UJT-1 was obtained from standard
processing of archival XMM-Newton data (Obs. Id. 0653690101) using the science
analysis system of this observatory (SAS). Assuming a simple power-law
spectrum, with the hydrogen column density set according to $A_{V}=4.5$ mag,
we estimated the 0.5-8 keV luminosity to be $3\times 10^{32}$ erg s-1 $(d/{\rm
4.5\leavevmode\nobreak\ kpc})^{2}$. Because our knowledge of the distance is
limited, this places UJT-1 at the high end of or well within average typical
TTS X-ray luminosities. Moreover, it indicates an object with a few solar
masses based on the known X-ray luminosity versus mass correlation for TTSs
(Preibisch et al., 2005).

##  6 Spectral energy distribution of the bow-shock tail and dust properties

Figure 7:  Averaged SED for the different regions along the limb-brightened
bow-shock shell. The proposed tentative fit is plotted as the continuous red
line from a dust grain model. The dotted line shows the fit of a single
radiation field to the largest wavelengths. Error bars correspond to one
standard deviation. See text for details.

The infrared emission from the whole structure is mainly due to heated dust in
the ISM that is swept up to the shocked layer. WISE and Spitzer observations
are taken at wavelengths corresponding to small particles (0.001 - 0.01
$\mu$m) that are stochastically heated, while AKARI and Planck observations
detect emission from colder, larger particles ($\geq 0.1$ $\mu$m) in thermal
equilibrium with the starlight. From these satellite images, we built the SED
of the bow shock from 3.35 to 550 $\mu$m (see Fig. 7) for ten different
regions along the limb-brightened bow shock, each of them with $\sim 2$% of
the whole solid angle with emission. We found a clear constancy in the shape
of the spectra, so that mean values are shown. To obtain an order of magnitude
of the dust properties, we first tried to fit the $\lambda\geq 100$ $\mu$m
part of the spectrum, which is attributed to the larger particles, because
they cause most of the infrared emission and dust mass. The simplified model
we used here (Dwek & Werner, 1981) is based on a modified blackbody with
emission efficiency $\sim\lambda^{-1}$ and a density of dust grains of 3 g
cm-3. For each of our ten regions, the fit gives a dust mass $M_{d}\sim 0.03$
$M_{\odot}$ and a radiative equilibrium temperature of $T_{d}\sim 35$ K. When
we take into account that the analyzed regions contain about 20 % of the total
shocked volume, the overall dust mass in the tail is estimated to be about
$M_{d}\sim 1.7$ $M_{\odot}$. When the distance to the star is reduced to 2
kpc, the total amount of dust mass becomes about $M_{d}\sim 0.34M_{\odot}$.

Using the same model (Dwek & Werner, 1981), we can estimate the infrared
luminosity for each of the ten regions we analyzed. As a result, a total
amount of $L_{IR}\sim 5.4\times 10^{37}$ erg s-1 $\left[d/{\rm
4.5\leavevmode\nobreak\ kpc}\right]^{2}$ is estimated for the whole tail
emission. We also tried to fit the complete SED to a more sophisticated dust
model (Draine & Li, 2007) that considers the contribution from a size
distribution of grains of different composition (including emission lines from
polycyclic aromatic hydrocarbons, or PAHs) exposed to a range of radiation
intensities. We find the best agreement for a dust distribution with 1.49% in
mass of PAH particles containing fewer than 1000 carbon atoms, and a range of
intensities from 1 to $10^{4}$ times the starlight in the solar neighborhood
ISM. However, Fig. 7 displays an excess in the $\sim 5$-$20$ $\mu$m range that
cannot be attributed to the dust model. This appears at all epochs of the bow-
shock shell and might be tentatively attributed to a new light dust component
that in this wavelength might crudely be approximated interval by a modified
blackbody with a $\lambda^{-2}$ emissivity law at $T\sim 300$ K (not shown in
the figure).

##  7 Wake growth scaling determination and period search

Figure 8:  Detailed view of the bow shock and model fit of its geometry. The
UJT-1 star is the brightest object on the left side of this IRAC Spitzer image
at the 3.6 $\mu$m wavelength. The fitted Wilkin profile is overplotted in
blue. The one-arcminute scale bar is equivalent to about 1.3 pc at the
proposed distance of 4.5 kpc.  Figure 9:  Scaling law for the growth of the
width of the bow-shock wake. The width is defined as the distance between the
two border points of the wake at the same distance from the star. The power-
law fit with exponent $0.33\pm 0.02$0 is shown as the solid orange line. Error
bars represent an estimated standard deviation of about 6 arcsecond.

The infrared wake behind the protostar grows with the distance from UJT-1, as
seen in the detailed view of Fig. 8. To determine the scaling law of this
behavior, we used the WISE band 3 image as the best image. Angular widths in
this section are expressed in arcseconds. The wake axis was aligned with the
observed position angle. Its border pixels were selected by visual inspection,
which gives similar values values as automated extraction, but they are more
accurate. The estimated uncertainty is about 6 arcseconds (a few pixels). We
determined the width $2R$ of the wake at a certain distance $z$ in from the
star by measuring the span between two pixels. Because we searched for a self-
similar scaling, only points far from the head of the bow shock ($z>300$
arcseconds) were considered. The determination coefficient of the resulting
power-law scaling $2R=a(z/z_{0})^{b}$ is $r^{2}=0.991$, with $a=28\pm 3$
arcseconds and $b=0.33\pm 0.02$ (see Fig. 9). Here the reference value $z_{0}$
is one arcsecond. These data were also used for a space periodicity analysis
with the phase dispersion minimization method (Stellingwerf, 1978). As a
result, we obtained a possible 4-5 arcminute wavelength periodicity, which for
the projected $\mathscr{v}_{*}$ and the proposed distance gives a shedding
frequency $f\sim 10^{-12}$ Hz.

##  8 Discussion

Based on the H$\alpha$ emission line with an equivalent width above 10 Å (in
absolute value) and the SED of the source (Fig. 6) with a spectral index
$-1.33\pm 0.03$ (in the 4-22 $\mu$m range), UJT-1 fulfills all the common
taxonomy criteria for a class II YSO, or classical TTS (Appenzeller & Mundt,
1989; Lada, 1987). Additional confirmation of the TTS nature comes from the
location of UJT-1 in the color-color diagram of Fig. 10 based on the work by
Meyer et al. (1997), which was created using the 2MASS photometry as starting
point. In this plot, the red cross displays the observed UJT-1 infrared
colors, and the blue cross corresponds to the estimated intrinsic colors after
correcting for the estimated interstellar absorption ($A_{V}=4.8$ mag). It is
reassuring that UJT-1 lies almost perfectly near the so-called TTS locus. For
further discussion, we can then assume typical TTS values for the mass and
radius of the star of $M_{*}\sim 2$ $M_{\odot}$ and $R_{*}\sim 2$ $R_{\odot}$.

Figure 10: Infrared color-color diagram showing the location of UJT-1 before
(red cross) and after (blue cross) correction for interstellar absorption. The
dashed line corresponds to the TTS locus as defined by Meyer et al. (1997).
Black squares and triangles indicate the location of ordinary dwarf and giant
stars according to Bessell & Brett (1988). The three parallel lines indicate
the direction of the reddening vectors in steps of $\Delta A_{V}=1$ mag
following the Rieke & Lebofsky (1985) extinction law.

###  8.1 Bona fide RATTS

At this point, from the obtained evidence, we propose that UJT-1 is an
excellent prototype of RATTS systems. This also supports the hypothesis that
they are a class. This is further reinforced by the remarkably long path
traced by the star along its trajectory, which is shown at infrared
wavelengths in Fig. 1. The path consists of a detached shell surrounding the
star together with a wavy, asymmetric, limb-brightened tail that extends
downward by almost one degree. An abrupt width change is observed in the
middle of the structure, and the tail is finally distorted and diffuses into
the interstellar medium (ISM). The closer view in Fig. 8 also highlights some
wiggles and disturbances. If it is located at 4.5 kpc and when we consider a
line-of-sight angle almost in the plane of sky, the observed wake would extend
over the considerable distance of$\leavevmode\nobreak\ 50$ pc. This exceeds
the size of the remarkable turbulent wake behind the famous star Mira by more
than a factor of 10 (Martin et al., 2007). The length of the UJT-1 bow shock
then translates into a crossing lifetime of $\sim 1$ Myr, which is compatible
with the age of a typical TTS (see also Appendix A about its original natal
site). The absence of lithium in NOT spectroscopy consistently suggests that
UJT-1 is not extremely young. It is highly remarkable that the effect of a
modest, individual star in the ISM appears to be traceable so many parsecs
away, and this is rivalled in the Milky Way only by the more energetic
relativistic jets of microquasar sources (Fabrika, 2004).

The strength of these long-range effects needs to be revised if UJT-1 is
closer to us. Even at 2 kpc, the observed wake would remain more than four
times longer than its Mira counterpart, however.

###  8.2 Mechanical scenario

This notorious large wake feature resembles the familiar turbulent wake
observed behind common bodies placed in a free stream, characterized by a
self-similar behavior beyond a certain downstream distance $z$. In particular,
this turbulent wake width scales according to Landau & M. (1987) as $z^{1/3}$.
Remarkably, this is the observed growing rate in UJT-1 tail (with a power-law
index $0.33\pm 0.02$; see Section 7). Thus, we might be tempted to assume that
it is a clear example of a classical turbulent wake. We might even think of a
possible TTS in a bright rimmed cloud (Hosoya et al., 2019). However, the
scenario is much more complex. The forming of a detached shell in the front of
the star indicates that it moves at supersonic velocities, but the star itself
also emits strong, supersonic winds. Therefore, from the point of view of the
star, the incoming ISM and the wind collide in a contact discontinuity in
which the gas and dust accumulate (van Marle et al., 2011). The heated ISM and
wind material expand outward from either side of the discontinuity, the former
up to a front bow shock ahead, the latter into the stellar wind to form a
reverse shock. This is a common scenario in the Universe and has been
reported, among others, in runaway massive, early-type stars (van Buren &
McCray, 1988), evolved supergiants (Jorissen et al., 2011), asymptotic giant
branch (AGB) stars (Martin et al., 2007; Noriega-Crespo et al., 1997), or even
pulsars (Kargaltsev et al., 2017). It is unprecedented in pre-main-sequence
stars like ours to this extent, however.

The mechanical scenario of a runaway star moving supersonically in a uniform
ISM depends on its velocity in the ISM rest frame, $\mathscr{v}_{*}$, the
stellar wind terminal velocity, $\mathscr{v}_{w}$, the stellar wind mass-loss
rate, $\dot{M}_{w}$, and the density of the ambient ISM, $\rho_{a}$. These
variables define a characteristic length, the stand-off distance from the star
to the bow-shock apex,
$R_{0}\sim(\mathscr{v}_{*}^{2}\dot{M}_{w}\rho_{a}^{-1}\mathscr{v}_{w}^{-2})^{1/2}$,
and the governing parameter $\eta=\mathscr{v}_{w}/\mathscr{v}_{*}$. It is also
possible to determine a characteristic dynamical timescale
$t_{dyn}=R_{0}/\mathscr{v}_{*}$. In addition, to determine the whole physical
scenario, we must consider the cooling properties of the ISM and wind shocked
layers (Comerón & Kaper, 1998), so that a new characteristic velocity arises,
the sound speed $c_{s}$, and its corresponding governing parameter, the Mach
number $\mathcal{M}=\mathscr{v}_{*}/c_{s}$. In the case of shocked gas, we
must use the post-shocked velocity. The characteristic timescale is therefore
better defined as $t_{dyn,s}=4R_{0}/\mathscr{v}_{*}$. Cooling timescales may
also be defined for both the shocked ISM and the wind (Appendix B). When any
of these cooling timescales is shorter than the shocked dynamical timescale,
we can assume that the corresponding layer is radiative, and it is adiabatic
in the opposite case.

In order to test these hypotheses, we need some numerical measurements and
estimates. As stated in Section 4, a distance of 4.5 kpc and a peculiar
velocity of $\mathscr{v}_{*}=45$ km s-1 is preferentially adopted. Then, the
stand-off distance from infrared images of UJT-1 (e.g., Fig. 1), subtending
$\sim 0.35$ arcminute, is equivalent to $R_{0}\sim 0.46$ pc given the line-of-
sight angle as well. . On the other hand, we consider a neutral Galactic ISM
with number density $n_{H}=2.5$ ${\rm cm}^{-3}$, which is slightly higher than
the typical 1 ${\rm cm}^{-3}$ because of the cloudy environment that is shown
in infrared maps. We then estimate $\rho_{a}=\mu n_{H}=5.0\times 10^{-24}$ g
cm-3, with $\mu=2.3\times 10^{-24}$ g being the adopted mean gas mass per
hydrogen atom. Finally, as measuring $\mathscr{v}_{w}$ for a TTS is
challenging, here we estimate $\mathscr{v}_{w}\sim 0.5\mathscr{v}_{esc}\sim
450$ km s-1, where $\mathscr{v}_{esc}$ is the escape velocity, while the
prefactor that accounts for the uncertainties in the acceleration processes is
taken for a main-sequence star of the same spectral type as UJT-1 (Eldridge et
al., 2006). This $\mathscr{v}_{w}$ value is about the same magnitude as the
winds of several hundredths km s-1 observed in TTSs (Kwan et al., 2007), and
it gives $\eta\sim 10$.

We thus obtain a cooling time for the shocked ISM
$\tau_{c,\leavevmode\nobreak\ sISM}\simeq 2.9\times 10^{9}$ s, and
$\tau_{c,\leavevmode\nobreak\ sw}\simeq 1.7\times 10^{14}$ s for the shocked
wind, while $\tau_{dyn,\leavevmode\nobreak\ s}\simeq 1.3\times 10^{12}$ s. As
a result, the bow-shock shell probably has an adiabatically shocked wind
layer, but a radiatively shocked ISM layer. This is also true when the star is
located at $d\simeq 2$ kpc, where the stand-off distance would become smaller
($R_{0}\sim 0.29$ pc) and $\eta=15$, but the cooling and dynamical times are
only slightly different. Although the real shape of the bow shock in general
has to be obtained numerically (Comerón & Kaper, 1998), in the particular case
of fully radiatively shocked layers, both are very thin, and it is possible to
use the Wilkin analytical solution for the shape of the bow shock (Wilkin,
1996). Under the same hypotheses, Wilkin also obtained a scaling law for the
radius of the shell $\tilde{R}=(3\pi\tilde{z})^{1/3}$ for a high enough
$\tilde{z}$. The tilde denotes dimensionless lengths relative to $R_{0}$. This
is remarkable because it is exactly the same self-similar solution for the
turbulent wake, and it explains why the tail in UJT-1 mimics this structure.
Because the actual shape in our bow shock is so similar to this $z^{1/3}$
scaling, we suspect that the Wilkin profile might be an acceptable solution in
our case in spite of the adiabatic wind layer. The fit result is shown in Fig.
8, with a 3D revolution form of the Wilkin profile superposed on a Spitzer
image, considering the line-of-sight angle stated previously. The bow shock
borders seem to fit the Wilkin prediction well, except for the very final part
of the tail, which may be due to buoyancy, wind drag, or the interaction with
the limits of a molecular cloud that might harbor the UJT-1 birthplace
(Appendix A). In any case, the fit in Fig. 8 enables us to further refine the
stand-off distance to $R_{0}=(0.46\pm 0.12)$ pc at the preferred distance,
accurately accounting for projection effects. Some fit residuals left still to
be accounted for, however.

###  8.3 Instabilities in the flow

For instance, instabilities may affect the shocked ISM layer because it is
radiative and therefore thin. The perturbations would originate at the apex of
the bow shock and would later be driven downward by the flow. Among the
instabilities that may be at work, the more typical are the Kelvin-Helmholtz
(KH) and Rayleigh-Taylor (RT) instabilities. Other possible instabilities in
bow shocks are the so-called transverse acceleration (TA) (Dgani et al., 1996)
and the nonlinear thin-shell (NTS) accelerations (Vishniac, 1994). We can use
timescale arguments to estimate the possibility that the instabilities to grow
quickly enough to affect the layer, as detailed in Appendix C. As a result, in
our case, KH and NTS instabilities may cause small disturbances in the apex of
the bow shock that are advected downward and amplified by TA instabilities,
but it seems difficult for RT instabilities to contribute. Notwithstanding, we
cannot strictly rule out the possibility that some of these wavy patterns may
be originated from the interaction of the shocked layers with small clumps of
different densities in a non-uniform ISM, as has been numerically tested in
similar scenarios (Toropina et al., 2019).

Figure 11:  Bow-shock limits according to the power law $z^{1/3}$. The dashed
cyan line corresponds to the current parameters of UJT-1, and the yellow line
shows a bow shock assuming an ISM density that an order of magnitude lower.
The background grayscale image is again the Spitzer view at 3.6 $\mu$m.

However, these instabilities are on the order of the width of the shocked
layer, so that larger departures from the theoretical shape observed in the
UJT-1 tail might be attributed to other causes. One possible mechanism is the
vortex shedding (see Appendix D), which has been proposed in some simulations
of AGB star wakes (Wareing et al., 2007). Although this cannot be strictly
ruled out from the calculations, a detailed view of the images does not
support the typical regular, periodic shedding of vortices, as was already
pointed out (Wareing et al., 2007). This irregular shedding might be
attributed to a fluctuating bow-shock shape, but the whole picture resembles a
turbulent behavior more closely, especially considering the uncertainties in
the viscosity estimate. Therefore, another mechanism for the large wavelength
distortions in the bow-shock shape is probably at work, especially at the
strong, abrupt change in width at approximately the middle of the tail (Fig.
11). Here, the bow-shock radius so suddenly expands by a factor of
$R_{1}/R_{2}\sim 2$ that the more natural explanation seems to be the result
of UJT-1 passing through the frontier between two regions of different
density. This rapid contraction is very similar to the contraction obtained in
simulations of AGB star bow shocks passing through different ISM densities
(Toropina et al., 2019; Yoon & Heinz, 2017). Assuming the bow shocks are in
pressure equilibrium, and that the ISM density is higher at the current
position of UJT-1, then $c_{s}$ there is lower, $\mathcal{M}$ is higher, and
the corresponding Mach cone is narrower, as it appears. In Fig. 11, two
different power-law $z^{1/3}$ Wilkin bow-shock profiles are fit to the
infrared Spitzer image. As the radius jump occurs so abruptly, it must evolve
as $\sim(n_{1}/n_{2})^{1/3}$, so that the required ISM density increase is
about one order of magnitude. Moreover, these anisotropic properties of the
ISM around UJT-1 might explain (Yoon & Heinz, 2017) the observed east-west
asymmetry in infrared emission. The stand-off distance obtained from the
Wilkin relation $\tilde{R}=(3\pi\tilde{z})^{1/3}$ above, together with the
fitting of the wake to $\tilde{R}\sim\tilde{z}^{1/3}$, results in $R_{0}$ of
about the values reported in the previous section.

###  8.4 Role of mass-loss variability

Figure 12:  Reconstruction of the past wind mass-loss in the TTS UJT-1 as
inferred from its trailing wake. The plot is scaled according to the currently
estimated value of $\dot{M}_{w0}=1.8\times 10^{-7}$ $M_{\odot}$ yr-1. Error
bars indicate one standard deviation. Mass-loss changes by a factor of $\sim
5$ over $\sim 0.5$ Myr are recovered.

On the other hand, intermediate-scale disturbances in the bow-shock profile
might also be attributed to actual changes in the stellar wind mass-loss rate.
Its current value can be estimated from the measured $R_{0}$ considering the
balance of ISM ram pressure and stellar wind flux momentum as

$\dot{M}_{w}=\dot{M}_{w0}\left[\frac{n_{H}}{2.5\leavevmode\nobreak\ {\rm
cm}^{-3}}\right]\left[\frac{R_{0}}{0.46\leavevmode\nobreak\ {\rm
pc}}\right]^{2}\left[\frac{\mathscr{v}_{*}}{45\leavevmode\nobreak\ {\rm
km\leavevmode\nobreak\ s}^{-1}}\right]\left[\frac{\eta}{10}\right]^{-1},$ (1)

where $\dot{M}_{w0}$ is the current value of the mass loss ($1.1\times
10^{-6}$ $M_{\odot}$ yr-1 for d = 4.5 kpc and $1.8\times 10^{-7}$ $M_{\odot}$
yr-1 for d = 2 kpc). These values are consistent with the upper limit of
observed outflows and neutral winds in TTSs (Hartigan et al., 1995; Ruiz et
al., 1992). Assuming that the ISM density remains constant up to the sudden
radius jump above and that the star and wind velocities remain constant as
well, we can estimate the mass-loss rate at any distance $z$ far enough from
the star for the Wilkin power-law equation (Wilkin, 1996) to apply. If
$\dot{M}_{w}$ is the current mass-loss rate and $R^{\prime}$, $R$ are the
measured and theoretical bow-shock radii at a certain distance, the
corresponding mass-loss is
$\dot{M}^{\prime}_{w}=\dot{M}_{w}(R^{\prime}/R)^{3}$. Therefore, in a
variable-wind scenario, we can reconstruct the stellar wind mass-loss history
over a huge time lapse encoded in the remarkably long wake of UJT-1. A plot of
the result is shown in Fig. 12, where we go back $\sim 5\times 10^{5}$ yr in
the past, witnessing the history of the mass-loss evolution of a pre-main-
sequence star with unprecedented coverage and detail.

###  8.5 Toward a consistent large-scale view

The most intriguing characteristic of the wake structure of UJT-1 is its very
extended length. This suggests that the star is passing through ISM regions
that are devoid of strong winds that could destroy its long-lasting
morphology. In addition, the observed infrared emission must remain
persistent. In bow shocks, this radiation is dominated by thermal dust that
was previously heated by absorption of starlight photons and collisions with
ions and electrons (Draine, 2003). In Section 6, we obtained SEDs from 3.35 to
550 $\mu$m for different regions along the bow-shock tail using available
data. The shape of the spectrum is highly consisten, revealing that the
smaller-grain dust-gas collisions are maintained in time. Large-grain ($\geq
0.1$ $\mu$m) emission with $T_{d}\sim 35$ K causes the far-infrared part of
the SEDs. It is noteworthy that the fitted temperature is slightly higher than
that expected at a distance $\sim R_{0}$ from our TTS ($R_{*}\sim 2$
$R_{\odot}$, and effective temperature $T_{\rm eff}\sim 5400$ K for the G5
spectral type) given by Dwek & Werner (1981)
$T_{d}=T_{eff}(R_{*}/R_{0})^{2/(4+s)}=8$-$25$ K, with $s$ being the slope of
the opacity in the infrared regime, which varies from 1 (silicates) to 2
(graphite). The dust mass in the wake is $M_{d}\sim 1.7$ $M_{\odot}$, similar
to the snowplow model estimate of $M_{d}\sim 1.9$ $M_{\odot}$ obtained
assuming that the wind sweeps up all the dust inside the bow-shock volume, as
calculated from the Wilkin analytical model for a total length of the shell of
$\tilde{L}\sim 100$, and a constant dust-to-gas ratio of 0.01, for which the
gas density and stand-off distance are also constant. If the star were an F0
dwarf with $T_{\rm eff}\sim 7400$ K and located at 2 kpc, the dust temperature
estimated with the Dwek & Werner (1981) relation above would be between 14 and
40 K, in agreement with the temperature obtained from fitting the SED. The
total dust mass in the wake would then be $M_{d}\sim 0.34$ $M_{\odot}$, again
in accordance with the mass obtained with the snowplow approximation,
$M_{d}\sim 0.26$ $M_{\odot}$. However, the obscured appearance of the UJT-1
wake might render the higher dust-mass values preferable.

Figure 13:  Neutral hydrogen cavity west of the UJT-1 stellar wake. Map of the
brightness temperature as detected in the EBHIS (Winkel et al., 2016) in the
channel corresponding to an LSR velocity of $-14$ km s-1. A deep minimum lies
west near the UJT-1 stellar wake. The blue-green map corresponds to the radio
emission at 20 cm from CGPS, which might be attributed to an unpublished old
supernova remnant that is likely to interact with dust particles in the wake.

The infrared luminosity of the tail we derived in Section 6 is not consistent
with the upper limit emission assuming a constant mass-loss rate and the local
ISM thermalizing the whole kinetic energy when it collides with the stellar
wind in the bow shock,
$L_{IR}<\frac{1}{2}\dot{M}_{w}(\mathscr{v}_{*}^{2}+_{w}^{2})\sim 6.9\times
10^{34}$ erg s-1, which changes to $L_{IR}<1.2\times 10^{34}$ erg s-1 when the
distance to the star is only 2 kpc. In spite of the crudeness of the
calculations, the UJT-1 tail appears to be heated from external sources. As no
clear O/B stars are identified in the near vicinity of the region, we
speculate that a close encounter with the shell of an unnoticed supernova
remnant causes the maintained infrared emission of the dust in the wake and
its estimated temperature, which is slightly above the 8-25 K that our TTS
could maintain at a distance $R_{0}$. The supernova remnant may be the
asymmetric shell that was detected in radio maps from the Canadian Galactic
Plane Survey (CGPS), which remarkably matches a deep minimum in HI emission
from the Effelsberg-Bonn HI Survey (EBHIS) at a similar kinematic distance as
the proposed UJT-1 natal cloud (Fig. 13). Therefore, the interaction might be
plausible, although it has to be very recent because the Wilkin $R\sim
z^{1/3}$ structure of the wake is not lost. On the other hand, the proposed
new small, light dust component at 5-20 $\mu$m found by fitting the Draine &
Li (2007) model should be consistent with metallic or carbon grain
composition, rather than dielectric or silicate composition (Bohren & R.,
1983), and should remain in the shocked layers of the bow shock for a long
time at $\sim 300$ K, which is the temperature expected in a dust tha tis
stochastically heated by hot gas ($\sim 10^{7}$ K as in the shocked wind
layer) (Dwek & Werner, 1981) .

##  9 Conclusions

We summarize our conclusions grouped into two classes: solid findings, and
tentative interpretations. As our first solid finding, we reported UJT-1 as an
outstanding example of a RATTS that moves supersonically and leaves a
conspicuous bow-shock shell behind with dimensions that are rarely found in
stellar contexts. This remains true despite the uncertainty in distance, which
ranges from 2 to 4.5 kpc. Second, this fast-moving TTS is seen through an
optical extinction of about five magnitudes. After we corrected for this
effect, its location in the infrared color-color diagram is shifted toward the
well-known TTS locus, thus confirming its YSO nature. Third, the discovered
system provides an excellent laboratory for studying large-scale interactions
of turbulent wind and the ISM. In particular, the bow-shock profile grows
remarkably close to the predictions of the Wilkin model, thus justifying a
connection between the star and the wake.

As more speculative distance- and model-dependent conclusions, we find that an
extrapolation of the UJT-1 motion backward apparently leads to a nearby
molecular cloud at roughly the high end of the distance range. In addition,
the possibility that intermediate-scale disturbances are due to different
types of instabilities, vortex shedding, or ISM clump interaction appears to
be conceivable, but is difficult to distinguish. Alternatively, a simpler
scenario is to invoke a TTS with a time-variable stellar wind. If this is the
case, the reported system has the potential to open a window of hundreds of
thousands of years over the mass-loss history in YSOs that otherwise would
remain inaccessible. Finally, we also propose that the large-scale conditions
of UJT-1 involve an external heating source that accounts for the hot-dust
component required to better fit the SED of the bow-shock tail. This is
tentatively attributed to a possible supernova remnant, suggested by a neutral
hydrogen cavity and extended radio features in the field, but it remains to be
confirmed.

###### Acknowledgements.

We acknowledge support by Agencia Estatal de Investigación of the Spanish
Ministerio de Ciencia, Innovación y Universidades grant
PID2019-105510GB-C32/AEI/ 10.13039/501100011033 and FEDER “Una manera de hacer
Europa”, as well as Programa Operativo FEDER 2014-2020 Consejería de Economía
y Conocimiento de la Junta de Andalucía (Refs. A1123060E00010 and FQM-322).
This publication makes use of data products from the Wide-field Infrared
Survey Explorer, which is a joint project of the University of California, Los
Ángeles, and the Jet Propulsion Laboratory/California Institute of Technology,
funded by the National Aeronautics and Space Administration. This work is
partly based on observations made in the Observatorios de Canarias del IAC
with the Nordic Optical Telescope and the Isaac Newton Telescope operated on
the island of La Palma by NOTSA and the Isaac Newton Group of Telescopes, and
on observations collected at the Centro Astronómico Hispano-Alemán (CAHA) at
Calar Alto, operated jointly by Junta de Andalucía and Consejo Superior de
Investigaciones Científicas (IAA-CSIC). Authors are also grateful to all other
astronomical data repositories used here that space limits do not allow to
enumerate. We additionally express our gratitude to our colleagues José M.
Torrelles (ICE-CSIC) and Luis F. Rodríguez (IRyA-UNAM) for valuable
discussions during this work.

## References

  * Appenzeller & Mundt (1989) Appenzeller, I. & Mundt, R. 1989, A&A Rev., 1, 291
  * Benaglia et al. (2010) Benaglia, P., Romero, G. E., Martí, J., Peri, C. S., & Araudo, A. T. 2010, A&A, 517, L10
  * Bessell & Brett (1988) Bessell, M. S. & Brett, J. M. 1988, PASP, 100, 1134
  * Bird (1994) Bird, G. A. 1994, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Oxford: Calendon Press)
  * Blondin & Koerwer (1998) Blondin, J. M. & Koerwer, J. F. 1998, New A, 3, 571
  * Bohren & R. (1983) Bohren, C. F. & R., H. C. 1983, Absorption and Scattering of Light by Small Particles (New York: Wiley)
  * Capitanio et al. (2017) Capitanio, L., Lallement, R., Vergely, J. L., Elyajouri, M., & Monreal-Ibero, A. 2017, A&A, 606, A65
  * Cenarro et al. (2001) Cenarro, A. J., Cardiel, N., Gorgas, J., et al. 2001, MNRAS, 326, 959
  * Comerón & Kaper (1998) Comerón, F. & Kaper, L. 1998, A&A, 338, 273
  * Comerón et al. (1993) Comerón, F., Torra, J., Jordi, C., & Gómez, A. E. 1993, A&AS, 101, 37
  * Dame et al. (2001) Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792
  * Dgani et al. (1996) Dgani, R., van Buren, D., & Noriega-Crespo, A. 1996, ApJ, 461, 927
  * Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241
  * Draine & Li (2007) Draine, B. T. & Li, A. 2007, ApJ, 657, 810
  * Dwek & Werner (1981) Dwek, E. & Werner, M. W. 1981, ApJ, 248, 138
  * Eldridge et al. (2006) Eldridge, J. J., Genet, F., Daigne, F., & Mochkovitch, R. 2006, MNRAS, 367, 186
  * Eldridge et al. (2011) Eldridge, J. J., Langer, N., & Tout, C. A. 2011, MNRAS, 414, 3501
  * Fabrika (2004) Fabrika, S.and Liftshitz, E. M. 2004, Astrophys. Space Phys. Res., 12, 1
  * Florsheim et al. (1965) Florsheim, B. H., Goldburg, A., & Washburn, W. K. 1965, AIAA Journal, 3, 1332
  * Hartigan et al. (1995) Hartigan, P., Edwards, S., & Ghandour, L. 1995, ApJ, 452, 736
  * Hattori et al. (2019) Hattori, K., Valluri, M., Castro, N., et al. 2019, ApJ, 873, 116
  * Heinze et al. (2018) Heinze, A. N., Tonry, J. L., Denneau, L., et al. 2018, AJ, 156, 241
  * Hoff et al. (1998) Hoff, W., Henning, T., & Pfau, W. 1998, A&A, 336, 242
  * Hosoya et al. (2019) Hosoya, K., Itoh, Y., Oasa, Y., Gupta, R., & Sen, A. K. 2019, International Journal of Astronomy and Astrophysics, 9, 154
  * Jorissen et al. (2011) Jorissen, A., Mayer, A., van Eck, S., et al. 2011, A&A, 532, A135
  * Kargaltsev et al. (2017) Kargaltsev, O., Pavlov, G. G., Klingler, N., & Rangelov, B. 2017, Journal of Plasma Physics, 83, 635830501
  * Kwan et al. (2007) Kwan, J., Edwards, S., & Fischer, W. 2007, ApJ, 657, 897
  * Lada (1987) Lada, C. J. 1987, in Star Forming Regions, ed. M. Peimbert & J. Jugaku, Vol. 115, 1
  * Landau & M. (1987) Landau, L. V. & M., L. E. 1987, Fluid Mechanics (Oxford: Pergamon Press)
  * Landolt (1992) Landolt, A. U. 1992, AJ, 104, 340
  * Li et al. (2021) Li, Y.-B., Luo, A. L., Lu, Y.-J., et al. 2021, ApJS, 252, 3
  * Martí et al. (2017) Martí, J., Luque-Escamilla, P. L., & García-Hernández, M. T. 2017, Bulgarian Astronomical Journal, 26, 91
  * Martin et al. (2007) Martin, D. C., Seibert, M., Neill, J. D., et al. 2007, Nature, 448, 780
  * McKee & Ostriker (2007) McKee, C. F. & Ostriker, E. C. 2007, ARA&A, 45, 565
  * Meyer et al. (2014) Meyer, D. M. A., Mackey, J., Langer, N., et al. 2014, MNRAS, 444, 2754
  * Meyer et al. (1997) Meyer, M. R., Calvet, N., & Hillenbrand, L. A. 1997, AJ, 114, 288
  * Moffat et al. (1998) Moffat, A. F. J., Marchenko, S. V., Seggewiss, W., et al. 1998, A&A, 331, 949
  * Neuhäuser (1997) Neuhäuser, R. 1997, Science, 276, 1363
  * Neuhäuser et al. (2020) Neuhäuser, R., Gießler, F., & Hambaryan, V. V. 2020, MNRAS, 498, 899
  * Neuhäuser et al. (1996) Neuhäuser, R., Magazzu, A., Sterzik, M. F., et al. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 109, Cool Stars, Stellar Systems, and the Sun, ed. R. Pallavicini & A. K. Dupree, 433
  * Neuhauser et al. (1997) Neuhauser, R., Torres, G., Frink, S., Covino, E., & Alcala, J. M. 1997, Mem. Soc. Astron. Italiana, 68, 1061
  * Noriega-Crespo et al. (1997) Noriega-Crespo, A., van Buren, D., Cao, Y., & Dgani, R. 1997, AJ, 114, 837
  * Preibisch et al. (2005) Preibisch, T., Kim, Y.-C., Favata, F., et al. 2005, ApJS, 160, 401
  * Rieke & Lebofsky (1985) Rieke, G. H. & Lebofsky, M. J. 1985, ApJ, 288, 618
  * Robitaille et al. (2006) Robitaille, T. P., Whitney, B. A., Indebetouw, R., Wood, K., & Denzmore, P. 2006, ApJS, 167, 256
  * Ruiz et al. (1992) Ruiz, A., Alonso, J. L., & Mirabel, I. F. 1992, ApJ, 394, L57
  * Sakamoto & Haniu (1990) Sakamoto, H. & Haniu, H. 1990, ASME Journal of Fluids Engineering, 112, 386
  * Stellingwerf (1978) Stellingwerf, R. F. 1978, ApJ, 224, 953
  * Sterzik & Durisen (1995) Sterzik, M. F. & Durisen, R. H. 1995, A&A, 304, L9
  * Toalá et al. (2016) Toalá, J. A., Oskinova, L. M., González-Galán, A., et al. 2016, ApJ, 821, 79
  * Toropina et al. (2019) Toropina, O. D., Romanova, M. M., & Lovelace, R. V. E. 2019, MNRAS, 484, 1475
  * van Buren & McCray (1988) van Buren, D. & McCray, R. 1988, ApJ, 329, L93
  * van Marle et al. (2011) van Marle, A. J., Meliani, Z., Keppens, R., & Decin, L. 2011, ApJ, 734, L26
  * Vishniac (1994) Vishniac, E. T. 1994, ApJ, 428, 186
  * Wareing et al. (2007) Wareing, C. J., Zijlstra, A. A., & O’Brien, T. J. 2007, ApJ, 660, L129
  * Wilkin (1996) Wilkin, F. P. 1996, ApJ, 459, L31
  * Winkel et al. (2016) Winkel, B., Kerp, J., Flöer, L., et al. 2016, A&A, 585, A41
  * Yoon & Heinz (2017) Yoon, D. & Heinz, S. 2017, MNRAS, 464, 3297

Figure 14:  Molecular cloud proposed as the birthplace of UJT-1. (Left). Map
of the antenna temperature due to CO emission toward this proposed RATTS
corresponding to an LSR velocity of $-9$ km s-1 as detected in Dame’s survey
(Dame et al. 2001). Equatorial coordinates are used. The white contours mark
the location of the stellar wake, and the dashed red line shows the inferred
direction of the peculiar stellar velocity. A conspicuous gas cloud, outlined
by a dashed circle, is visible about 40 arcminutes away in the direction
opposite to the velocity. The 30 arcminutes horizontal scale bar is equivalent
to about 40 pc at the proposed distance of 4.5 kpc. (Right). Kinematic plot
showing the corresponding distance as a function of the LSR velocity within
the $-13$ to $-6$ km s-1 range (dashed lines) in which the cloud CO emission
is detected. A value close to 4 kpc is favored, in agreement with the distance
range to UJT-1 based on photometric and spectroscopic considerations.

##  Appendix A Origin of UJT-1

.

Our previous kinematic analysis enables us to guess the possible birthplace of
this protostar by extrapolating its peculiar velocity backward. We first did
this on the Dame carbon monoxide (CO) survey (Dame et al. 2001) in our Galaxy
(Fig. A14, left panel). Here, it is remarkable that a molecular cloud stands
well aligned less than one degree away in the opposite direction of motion of
the star, corrected for Galactic rotation. The cloud is located at J2000.0
coordinates $20^{h}32.7^{m}$ and $+44^{\circ}41^{\prime}$, at about 40
arcminutes from the proposed RATT, with an angular diameter of $\sim
10^{\prime}$. Moreover, its line emission is detected within the LSR range of
$-13$ to $-6$ km s-1 peaking at $-9$ km s-1. As shown in the right panel of
Fig. A14, this corresponds to a kinematic distance close to 4 kpc, which
agrees with the favored distance to UJT-1 of 4.5 kpc based on photometric and
spectroscopic arguments. When we propose this cloud association, we assume
that our simple Galactic rotation model is not very strongly affected by the
anomalous motions in the Cygnus superbubble region (Comerón et al. 1993). If
this is not the case, a closer distance estimate cannot be strictly excluded.

If the two coincidences in position and velocity are correct, they strongly
flag the cloud as the likely natal place from which the star was ejected.
Based on the peculiar velocity and the projected path, the travel time to its
current position is estimated as $\sim 1$ Myr. Additional confidence in the
LSR velocity and kinematic distance range comes from the Effelsberg-Bonn HI
Survey (EBHIS51) map in Fig. 13. Here, the LSR velocity of $-14$ km s-1 and
nearby channels depicts a clear minimum in the neutral hydrogen emission
extending along the western UJT-1 wake. This minimum seems to coincide with an
asymmetric shell shape in extended radio emission at 20 cm from the Canadian
Galactic Plane Survey (CGPS), mimicking an old supernova remnant that could
lie behind the persistent infrared dust emission in the wake of UJT-1.
However, the nature of this remnant-like feature remains to be confirmed.

##  Appendix B Calculating the cooling timescales

.

Let $E=\frac{1}{\gamma-1}n_{H}k_{B}T_{i}$ to be the local thermal energy
density, with $\gamma$ being the adiabatic index, $n_{H}$ the density number,
and $k_{B}$ the Boltzmann constant (the subscript $i$ means a shocked wind or
shocked ISM). Because no ionization field is expected in the vicinity of a
low-mass star like ours, the plasma is assumed to be in collisional ionization
equilibrium (CIE). Therefore, $\dot{E}=n^{2}_{H}f_{c}(T_{i})$ is the local
cooling rate of the gas, with $f_{c}$ being the commonly used cooling
function. Then, the cooling timescale is defined as

$\tau_{c,i}=\frac{E}{\dot{E}}=\frac{1}{\gamma-1}\frac{k_{B}T_{i}}{f_{c}(T_{i})n_{H}}.$
(2)

The cooling function may be approximated by a piecewise power law,
$f_{c}(T)=C(T)T^{\beta(T)}$. For monoatomic gas and a strong adiabatic shock,
we can write the cooling timescale in terms of the pre-shock conditions.
Assuming that the ISM and wind ram pressures are in balance,

$\tau_{c,ISM}=\frac{3}{8}\frac{k_{B}}{n_{H}c_{sISM}}\left(\frac{3\mu\mathscr{v}_{*}^{2}}{16k_{B}}\right)^{1-\beta_{sISM}},$
(3)
$\tau_{c,w}=\frac{3}{8}\frac{k_{B}}{n_{H}c_{sw}}\eta^{2}\left(\frac{3\mu\mathscr{v}_{w}^{2}}{16k_{B}}\right)^{1-\beta_{sw}},$
(4)

where $w$ and $ISM$ stand for wind and ambient medium, $sw$ and $sISM$ stand
for the shocked wind and ISM, and $C_{sISM}$ and $\beta_{sISM}$ are the
function of $T_{sISM}$, while $C_{sw}$ and $\beta_{sw}$ are the function of
$T_{sw}$. We adopt a pre-shock temperature of 3300 K, which corresponds to the
equilibrium temperature of the CIE cooling curve we used to estimate $C$ and
$\beta$ assuming solar abundances for H, He, and Z (Meyer et al. 2014).

##  Appendix C Calculating the instability growth timescale

Instability modes with wavelengths of about $\lambda$ are expected to grow to
finite amplitudes if the dynamical timescale
$\tau_{dyn}=R_{0}/\mathscr{v}_{*}\simeq 3.2\times 10^{11}$ s exceeds the
instability timescale for this wavelength. This is deduced from the solution
of the linear analysis and eigenvalue problem on each instability. For
instance, KH instabilities that are caused by the shear between the gases
coming from the stellar wind and the ISM have a typical growth timescale
$\tau_{\rm
KH}\sim(1+\chi)\chi^{1/2}\lambda(\mathscr{v}_{w}-\mathscr{v}_{*})^{-1}=(1+\chi)\chi^{-1/2}\tilde{\lambda}(\eta-1)^{-1}\tau_{dyn}$,
where $\chi=\rho_{sISM}/\rho_{sw}\sim 2.6\times 10^{3}$ is the density jump in
the contact discontinuity, $\lambda$ is the perturbation wavelength, and the
tilde is used for dimensionless lengths relative to $R_{0}$. Substitution of
numerical values gives $\tau_{\rm KH}\sim 5.7\tilde{\lambda}\tau_{dyn}$, so
that only $\tilde{\lambda}<1.8\times 10^{-1}$ can grow fast enough. This
excludes the case of $\lambda\sim R_{0}$, but allows perturbations of about
the layer width at the apex, $h_{0}$, where the instabilities first appear
(Comerón & Kaper 1998), because $\tilde{h_{0}}\sim\mathcal{M}^{-2}\sim
10^{-2}$, where $\mathcal{M}\sim 10$ is the Mach number for the star. However,
RT instabilities, which take place when a lighter fluid supports a heavy
fluid, apparently do not occur at first because $\chi\gg 1$ and the effective
gravity $g_{eff}$ is due to centrifugal force and directed outward.
Nevertheless, if the shocked layer were bent inward by a perturbation
triggered by a KH instability, the curvature would induce a centrifugal
$g_{eff}$ directed inward, with the dense shocked ISM over the less dense
shocked wind. Then the corresponding growth timescale would be $\tau_{\rm
RT}\sim\left[2\pi\lambda(\chi+1)(\chi-1)^{-1}g_{eff}^{-1}\right]^{1/2}$. We
again assume $\lambda\sim h_{0}$ and estimate the effective gravity by
considering that the ram pressure
$\rho_{a}\mathscr{v}_{*}^{2}\sim\rho_{ISM}g_{eff}R_{c}=\mathcal{M}^{2}\rho_{a}g_{eff}R_{c}$,
with $R_{c}$ being the curvature radius. As
$\tilde{h_{0}}\sim\mathcal{M}^{-2}$ and $\chi\gg 1$, then we can write
$\tau_{\rm RT}\sim\tau_{dyn}(2\pi\tilde{R_{c}})^{1/2}$, so that the modes with
wavelength $\lambda\sim h_{0}$ can grow enough due to RT instabilities if
$\tilde{R_{c}}<(2\pi)^{-1/2}$, which seems reasonable at first glance.
However, as KH instabilities of about $R_{0}$ are not expected to grow, it is
difficult to obtain such a large $R_{c}$. On the other hand, NTS instabilities
that arise from the shear in shock-bounded slabs produced by deviations from
equilibrium might trigger these waves. Their growth timescale is $\tau_{\rm
NTS}\sim 2\pi(c_{s}k(kl)^{1/2})^{-1}$, with $l$ being the slab displacement
and $k$ the wavenumber (Vishniac 1994). As the most unstable $\lambda$
corresponds to $kl\sim 1$, then $\tau_{\rm
NTS}\sim\lambda/c_{s}=\tilde{\lambda}\mathcal{M}\tau_{dyn}$ (which, as may be
seen, is one order of magnitude greater than $\tau_{\rm KH}$), from which we
infer that this mechanism might maintain perturbations with a wavelength
$\tilde{\lambda}<\mathcal{M}^{-1}$, which again excludes $\lambda\sim R_{0}$,
but allows instabilities to grow near the apex. TA instabilities occur as a
consequence of the acceleration in the flow normal to the shock due to its
curvature. Their growth timescale is $\tau_{\rm TA}\sim(2\pi\lambda
l_{c}^{-1})^{1/2}R_{0}\mathscr{v}_{w}^{-1}=(2\pi\tilde{\lambda}\tilde{l_{c}}^{-1})^{1/2}\tau_{dyn}\eta^{-1}$,
where $l_{c}$ is the distance along the bow shock measured from the stagnation
point at the apex (Dgani et al. 1996). Thus, TA instabilities can grow fast
enough in our case if $\tilde{\lambda}<\eta^{2}/(2\pi)\tilde{l_{c}}$, so that
they are more efficient for larger $l_{c}$, and they can hardly explain the
first instabilities in the vicinity of the apex ($l_{c}\sim 0$). Therefore, TA
instabilities at best amplify the distortions generated by NTS and/or KH
instabilities as they are carried on to the wings of the bow shock. This may
be shown bettern when we consider $\tau_{\rm TA}/\tau_{\rm
NTS}\sim(2\pi{\tilde{\lambda}}^{-1}{\tilde{l_{c}}}^{-1})^{1/2}\eta^{-1}\mathcal{M}^{-1}$,
which indicates that TA instabilities grow faster than NTS instabilities only
for $\tilde{l_{c}}>2\pi\eta^{-2}\sim 10^{-1}$, so that the latter are dominant
in the apex region. It is noteworthy that although NTS and TA instabilities
were defined under assumptions different from the bow-shock configuration,
dedicated numerical simulations show that these instabilities can also cause
wiggles in the bow-shock shell (Blondin & Koerwer 1998). All these conclusions
still hold for the case of the shorter distance $d=2$ kpc because the order of
magnitude of the involved timescales does not change.

##  Appendix D Vortex shedding as a possible cause of undulations in the wake

Vortex shedding with frequency $f$ from a sphere of diameter $D$ moving at
velocity $v$ in a fluid with density $\rho$ and dynamical viscosity $\mu_{f}$
is fully characterized by the Reynolds number $Re=\rho vD\mu_{f}^{-1}$,
together with the Strouhal number $St=fDv^{-1}$ (Landau & M. 1987). In the
laboratory (Sakamoto & Haniu 1990), vortex shedding appears in a range of
$Re=3\times 10^{2}$-$3.7\times 10^{5}$, with $St\sim 0.2$ for the highest
$Re$, although this last value can be higher ($\gtrsim 0.6$) in the hypersonic
regime (Florsheim et al. 1965). In our case, $f$ was estimated above to be
$\sim 10^{-12}$ Hz, while the viscosity $\mu_{f}$ can be obtained from the
Chapman-Eskong approximation (Bird 1994). Taking $R_{0}$ as the characteristic
length, the substitution of the estimated parameters renders $Re\sim 4.1\times
10^{4}$ and $St\sim 0.45$. While $Re$ is compatible with the hypersonic
expectations, the $St$ value is only slightly lower. Similar results are
obtained when the distance to the star is $d=2$ kpc.