# The surprising evolution of the shadow on the TW Hya disk111Based on
observations made with the NASA/ESA Hubble Space Telescope, obtained at the
Space Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.
These observations are associated with programs #16228 and #13573.

John Debes AURA for ESA, Space Telescope Science Institute, 3700 San Martin
Dr., Baltimore, MD 21218, USA  Rebecca Nealon Department of Physics,
University of Warwick, Coventry CV4 7AL, UK  Centre for Exoplanets and
Habitability, University of Warwick, Coventry CV4 7AL, UK  Richard Alexander
School of Physics and Astronomy, University of Leicester, University Rd.,
Leicester LE1 7RH, UK  Alycia J. Weinberger Earth and Planets Laboratory,
Carnegie Institution for Science, 5241 Broad Branch Rd NW, Washington, DC
20015  Schuyler Grace Wolff Department of Astronomy, The University of
Arizona, 933 N Cherry Ave, Tucson, AZ, 85750, USA  Dean Hines Space Telescope
Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA  Joel Kastner
School of Physics & Astronomy, Rochester Institute of Technology, 1 Lomb
Memorial Dr., Rochester, NY 14623, USA  Hannah Jang-Condell NASA Headquarters,
300 Hidden Figures Way SW, Washington, DC 20546, USA  Christophe Pinte School
of Physics and Astronomy, Monash University, Clayton Vic 8300, Australia
Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France  Peter Plavchan
Department of Physics & Astronomy, George Mason University, 4400 University
Dr., Fairfax, VA 22030  Laurent Pueyo Space Telescope Science Institute, 3700
San Martin Dr., Baltimore, MD 21218, USA

###### Abstract

We report new total intensity visible light high contrast imaging of the TW
Hya disk taken with the Space Telescope Imaging Spectrograph (STIS) on the
Hubble Space Telescope (HST). This represents the first published images of
the disk with STIS since 2016, when a moving shadow on the disk surface was
reported. We continue to see the shadow moving in a counter-clockwise fashion,
but in these new images the shadow has evolved into two separate shadows,
implying a change in behavior for the occulting structure. Based on radiative
transfer models of optically thick disk structures casting shadows, we infer
that a plausible explanation for the change is that there are now two
misaligned components of the inner disk. The first of these disks is located
between 5-6au with an inclination of 5.5°and PA of 170°, the second between
6-7au with and inclination of 7°and PA of 50°. Finally, we speculate on the
implications of the new shadow structure and determine that additional
observations are needed to disentangle the nature of TW Hya’s inner disk
architecture.

Circumstellar disks (235), Exoplanet formation (492), Coronagraphic imaging
(313), Hubble Space Telescope (761), Protoplanetary disks (1300)

##  1 Introduction

TW Hya is a nearby, well-studied face-on protoplanetary disk (d=60
pc;i$\sim$5° Gaia Collaboration et al., 2021; Qi et al., 2004; Huang et al.,
2018). It has been observed from X-ray to radio wavelengths. Estimates of TW
Hya’s stellar mass generally lie in the range 0.6-0.8 M⊙, with recent
(dynamical) estimates favoring the upper end of this range (e.g. Teague et
al., 2022). Given its age ($\sim$10 Myr), this makes the TW Hya disk an
interesting test case for planet formation in a large, relatively isolated,
gas-rich disk orbiting a star somewhat less massive than the young Sun, but
more massive than transiting young M dwarf systems such as AU Mic. (Teague et
al., 2019; Sokal et al., 2018). HST has imaged the disk numerous times,
including with WFPC2 (IWA$\approx$0.7″; Krist et al., 2000), NICMOS Near-IR
(NIR) F110W and F160W (IWA$\approx$0.4″; Weinberger et al., 2002), and STIS
visible light imaging/spectroscopy (IWA$\approx$0.2″; Roberge et al., 2005;
Debes et al., 2017, hereafter D17). The disk is polarized, allowing visible
and NIR images (Akiyama et al., 2015; van Boekel et al., 2017; Poteet et al.,
2018). In the sub-mm, dust continuum emission extends from 3-$\sim$100
au(Andrews et al., 2016; Huang et al., 2018; Ilee et al., 2022), with multiple
gaps or rings present that are close in radius to structures seen in scattered
light and consistent with active planet formation (Jang-Condell & Turner,
2012; Andrews et al., 2016; Ueda et al., 2020; Macías et al., 2021, and
references therein). CO gas emission extends to $\sim$200 au, and scattered
light from small dust extends as far out as 452 au(Huang et al., 2018, D17).

The TW Hya disk has an azimuthal surface brightness asymmetry between 30-150
au (Roberge et al., 2005; Debes et al., 2013). The asymmetry has been
attributed to inclination effects for flared disks (Debes et al., 2013), the
forward scattering properties of dust (Roberge et al., 2005), or the presence
of a warp (Roberge et al., 2005; Rosenfeld et al., 2012). Images of the
asymmetry with STIS and NICMOS over six epochs between 0.6-2$\mu m$ show that
the asymmetry is not stationary; the apparent changes in disk illumination can
be modeled with a rotating shadow with a constant angular velocity of 22.7∘/yr
(P=15.9 yr). Assuming the shadow arises from an obstructing feature moving
with Keplerian velocity constrains the location of the originating structure
to $<5.6(\frac{M_{\star}}{0.7M_{\odot}})$ au (D17). ALMA CO maps of the outer
disk also show the shadow (Teague et al., 2022). The position of the shadow in
CO taken in 2019 is smaller than expected based on the predictions of D17.

The TW Hya shadow is one of a handful of other disks that show evidence of
inner occulting structures, often interpreted as mis-aligned inner disks
casting narrow or broad shadows depending on the inclination difference with
the outer disk (e.g. PDS 66, HD 34700, HD 16192, HD 142527, HD 139614, and
RXJ1604.3-2130A Wolff et al., 2016; Uyama et al., 2020; Bertrang et al., 2018;
Marino et al., 2015; Muro-Arena et al., 2020; Pinilla et al., 2018). Some of
these disks also show evidence for time variability, but robustly sampling the
variability timescales is difficult for high contrast imaging targets. The
robustness of the shadow detection and coverage of its characteristic
timescale of motion make the shadow on TW Hya’s disk one of the best
characterized shadows to date that may be caused by active planet formation.

Nealon et al. (2019) investigated the behavior of shadows due to misaligned
disks that precess and warps associated with planetary companions orbiting the
host star, in order to understand TW Hya’s shadow. They found that massive
planets with inclinations relative to the outer disk of a few degrees or more
were sufficient to create shadows of the right amplitude in the outer disk
either through a warp exterior to the planet’s orbit or from a misaligned disk
interior to the planet. This work motivated a multi-Cycle HST program with the
STIS coronagraph in order to disentangle these two origins. TW Hya has only
been observed for about half the observed periodicity, and a misaligned,
precessing disk will continue to show asymmetries in the disk for a whole
precession period, while features associated with a misaligned planet will
disappear for half an orbital period. It is possible that both structures are
present and contributing to the shadow behavior.

We report here the first epoch results from our monitoring program of TW Hya’s
disk to further track the evolution of the asymmetry and putative shadow. In
§2 we describe the new STIS observations. In §3 we report our findings of the
latest epoch which show a marked departure in the behavior of the shadow. In
§4 we discuss the possible explanation of this evolution. Finally we present
our Conclusions in §5.

Figure 1: Logarithmic images of TW Hya from various epochs. North is up and
East to the left of the images. (upper left) Image of the first visit of
Program 16228 taken on 07 June 2021. Black areas represent missing data due to
diffraction spikes and occulting masks. (upper right) The second visit of
Program 16228, taken at a different spacecraft orientation. (lower left) Fully
combined images of TW Hya from 2016 (lower right) Combination of Visit 1 and
Visit 2 from Program 16228. Figure 2: STIS visible light images of the TW Hya
disk in 2016 and 2021. (top) Each panel is labeled with inidividual PSF
subtracted visits of TW Hya taken in 2021 behind BAR5 and WEDGEA1.0 and
multiplied by $R^{2}$ to highlight the rings and gaps present in the system.
Diffraction spikes and areas blocked by the occulters are masked out. We
recover the overall nature of the disk irrespective of the spacecraft
orientation. (bottom) Comparison of the disk between 2016 and 2021. In 2021,
the smooth shadow features (seen at a PA of $\sim$50$\arcdeg$ in 2016 and
marked as ”A”) have split into two shadow features marked by the ”B” at 49
$\arcdeg$ and ”C” at 175 $\arcdeg$.

##  2 Observations

TW Hya was observed on 7 June 2021 with the STIS Instrument (Woodgate et al.,
1998), using the CCD and 50CORON aperture as part of GO#16228 (PI: Debes). The
same PSF reference for TW Hya, HD 85512, was used as in Roberge et al. (2005)
and D17. TW Hya was observed for two consecutive orbits with spacecraft
orientations separated by 17$\arcdeg$, followed by an orbit on HD 85512, and a
final orbit on TW Hya separated by 32$\arcdeg$ from the first TW Hya orbit.
Each star was first placed behind the BAR5 aperture location which has a half-
width of 0$\farcs$15 with inner working angles of $\sim$0$\farcs$2 (Schneider
et al., 2017; Debes et al., 2019). For TW Hya, 5$\times$110 s exposures were
obtained per orbit behind BAR5, with sufficiently short exposure times chosen
to avoid saturating the detector at the mask edge. The PSF reference was
dithered perpendicular to BAR5 by $\pm$12.5 mas to account for the 0.25 pixel
target acquisition non-repeatability (Schneider et al., 2017). We obtained
5$\times$11.2 s exposures at each dither location. Next, each star was placed
behind the WEDGEA1.0 occulter to ensure sufficient signal-to-noise (S/N) in
the outer disk without saturating the detector. The WEDGEA1.0 exposures were
3$\times$425 s per orbit with TW Hya and 22$\times$32 s exposures for HD
85512. The total integration time for TW Hya was 1100 s behind BAR5 and 2550 s
behind WEDGEA1.0. This procedure replicates the other STIS images obtained in
2015 and 2016 (D17).

Unfortunately, a delay in guide star acquisition rendered the final orbit of
TW Hya unusuable, since the star was not placed behind the proper apertures.
The purpose of the multiple orbits of TW Hya are primarily to ensure nearly
360∘ coverage on the disk exterior to 0$\farcs$5\. The result of the failed
orbit is a modest decrease in disk S/N and a smaller angular coverage at small
angular radii as seen in Figure 1. However, as we show below, this did not
significantly impact our results when comparing to previous epochs at radii
$>$40 au (0.67$\arcsec$).

We followed the PSF subtraction laid out in D17, which we briefly review.
First, the stellar center was determined for both target and PSF reference. We
then iteratively registered and scaled the PSF reference, minimizing
subtraction residuals within a mask that included the diffraction spikes of
the star and subtracted off the PSF reference. We note that the ratio of TW
Hya’s flux to HD 85512 across the STIS CCD bandpass was higher in 2021 between
0.2-1.0$\micron$ by 24-31%. Compared to 2015 and 2016 where the best fit ratio
of TW Hya’s flux to HD 85512 was 0.052 and 0.055, we find that the best
scaling between the two stars is now 0.0683.

TW Hya has been known to be variable by $\sim$20-30% in the visible on periods
close to $\sim$3 days (Mekkaden, 1998). To investigate whether this magnitude
of brightening might be expected in the STIS bandpass, we downloaded existing
contemporaneous STIS/G430L and STIS/G750L spectra taken with the 52x2 slit of
TW Hya available in MAST and estimated the count rate on the STIS detector for
the 50CORON mode for each spectrum using the STIS ETC. The spectra came from
four epochs (28 Jan 2010, 4 Feb 2010, 28 May 2010, and 18 Apr 2015) as part of
programs 11608 (PI: Calvet) and 13775 (PI: Espaillat). The varying flux in the
blue continuum and in emission line intensity of the different spectra account
for peak-to-peak changes in predicted count rates on the STIS detector across
the four epochs of 22%, in line with the variability observed. Additionally,
we inspected the acquisition images for TW Hya (which used the F28X50LP
filter) between 2016 and 2021. Even though the filter throughput for the
acquisition images is for flux beyond 5000 Å, aperture photometry of the
acquisition images showed a 22% increase in TW Hya’s brightness between the
two epochs. The change in SED for TW Hya adds modest variability in how well
the target PSF reference matches in color to TW Hya, which could result in
more or less PSF subtraction residuals.

With the scalings, we recover the disk again, seen in Figure 1. We show
reductions of each individual visit to demonstrate the impact of PSF
subtraction residuals. Residual features are present only in the detector
frame of the observations and will appear to rotate when in a sky oriented
frame. While some modest features appear to be due to PSF residuals, they are
at the level of our estimated noise per pixel and do not impact our
measurements signficantly for the analysis described in §3.

As noted in D17 and van Boekel et al. (2017, hereafter VB17), when one removes
the 1/R2 dependence of disk illumination from the star, a gap is present at
$\sim$80 au, with an inner ring peaking at $\sim$30 au, with another gap
interior to that as shown in the 2016 image in Figure 1. There is an azimuthal
surface brightness asymmetry, which we interpret as due to a shadow (D17). In
past epochs, this shadow asymmetry was well described by a smooth sinusoidal
curve. The position angles (PA) of the peak and trough of the asymmetry have,
to date, moved with constant angular velocity across several wavelengths in
the visible and NIR. Using the 15.9 yr predicted period from D17 and setting
$t_{\mathrm{o}}$ to equal the first observation of the TW Hya disk, we find
that roughly half of the predicted orbit for the shadow has not been covered,
and the June 2021 observations mark a similar phase of the orbit to archival
NICMOS polarimetric observations (Poteet et al., 2018) and archival medium
band NICMOS total intensity observations of the disk (D17) taken in 2004 and
2005, both of which are of lower S/N compared to the observations with STIS.

We created an azimuthal average surface brightness profile of the 2021
dataset, converted it into physical flux units per pixel, and compared the
resulting fluxes with azimuthal averages of the 2016 observations of the TW
Hya disk. Due to the broad bandpass of the 50CORON mode, we estimate the
conversion from counts s-1 on the detector to surface brightness by using the
stsynphot package and the average spectrum of TW Hya obtained with STIS. We
find that the conversion is 1 mJy arcsec-2=3.54 DN/s. We find that the flux
profile of TW Hya is slightly different in 2021 compared to 2016, primarily
fainter by 20% interior to 1$\arcsec$ as well as beyond 1.5$\arcsec$.

##  3 Azimuthal Surface Brightness Evolution

Unlike the observations taken in 2016, the inner ring shows two narrow dark
lanes at PA of $\sim$50∘ and $\sim$180∘. The peak of the azimuthal asymmetry
has moved counter-clockwise toward due North as predicted by the shadow
rotation model presented in D17. To quantify these changes, we replicate the
methods of D17 by measuring the azimuthal surface brightness (ASB) of the disk
as a function of PA and radius for the 2021 epoch and compare those to the
previously published profiles in 2016.

To obtain azimuthal asymmetries, we divided each ASB profile for a given
radius by the mean SB at that radius to allow for unbiased radial averaging.
We chose angular radii that matched the measurements made in D17, and
converted to physical differences with the new Gaia EDR3 parallax (Gaia
Collaboration et al., 2016, 2021): seven radial locations centered at 0.46″
(28 au), 0.66″ (40 au), 0.89″ (54 au), 1.14″ (68 au), 1.47″ (88 au), 1.83″
(110 au), and 2.36″ (142 au), with widths of 0.2″ (12 au), 0.2″ (12 au), 0.25″
(15 au), 0.25″ (15 au), 0.41″ (25 au), 0.30″ (18 au), and 0.56″ (34 au)
respectively. The locations are within and near the two gaps seen in scattered
light at 22 and 88 au with sufficient S/N. We estimated the root-mean square
uncertainty per azimuthal bin by calculating both the average uncertainty per
pixel as estimated by the total counts within each pixel or by the standard
deviation of counts within an azimuthal bin, whichever was larger. At these
radii, the profiles in 2016 beyond 40 au have a median S/N of 19 per azimuthal
bin, while for 2021 the median S/N per azimuthal bin is 18.

Figure 3 shows the results for the 2016 and 2021 epochs. The simple cosine fit
used in D17 is no longer sufficient to explain the observed ASB profile in
2021, thanks to the apparent presence of two shadow features. Figure 3 shows
that the two shadows diminish in depth at larger radii, such that they are no
longer distinguishable as two components and that the peak surface brightness
increases in PA with radius. We thus fit a cosine at 141 au and the PA of the
peak surface brightness is 330${}^{\circ}\pm 10^{\circ}$, in tension ($\sim
2.3\sigma$) with the 2005 measurement from the NICMOS medium band data at a
similar phase of the predicted motion reported in D17. We note that the 2004
NICMOS total intensity data derived from the POL filters in D17 and P19 showed
a similar lower than expected PA for the given period.

There are three possibilities that can explain the new observations. First,
there are new features caused by the recent appearance of additional shadowing
structures unrelated to what has been seen before, such as new dust clumps.
Secondly, warps due to companions that cast shadows can appear when the planet
is on the near side of the disk from the perspective of the observer (Nealon
et al., 2019). This is possible since we are seeing the phenomena at a PA that
is close to the major axis of the disk at 152$\arcdeg$ as measured by CO data
(Huang et al., 2018). Thirdly, it is possible that these two shadows arise
from two mutually inclined disks that have precession timescales which are
slightly different and casting overlapping shadows, thus mimicing a single
shadow over the previously observed epochs.

The ASB profiles can be fit by combinations of shadows and cosine functions.
At present it is difficult to find a unique solution between these
possibilities. We therefore focus on the simplest case–that there are two
shadows that previously overlapped but are moving at different angular
velocities and partially separated in 2021. Thus, we can approximate the shape
of the shadows as a linear combination of Gaussian curves.

In order to fit the two shadow features simultaneously, we use CURVEFIT.pro,
an IDL nonlinear least squares fitting routine assuming two Gaussian shadows,
with constant offset ($C$), shadow depth ($D_{\mathrm{i}}$), shadow PA
($\theta_{\mathrm{o},i}$), and shadow angular width ($\sigma_{\mathrm{i}}$):

$F(\theta)=C+\sum_{i=1}^{2}{}D_{\mathrm{i}}\exp{\frac{-(\theta-\theta_{\mathrm{o},i})^{2}}{2\sigma_{\mathrm{i}}^{2}}}$
(1)

where we convert the shadow width to an angular full-width at half maximum
$FWHM_{\mathrm{i}}\approx 2.3548\sigma_{\mathrm{i}}$. The fits are shown in
Figure 3 and show no significant residuals.

Figure 3: Comparison of the 2016 and 2021 epochs of TW Hya’s azimuthal surface
brightness distribution as a function of radius. On the left, the 2016 epoch
shows that the azimuthal brightness is well defined by a cosine or sine curve
with a clear peak and trough position. The 2021 epoch to the right shows a
departure from this behavior. The presence of two quasi-Gaussian shadows are
present. The red lines show cosine fits to the data taken in 2016 from D17 and
two-Gaussian fits to the data taken in 2021. Figure 4: Comparison of the
azimuthally averaged surface brightness profiles of TW Hya from 2016 and 2021.
At small radii and larger radii, the surface brightness has slight differences
of 20%

Figure 5 shows the PA, widths, and depths of the shadows as a function of
radius. Shadow A has a median PA of 49$\arcdeg$, while Shadow B has a median
PA of 175$\arcdeg$. The shadows show a slight increase in PA in the
counterclockwise direction as radius increases, which may be indicative of
warping in the disk surface.

Figure 5: Fit parameters for shadow A and B as a function of radius. Shadow A
is represented by black squares while shadow B is represented by red squares.
(top) Relative PA from the median. (middle) The measured FWHMs of the two
shadows do not show significant changes with radius. Note that at the largest
radius, the uncertainty in FWHM is large since shadow position and width are
degenerate when the shadows overlap. The shadows decrease in depth as a
function of radius.

The FWHM of the shadows is broad: the median FWHM of shadow A is 87$\arcdeg$,
shadow B is 93$\arcdeg$. These widths can explain why the shadow might have
appeared to be singular previously if the separation of the two shadows was
less than their FWHMs. Finally, the depth of the shadow changes with radius,
starting at depths of near 80% and decreasing to 20% by 140 au.

Figure 6: Schematic of our proposed model to explain the evolution of the TW
Hya shadow. We consider two mutually inclined rings between 5-6 au and 6-7 au
with slightly different inclinations relative to the outer disk as an origin
to the shadows we observe in 2021. Figure 7: Comparison between azimuthal
profiles at different radii taken in 2021 and different radiative transfer
models of shadows from inclined rings. The black squares are the data with
error bars, while the orange dashed line is the profile from a single inclined
ring casting a shadow. The red solid line is the profile derived from two
inclined rings but with differing PAs. The two-ring model better reproduces
the observed profiles than a single ring model.

##  4 Implications for the Origin of the TW Hya Shadows

We review the following characteristics and evolution of the disk shadow in
order to provide a summary of the current observational constraints. We will
restrict ourselves to those behaviors that are seen in the STIS data, since
this provides the set of constraints that cannot be explained by changes in
behavior of the shadow with wavelength, instrumental setup/systematics, or
disk polarization.

From 2000-2016 the shadow covered 180$\arcdeg$ of azimuth and the motion of
both edges of the shadow were consistent with constant angular motion. The
shadow depth was constant over this time period and with radius. The PA of the
shadow at r$<$50 au did not always match that of outer parts of the disk. The
disk shadow detected prior to 2016 may have evolved into two shadows by 2021.
The behavior of the shadow between 2000, 2015, and 2016 suggested that the
shadow rotated counter-clockwise with a period of $\sim$15.9 years. The
behavior of the shadow in 2021 may mean that either rotation or oscillatory
behavior of the shadow is possible, but the overall period still seems to be
$\sim$15 yr, limited somewhat by sparse time sampling of the shadow’s motion
with STIS-only epochs. In 2021, there is a trend where the depth of the
shadows changes quasi-monotonically with radius, with the shadow becoming less
pronounced further out in the disk. This behavior is not as pronounced in
other STIS epochs.

In order to put a quasi-physical explanation to this new shadow behavior, we
begin in 4.1 to model the disk without any shadowing structures as a baseline
disk model to compare to two different scenarios, a single inclined ring
casting a shadow and two mutually inclined rings casting shadows. In 4.2 we
attempt to fit the data with one inclined ring. This is because a tilted ring
is physically motivated (e.g. D17), conserves angular momentum, in some cases
has been observed (Marino et al., 2015) and has been inferred in other disks
(Muro-Arena et al., 2020). We then use the rapid radiative transfer code
MCFOST (Pinte et al., 2006, 2009) to model the shadow cast by an inclined ring
to try and change the location, shape and inclination of this ring to match
the flux profiles in Figure 3. Ultimately, we show the best fitting inclined
inner ring cannot fit the observed ASB and investigate the impact of adding a
second inclined ring in 4.3. There we find a good match to both the SED of the
disk and the visible scattered light images, which suggests that this scenario
is the most likely explanation for the behavior we observe.

Figure 8: (top) Black points represent the SED of TW Hya. The red curve is the
model SED from our fiducial radiative transfer model. (bottom) Comparisons
between the fiducial model SED and those of the single and double ring models.
All models are similar to within 5-10% of the fiducial model.

###  4.1 A fiducial Model for the TW Hya disk

The first step is to create a disk model that is broadly consistent with both
the visible total intensity scattered light morphology and the spectral energy
distribution measured for TW Hya without any shadows. To that end, we
recreated the model proposed by VB17 with MCFOST. VB17 was focused on
simultaneously fitting the polarized intensity radial surface brightness
profile of the disk near-IR wavelengths with the TW Hya SED. This model was
also derived from work done by Menu et al. (2014), which was a model focused
on self-consistently fitting multi-wavelength interferometry data with the TW
Hya SED. Both models relied on a 1-D axisymmetric disk structure with a
vertical scale height dictated by hydrostatic equilibrium. Both models assumed
that the inner gas disk was truncated with a rounded inner taper and that the
largest dust grains were decoupled from the gas radial density profile. VB17
additionally added an outer gas density taper at a larger radius than Menu et
al. (2014), used the ALMA 870 µm radial surface brightness profile (Andrews et
al., 2016) as a proxy for the large dust radial profile, and added gaps in the
form of density decrements at particular radii to reproduce the polarized
intensity radial scattered light surface brightness profiles observed with
SPHERE. For the purposes of our models, we assumed that the outer disk has an
inclination of 0${\arcdeg}$ and a position angle on the sky of 155°. While the
true inclination of the outer disk is likely closer to 5$\arcdeg$, the impact
to both the scattered light images and SED are negligible.

While we consider a very similar disk structure, we make some alternate
choices that appear to retain the basic observational features of the TW Hya
disk, while making the model more amenable to calculating the radiative
transfer of tilted inner rings. We are not trying to make a disk model that
matches all available observations but a plausible model that can be used to
investigate the impact of tilted inner rings on the observed surface
brightness. We leave a deeper investigation into the degeneracies and
constraints imposed by disk shadows for later work. In general, we find decent
agreement with those of VB17, though we point out differences where they
occur.

####  4.1.1 Radial Surface Density of Gas and Small Dust

The previous two models of the TW Hya disk calculated an analytical radial
surface density profile for the gas while assuming a gas to dust ratio for the
small dust grains and calculated a self-consistent scale height for the disk
as a function of radius. Our approach is to calculate a fully analytical gas
density distribution with the following form:

$\rho(r,z)=\Gamma_{in}(r)\frac{\Sigma(r)}{\sqrt{2\pi}H(r)}f(r)\Gamma_{out}\int_{-\infty}^{\infty}\exp{\left(-\frac{z^{2}}{2H(r)^{2}}\right)}dz$
(2)

where $\Gamma_{in}(r)$ is the interior density taper, $\Sigma(r)$ is the
surface density profile, $f(r)$ is a multiplicative combination of density
gaps to mimic the gaps observed in various scattered light datasets,
$\Gamma_{out}$ is the exterior density taper, and $H(r)$ is the vertical scale
height of the gas (VB17, and references therein).

The interior taper replicates the format in VB17 and is:

$\Gamma_{in}(r<r_{in})=\exp{\left[\left(\frac{1-r/r_{in}}{w}\right)^{3}\right]}.$
(3)

The surface density profile is a power-law:

$\Sigma(r)=\left(\frac{r}{r_{in}}\right)^{-p}$ (4)

where we choose p=-0.7 and $r_{in}$ is the location where the interior taper
takes effect.

The scattered light profile of TW Hya clearly shows the evidence of gaps or
other perturbations in the upper layers of the disk, some of which correspond
to features also seen in the sub-mm where the disk is optically thin beyond a
few 10s of au(Andrews et al., 2016; Debes et al., 2013, VB17). We reproduce
these features by including Gaussian density perturbations following VB17:

$\displaystyle f_{i}=$ $\displaystyle
1-d_{i}\exp{\left[-\frac{(r-r_{c,i})^{2}}{2\sigma_{in,i}^{2}}\right]}$
$\displaystyle\mathrm{for}\ r<r_{c,i}$ (5) $\displaystyle f_{i}=$
$\displaystyle
1-d_{i}\exp{\left[-\frac{(r-r_{c,i})^{2}}{2\sigma_{out,i}^{2}}\right]}$
$\displaystyle\mathrm{for}\ r>r_{c,i}$ (6)

While we assume the perturbations are in density, we do not include similar
perturbations in scale height which might also occur in a system undergoing
gap formation (Bi et al., 2021).

We include a taper to the outer part of the gas density to match the sharp
drop in surface brightness seen in scattered light and in deep observations of
the outer gas disk in the sub-mm (D17, Ilee et al., 2022). We assume the
following functional form for the taper that occurs at $r>r_{out}$:

$\Gamma_{out}=\exp{\left[\frac{-(r-r_{out})^{2}}{2\sigma_{out}}\right]}\mathrm{for}\
r>r_{out}\\\ $ (8)

Finally, we assume an analytical scale height such that $\frac{H}{r}$=0.05 at
0.5 au and that the disk has a flared geometry:

$H(r)=0.05\left(\frac{r}{0.4au}\right)^{1+q}$ (9)

where $q=0.15$, shallower than the flaring expected for a purely isothermal
disk (Chiang & Goldreich, 1997). We choose this analytical form to best match
the SED and the surface brightness of the observed disk. To realize the
density distribution we create a cylindrical grid in the density with
4.8$\times 10^{5}$ cells logarithmically spaced in radius.

####  4.1.2 Stellar Parameters

The fundamental stellar parameters of TW Hya cover a wide range in mass and
Teff, primarily because of the intrinsic difficulty of determining the
parameters of young stars. For the purposes of the fiducial model, we assume
similar values to those assumed in VB17, accounting for the larger distance to
TW Hya determined by Gaia DR3(Gaia Collaboration et al., 2021). In MCFOST, we
chose Teff=4000 K and $\log g$=4.0, assuming a mass of 0.6 $M_{\odot}$. These
parameters are quite close to those determined by recent studies of TW Hya’s
temperature and gravity (Sokal et al., 2018). They are also consistent with
the inferred M${}_{\star}sin~{}i$ from CO measurements of Keplerian velocities
(Huang et al., 2018; Teague et al., 2018).

Table 1: Table of parameters used to describe the density distribution of gas in the TW Hya disk Parameter | VB17 | Fiducial | One Ring | Two Rings
---|---|---|---|---
p | -0.75 | -0.70 | -0.70 | -0.70
w | 0.45 au | 0.50 au | 0.50 au | 0.50 au
rin | 2.7 au | 3 au | 3 au | 3 au
q | …aaq was not a free parameter in VB17’s model but was calculated via hydrostatic equilibrium | 0.15 | 0.15 | 0.15
rout | 104 au | 100 au | 100 au | 100 au
$\sigma_{out}$ | 50 au | 50 au | 50 au | 50 au
Table 2:  Table of parameters used to describe the density gaps in the TW Hya disk. Parameter | VB17aacorrected for Gaia DR3 distance | Fiducial | One Ring | Two Rings
---|---|---|---|---
d1 | 0.44 | 0.44 | 0.44 | 0.42
rc,1 | 94 au | 85 au | 85 au | 85 au
$\sigma_{in,1}$ | 16.7 au | 12 au | 12 au | 10 au
$\sigma_{out,1}$ | 17.8 au | 30 au | 30 au | 30 au
d2 | 0.61 | 0.65 | 0.65 | 0.9
rc,2 | 23 au | 24 au | 24 au | 24 au
$\sigma_{in,2}$ | 4.2 au | 8 au | 8 au | 10 au
$\sigma_{out,2}$ | 20 au | 8 au | 8 au | 8 au
d3 | 0.81 | 0.81 | 0.81 | 0.92
rc,3 | 6.7 au | 6.7 au | 6.7 au | 6.7 au
$\sigma_{in,3}$ | 3.3 au | 3.3 au | 3.3 au | 3.3 au
$\sigma_{out,3}$ | 14.4 au | 18.7 au | 18.7 au | 18 au

####  4.1.3 UV Excess

For completeness, we include a UV excess to TW Hya that fits the visible
observed photometry. In reality, TW Hya’s UV flux can vary by factors of two
or more. We verified that the addition or subtraction of the UV excess in
MCFOST has a minimal impact on the predicted scattered light surface
brightness and the SED.

Figure 9: Comparison between the 2021 STIS image of the TW Hya disk (upper
left) compared to different models, including a fiducial model with no tilted
inner rings, a singular tilted inner ring, and two tilted inner rings.

####  4.1.4 Grain Properties

We initially replicated the composition of dust used in VB17, namely using a
mixture of 80% amorphous magnesium iron olivine-style silicates (MgFeSiO4)
(Dorschner et al., 1995) and 20% amorphous carbon (Rouleau & Martin, 1991). We
note that VB17 made use of different amorphous carbon optical constants but it
does not appear to greatly impact the results (Preibisch et al., 1993). We
selected an MRN distribution for the dust, using a minimum grain size of 0.1
$\micron$ and maximum grain size 0f 10000 $\micron$. This effectively assumes
that the large grains trace the smaller grains, which is not supported either
by interferometry in the mid-IR or in the sub-mm, nor follows the models of
VB17 or Menu et al. (2014), where large grains followed a separate density
pattern. That said, since our model focused solely on reproducing a 1-D SED
and 2-D scattered light model, we assume that this choice has only a minor
bearing on our results with respect to reproducing shadows. We assume that the
gas to dust ratio in the disk is 100. The composition of dust considered by
VB17 results in a surface brightness that is too faint to match the STIS
observations. This can be ameliorated by increasing the scale height of the
disk, but resulted in a poorer fit to the SED. After some trial and error we
use a composition more similar to that of Debes et al. (2013), which is based
on SED fitting by Calvet et al. (2002). For our fiducial model, we assume a
mixture of 60% Olivine and 40% water ice, but we note that this is likely a
non-unique solution.

####  4.1.5 MCFOST parameters

We input the analytical 3-D density profile as a grid into MCFOST and then
selected 107 for the number of photon packages used for the scattered light
images. Additionally, we selected a pixel size equivalent to the platescale of
a STIS pixel, we fixed the disk mass at 5$\times$10-4$M_{\odot}$, the inner
radius at 0.38 au and the outer radius 222 au to match those of VB17
accounting for the updated DR3 Gaia distance. We assumed Mie scattering for
our calculations and a single wavelength at 7000 Å that mimics the expected
central wavelength of the observations. The STIS CCD bandpass extends from
2000-10000 Å and thus the central wavelength of an observation is dependent on
the underlying SED of the source, which to zeroth order is equivalent to the
SED of TW Hya, since the dust scatters light nearly constantly across visible
wavelengths (Roberge et al., 2005). Additonally, we calculate a model PSF from
the TinyTim software 222https://www.stsci.edu/hst/instrumentation/focus-and-
pointing/focus/tiny-tim-hst-psf-modeling, assuming K7 stellar spectral type
and convolve it with the model images to compare with the STIS data.

####  4.1.6 Fiducial Model fidelity

Figure 8 shows a comparison of the fiducial model SED against the observed
photometry reported in Menu et al. (2014). The model SEDs for the fiducial
model and the subsequent tilted ring models we consider are virtually
identical, demonstrating that such a structure has a relatively minimal impact
on the one dimensional unresolved SED. Our fiducial SED fits about as well as
those of Menu et al. (2014) and VB17, though we did not attempt to optimize
the fit–the reduced $\chi^{2}$ for the fiducial fit is 3.22. Overall, the SED
is close to that observed, with an overprediction of flux between 50-200 µm,
and slight underprediction of flux in the sub-mm regime. Additionally, Figure
9 shows a comparison between the observed STIS image in 2021 and the fiducial
model. The agreement between both the SED and the general surface brightness
of the disk is excellent. Since the fiducial disk is face-on and with no
shadows, the ASB is constant.

###  4.2 Shadows from Single Tilted Rings

We next consider the resulting shadows caused by a single tilted ring at small
inclination. To tilt a restricted section of the disc, we incline the midplane
as required and ensure we are measuring the perpendicular distance from the
midplane to determine the density at a given location. If the inclination of
the ring is too small or the flaring of the disk large the shadow will have a
small depth and not extend over a large range of radii in the disk. As
inclination increases for fixed flaring, the shadow deepens and widens. At
moderate inclination the single shadow breaks into two as the star illuminates
the outer disk in between the line of nodes for the inner disk and approaches
the situation considered by Min et al. (2017). However, in these cases the
depth of the remaining shadows is much deeper than observed for TW Hya and
does not match the detailed ASB. For that reason we consider a small
inclination to the inner ring as the preferred scenario. For our fiducial
model, we find a ring at a position angle of 100$\arcdeg$ at a radius of
between 5 and 6 au with an inclination of 7$\arcdeg$ that roughly replicates
the depth, but not the detailed shape of the June 2021 ASB in Figure 7.
Further, we note that the width of the single disk shadow feature is narrower
than what was observed previously for TW Hya, implying that in past epochs the
shadow feature seen was already indicating the presence of more than one
tilted disk.

###  4.3 Shadows from Two Tilted Rings

We repeat this process but consider two independently oriented tilted rings
(See Figure 6). Here we choose two rings misaligned by 5.5° and 7° with
respect to the outer disk plane, the first with a PA=170$\arcdeg$ extending
from 5-6 au and the second from 6-7 au with PA=50$\arcdeg$. These distances
are consistent with the $\sim$16 yr timescale of the existing shadow motions
assuming the motions are Keplerian and the magnitude of the inclination is
chosen to roughly match the depth of the observed shadows, but other
configurations are likely also possible. Investigating the degneracies in the
ring location and width is beyond the scope of this paper. If we place the
rings too far out $r_{c}>15~{}au$, we cannot easily reproduce the surface
brightness beyond 30 au.

We find that the inclusion of a second ring forces slight alterations to the
fiducial inner gap depths to retain a good fit to the azimuthally averaged
surface brightness profile. In particular, the 24 au gap needs $d=0.9$ as
opposed to $d=0.65$ in the fiducial profile and the 6.68 au gap needs $d=0.92$
as opposed to $d=0.81$.

With the inclusion of the second ring, the resulting model image matches the
June 2021 image as well as the ASB (See Figures 9 and 7). We therefore
demonstrate that two mutually inclined rings at different position angles well
represent the observed disk images and represent the best explanation for the
change in behavior. At earlier times, the PAs of the two rings were closer
together and their resultant shadows overlapped, causing an apparent single
shadow.

###  4.4 Observational implications for two tilted rings

The implied radial distances of $\sim$5-7 au for the two ring model are
accessible via interferometry at visible and NIR wavelengths as well as with
direct imaging. Such observations can provide a potentially useful constraint
on the inner shadowing structure at wavelengths close to the STIS bandpass
with our images having an effective wavelength close to the Johnson-Cousins
$R$ band, but such constraints are heterogeneous in the literature. K-band
GRAVITY observations of the inner disk edge (R$\sim$0.04 AU) are consistent
with inclinations of $<20$°(Gravity Collaboration et al., 2021) a multi-
wavelength treatment of interferometry spanning NIR to radio wavelengths
suggested a disk structure with an inner radius of $\sim$0.8 au (corrected for
the updated distance to TW Hya) and a puffed up disk rim located at $\sim$3.6
au (Menu et al., 2014).

Polarized intensity imaging with SPHERE/ZIMPOL showed gaps/rings which could
be coincident with the shadowing structures we consider, as did ALMA (Andrews
et al., 2016, VB17). Interior to 28 au, the ZIMPOL observations had a
tentative detecton of an inner disk between 2-6 au and a ring at $\sim$16 au.
ALMA showed emission maxima at $\sim$3, 11, and 16 au with some additional
substructure in a plateau around 29-38 au. Most of the locations interior to
15 au are consistent with the shadows, although more stringent constraints can
be placed if the exact motion of the shadows is determined.

Additional monitoring of the shadows’ motion is needed to understand whether
the 16 yr periodicity seen is still relevant and whether it is due to orbital
motion, or mutually interacting misaligned rings. Any scenario that involves
precession will require fairly massive planetary or substellar companions and
short orbital periods to show periodicity on 16 yr timescales (Nealon et al.,
2020).

##  5 Conclusions

We report new images of the TW Hya protoplanetary disk in visible total
intensity scattered light. We have shown that the disk’s azimuthal brightness
asymmetries at r$>$40 AU, previously interpreted as arising from a misaligned
inner disk and casting a single shadow, now appears to be comprised of two
shadows, implying that a single disk at r$<$6 au is no longer a favored
explanation as the originating shadowing structure. This is the first
detection of a shadow splitting into more than one feature on a protoplanetary
disk.

We instead consider two inclined rings at roughly the same orbital separation.
This suggested explanation, when coupled with a series of gaps in dust density
and moderate disk flaring matches both the one dimensional SED as well as the
detailed features seen in scattered light. Nonetheless, the models presented
here are not necessarily unique since it is beyond the scope of this paper to
fully sample parameter space. It is true, however, that a detailed joint fit
of the SED, resolved sub-mm continuum images, interferometry data, and
scattered light in both total and polarized intensity represents a significant
number of constraints on both the disk surface structure and the presence of
interior tilted rings. The presence of shadows adds unique constraints to such
modeling, hopefully in a way that provides a more unique model for TW Hya’s
disk.

The recent behavior of the shadow throws into question the shadow rotation
period of 15.9 yr inferred by D17. This period was originally assumed to be a
precession period of the inclined disk, which required a fairly massive
companion to be consistent with observations. Additionally, the recent
discovery of a similar shadow in sub-mm CO line emission (Teague et al., 2022)
suggests that instead of pure rigid rotation, the shadow(s) may oscillate in
PA due to the mutual precession of the disks. Future observations will reveal
whether the behavior of the shadow is truly periodic, and what that period
relates to physically. If the shadow angular motions can be better
disentangled, it could open up more powerful predictions for physical
processes driving the shadow motions.

The presence of multiple tilted rings is highly suggestive of multiple
planetary companions in the innermost regions of the TW Hya disk. While
planets have been searched for in NIR data at around 6 au by VB17, no
candidates were found down to a few tens of Earth masses in the existing disk
gaps, assuming the level of disk extinction and planet evolutionary state is
well known. Additional unsuccessful planet searches have been conducted with
the Keck Vector Vortex coronagraph (Ruane et al., 2017) and in H $\alpha$
(Cugno et al., 2019; Follette et al., 2022). Intriguing compact sources in the
sub-mm have also been pointed to as evidence for planetary companions in the
outer disk (e.g. Tsukagoshi et al., 2019). Additionally, the limited radial
velocity data taken for TW Hya in the NIR is not sufficient to test the
prediction of giant planet mass companions at distances of 5-10 au by itself,
and visible radial velocities are plagued by the significant accretion and
stellar activity of TW Hya (Huélamo et al., 2008). The Gaia mission’s DR3
astrometric residual noise is not unusual compared to other stars of similar
magnitude. Limits from a combination of Hipparcos and Gaia EDR3 proper motion
anomaly measures show that the motion is consistent with no observable
acceleration and place limits of 1.32, 1.29, 1.30, and 6.69 MJup for
companions at orbits of 3, 5, 10, and 30 au, suggesting that if any planets
are causing the structures that they are likely less than these masses
(Brandt, 2021; Kervella et al., 2022). Given the inferred orbital separation
of the shadowing structures, the favorable face-on nature of the disk, and the
sensitivity that the full Gaia astrometric observations will have, it is
possible an astrometric signature of the perturbing companions could be
detected in the future. The combination of a measured planet mass and
resulting shadowing structure would mean that TW Hya will remain a key object
for the study of planet formation and disk-planet interactions.

Some of the data presented in this paper were obtained from the Mikulski
Archive for Space Telescopes (MAST) at the Space Telescope Science Institute.
The specific observations analyzed can be accessed via
https://doi.org/10.17909/mtx2-vz53 (catalog 10.17909/mtx2-vz53). We thank the
anonymous referee for their thoughtful comments. We additionally thank Emily
Rickman for helpful discussions related to companion limits to TW Hya and
Inbok Yeah for discussions related to MCFOST modeling. Support for program
#16228 was provided by NASA through a grant from the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., under NASA contract NAS5-26555. This work is also based on
observations made with the NASA/ESA Hubble Space Telescope, obtained from the
Data Archive at the Space Telescope Science Institute, which is operated by
the Association of Universities for Research in Astronomy, Inc., under NASA
contract NAS5-26555. These observations are associated with programs
#11608,13753,and 13665. The authors wish to thank the STScI Program
Coordinator, S. Meyett, for her tireless work in getting the visits for
Program 16228 scheduled successfully. RN acknowledges funding from UKRI/EPSRC
through a Stephen Hawking Fellowship (EP/T017287/1). RA gratefully
acknowledges funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No 681601), and from the Science & Technology Facilities Council
(STFC) through Consolidated Grant ST/W000857/1. This work has made use of data
from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

## References

  * Akiyama et al. (2015) Akiyama, E., et al. 2015, ApJ, 802, L17
  * Andrews et al. (2016) Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40. doi:10.3847/2041-8205/820/2/L40
  * Benisty et al. (2017) Benisty, M., Stolker, T., Pohl, A., et al. 2017, A&A, 597, A42. doi:10.1051/0004-6361/201629798
  * Benisty et al. (2018) Benisty, M., Juhász, A., Facchini, S., et al. 2018, A&A, 619, A171. doi:10.1051/0004-6361/201833913
  * Bertrang et al. (2018) Bertrang, G. H.-M., Avenhaus, H., Casassus, S., et al. 2018, MNRAS, 474, 5105. doi:10.1093/mnras/stx3052
  * Bi et al. (2021) Bi, J., Lin, M.-K., & Dong, R. 2021, ApJ, 912, 107. doi:10.3847/1538-4357/abef6b
  * Brandt (2021) Brandt, T. D. 2021, ApJS, 254, 42. doi:10.3847/1538-4365/abf93c
  * Calvet et al. (2002) Calvet, N., D’Alessio, P., Hartmann, L., et al. 2002, ApJ, 568, 1008. doi:10.1086/339061
  * Chiang & Goldreich (1997) Chiang, E. I. & Goldreich, P. 1997, ApJ, 490, 368. doi:10.1086/304869
  * Cugno et al. (2019) Cugno, G., Quanz, S. P., Hunziker, S., et al. 2019, A&A, 622, A156. doi:10.1051/0004-6361/201834170
  * Debes et al. (2013) Debes, J. H., Jang-Condell, H., Weinberger, A. J., Roberge, A., & Schneider, G. 2013, ApJ, 771, 45
  * Debes et al. (2017) Debes, J. H., Poteet, C. A., Jang-Condell, H., et al. 2017, ApJ, 835, 205. doi:10.3847/1538-4357/835/2/205
  * Debes et al. (2019) Debes, J. H., Ren, B., & Schneider, G. 2019, Journal of Astronomical Telescopes, Instruments, and Systems, 5, 035003. doi:10.1117/1.JATIS.5.3.035003
  * Dorschner et al. (1995) Dorschner, J., Begemann, B., Henning, T., et al. 1995, A&A, 300, 503
  * Follette et al. (2022) Follette, K. B., Close, L. M., Males, J. R., et al. 2022, arXiv:2211.02109. doi:10.48550/arXiv.2211.02109
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1. doi:10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1. doi:10.1051/0004-6361/202039657
  * Gravity Collaboration et al. (2021) Gravity Collaboration, Perraut, K., Labadie, L., et al. 2021, A&A, 655, A73. doi:10.1051/0004-6361/202141624
  * Herczeg & Hillenbrand (2014) Herczeg, G. J. & Hillenbrand, L. A. 2014, ApJ, 786, 97
  * Huélamo et al. (2008) Huélamo, N., Figueira, P., Bonfils, X., et al. 2008, A&A, 489, L9
  * Huang et al. (2018) Huang, J., Andrews, S. M., Cleeves, L. I., et al. 2018, ApJ, 852, 122. doi:10.3847/1538-4357/aaa1e7
  * Hughes et al. (2007) Hughes, A. M., Wilner, D. J., Calvet, N., D’Alessio, P., Claussen, M. J., & Hogerheijde, M. R. 2007, The Astrophysical Journal, 664, 536, 2007: The American Astronomical Society
  * Hughes et al. (2011) Hughes, A. M., Wilner, D. J., Andrews, S. M., et al. 2011, ApJ, 727, 85. doi:10.1088/0004-637X/727/2/85
  * Jang-Condell & Turner (2012) Jang-Condell, H., & Turner, N. J. 2012, ApJ, 749, 153
  * Ilee et al. (2022) Ilee, J. D., Walsh, C., Jennings, J., et al. 2022, arXiv:2205.01396
  * Kervella et al. (2022) Kervella, P., Arenou, F., & Thévenin, F. 2022, A&A, 657, A7. doi:10.1051/0004-6361/202142146
  * Krist et al. (2000) Krist, J. E., Stapelfeldt, K. R., Ménard, F., Padgett, D. L., & Burrows, C. J. 2000, ApJ, 538, 793
  * Marino et al. (2015) Marino, S., Perez, S., & Casassus, S. 2015, ApJ, 798, L44. doi:10.1088/2041-8205/798/2/L44
  * Mekkaden (1998) Mekkaden, M. V. 1998, A&A, 340, 135
  * Nealon et al. (2019) Nealon, R., Pinte, C., Alexander, R., et al. 2019, MNRAS, 484, 4951. doi:10.1093/mnras/stz346
  * Nealon et al. (2020) Nealon, R., Price, D. J., & Pinte, C. 2020, MNRAS, 493, L143. doi:10.1093/mnrasl/slaa026
  * Macías et al. (2021) Macías, E., Guerra-Alvarado, O., Carrasco-González, C., et al. 2021, A&A, 648, A33. doi:10.1051/0004-6361/202039812
  * Menu et al. (2014) Menu, J., van Boekel, R., Henning, T., et al. 2014, A&A, 564, A93. doi:10.1051/0004-6361/201322961
  * Min et al. (2017) Min, M., Stolker, T., Dominik, C., et al. 2017, A&A, 604, L10. doi:10.1051/0004-6361/201730949
  * Muro-Arena et al. (2020) Muro-Arena, G. A., Benisty, M., Ginski, C., et al. 2020, A&A, 635, A121. doi:10.1051/0004-6361/201936509
  * Pinilla et al. (2018) Pinilla, P., Benisty, M., de Boer, J., et al. 2018, ApJ, 868, 85. doi:10.3847/1538-4357/aae824
  * Pinte et al. (2006) Pinte, C., Ménard, F., Duchêne, G., et al. 2006, A&A, 459, 797. doi:10.1051/0004-6361:20053275
  * Pinte et al. (2009) Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967. doi:10.1051/0004-6361/200811555
  * Poteet et al. (2018) Poteet, C. A., Chen, C. H., Hines, D. C., et al. 2018, ApJ, 860, 115. doi:10.3847/1538-4357/aac2e4
  * Preibisch et al. (1993) Preibisch, T., Ossenkopf, V., Yorke, H. W., et al. 1993, A&A, 279, 577
  * Qi et al. (2004) Qi, C., Ho, P. T. P., Wilner, D. J., Takakuwa, S., Hirano, N., Ohashi, N., Bourke, T. L., Zhang, Q., Blake, G. A., Hogerheijde, M., Saito, M., Choi, M., & Yang, J. 2004, ApJ, 616, L11
  * Qi et al. (2013) Qi, C., Öberg, K. I., Wilner, D. J., D’Alessio, P., Bergin, E., Andrews, S. M., Blake, G. A., Hogerheijde, M. R., & van Dishoeck, E. F. 2013, Science, 341, 630
  * Roberge et al. (2005) Roberge, A., Weinberger, A. J., & Malumuth, E. M. 2005, ApJ, 622, 1171
  * Rosenfeld et al. (2012) Rosenfeld, K. A., Qi, C., Andrews, S. M., et al. 2012, ApJ, 757, 129
  * Rouleau & Martin (1991) Rouleau, F. & Martin, P. G. 1991, ApJ, 377, 526. doi:10.1086/170382
  * Ruane et al. (2017) Ruane, G., Mawet, D., Kastner, J., et al. 2017, AJ, 154, 73. doi:10.3847/1538-3881/aa7b81
  * Schneider et al. (2017) Schneider, G., Gaspar, A., Debes, J., et al. 2017, Instrument Science Report STIS 2017-3, 35 pages
  * Sokal et al. (2018) Sokal, K. R., Deen, C. P., Mace, G. N., et al. 2018, ApJ, 853, 120. doi:10.3847/1538-4357/aaa1e4
  * Teague et al. (2018) Teague, R., Henning, T., Guilloteau, S., et al. 2018, ApJ, 864, 133. doi:10.3847/1538-4357/aad80e
  * Teague et al. (2019) Teague, R., Bae, J., Huang, J., et al. 2019, ApJ, 884, L56. doi:10.3847/2041-8213/ab4a83
  * Teague et al. (2022) Teague, R., Bae, J., Benisty, M., et al. 2022, arXiv:2204.06691
  * Tsukagoshi et al. (2019) Tsukagoshi, T., Muto, T., Nomura, H., et al. 2019, ApJ, 878, L8. doi:10.3847/2041-8213/ab224c
  * Ueda et al. (2020) Ueda, T., Kataoka, A., & Tsukagoshi, T. 2020, ApJ, 893, 125. doi:10.3847/1538-4357/ab8223
  * Uyama et al. (2020) Uyama, T., Currie, T., Christiaens, V., et al. 2020, ApJ, 900, 135. doi:10.3847/1538-4357/aba8f6
  * Vacca & Sandell (2011) Vacca, W. D. & Sandell, G. 2011, ApJ, 732, 8
  * van Boekel et al. (2017) van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132. doi:10.3847/1538-4357/aa5d68
  * Weinberger et al. (2002) Weinberger, A. J., Becklin, E. E., Schneider, G., Chiang, E. I., Lowrance, P. J., Silverstone, M. D., Zuckerman, B., Hines, D. C., & Smith, B. A. 2002, The Astrophysical Journal, 566, 409, 2002: The American Astronomical Society
  * Wolff et al. (2016) Wolff, S. G., Perrin, M., Millar-Blanchaer, M. A., et al. 2016, ApJ, 818, L15. doi:10.3847/2041-8205/818/1/L15
  * Woodgate et al. (1998) Woodgate, B. E., Kimble, R. A., Bowers, C. W., et al. 1998, PASP, 110, 1183. doi:10.1086/316243


# Efficient adapter transfer of Self-Supervised speech models for Automatic
Speech Recognition

###### Abstract

Self-supervised learning (SSL) is a powerful tool that allows learning of
underlying representations from unlabeled data. Transformer based models such
as wav2vec 2.0 and HuBERT are leading the field in the speech domain.
Generally these models are fine-tuned on a small amount of labeled data for a
downstream task such as Automatic Speech Recognition (ASR). This involves re-
training the majority of the model for each task. Adapters are small
lightweight modules which are commonly used in Natural Language Processing
(NLP) to adapt pre-trained models to new tasks. In this paper we propose
applying adapters to wav2vec 2.0 to reduce the number of parameters required
for downstream ASR tasks, and increase scalability of the model to multiple
tasks or languages. Using adapters we can perform ASR while training fewer
than $10$% of parameters per task compared to full fine-tuning with little
degradation of performance. Ablations show that applying adapters into just
the top few layers of the pre-trained network gives similar performance to
full transfer, supporting the theory that higher pre-trained layers encode
more phonemic information, and further optimizing efficiency.

Index Terms—  Automatic Speech Recognition, Self-Supervision, Adapters,
Transfer Learning

##  1 Introduction

Self-supervision is now standard practice in Natural Language Processing
(NLP), with models such as BERT [1] achieving state of the art results. More
recently, self-supervised learning (SSL) approaches have been applied to
speech tasks with great success. Recent works such as wav2vec 2.0 [2] and
HuBERT [3] show that self-supervision in speech can achieve state of the art
results.

Most SSL models rely on unsupervised pre-training followed by supervised fine-
tuning on a downstream task. The unsupervised pre-training allows a model to
benefit from the large amounts of unlabeled data which are readily available.
Fine-tuning uses a comparatively small amount of labeled data and retrains the
majority of the self-supervised model to the desired downstream task.

Fig. 1: A comparison of trained model parameters vs WER. Pre-trained self-
supervised models are either fine-tuned or adapted for ASR using 10 hours of
supervised data and evaluated on LibriSpeech dev-clean. The Y axis shows
various iterations of fine-tuning/adapters with different numbers of layers
trained, and thus different numbers of trainable parameters. Adapters achieve
a similar performance in WER to fine-tuning with only a fraction of
parameters. Training just layer normalizations performs poorly.

While these approaches gain excellent results, fine-tuning the model is
computationally expensive and does not scale well to multiple tasks, such as
in the case of multi-lingual Automatic Speech Recognition (ASR); a complete
set of fine-tuned parameters must be learnt and stored per downstream task.
This is commonly in the scale of $O(10^{8})$ parameters per task. Once the
model is fine-tuned for one task, the entire model is fixed for that task, and
the base model must be reloaded to transfer to future tasks.

Adapters are small trainable modules that can be applied into the layers of a
frozen pre-trained network and tuned to a particular task. Recently these have
been applied to pre-trained unsupervised models in NLP [4, 5]. The benefit of
adapters is that they allow adaptation to new tasks with a relatively small
number of parameters per task. As the existing model parameters remain frozen,
the original model can support multiple downstream tasks, only a set of
adapter parameters are required per task. This makes adapters more efficient
to train, and highly scalable to multiple tasks.

Fig. 2: Left: structure of wav2vec 2.0 model with task specific adapter
modules. f($\cdot$) is a convolutional encoder which is followed by 12
transformer encoder blocks. For downstream ASR a linear classifier,
CTC($\cdot$), is applied to the output of the transformer blocks. For adapter
transfer, adapter modules are inserted into the model, and during adaptation
to the downstream ASR task only the adapters, layer normalization layers and
linear classifier are trained, the rest of the network is frozen. Right:
structure of adapter module consisting of a down-projection, a non-linearity
and an up-projection, with a skip connection.

The task of ASR is very well researched for high resource languages such as
English. However, it is more difficult to gain good performance on languages
where there is limited paired text and audio training data. wav2vec 2.0 shows
good performance with limited labeled data, and multi-lingual wav2vec 2.0 [6]
shows good results across a range of languages including languages unseen
during pre-training. Self-supervised speech models are able to leverage
generic features from pre-training to adapt well to scenarios with limited
labeled data and are thus ideal for multi-lingual ASR. However, it is
impractical to train and store fine-tuned parameters for every language.
Furthermore, fine-tuning overwrites pre-trained model parameters so may not be
the optimal way of utilizing pre-training knowledge.

In this paper we apply adapters to a pre-trained wav2vec 2.0 model in order to
transfer the pre-trained representations to the task of downstream ASR.
Adapters suggest a way of utilizing these pre-trained representations in a
much more efficient manner, with many fewer parameters required per downstream
task than usual fine-tuning methods. This has applications for multi-task
scenarios such as multi-lingual ASR, and also increases the accessibility of
self-supervised speech models for research. To the best of our knowledge this
is the first time adapter modules have been applied to a self-supervised
speech model for ASR.

##  2 Background

###  2.1 Self-supervised ASR

Labeled speech data is relatively scarce and expensive to generate. Self-
supervised methods are able to take advantage of the huge amount of unlabeled
speech data available to learn generic speech features.

There have been various proposed methods such as CPC [7], APC [8], and wav2vec
[9]. All aim to learn an encoded representation of a raw speech waveform using
various self-supervised tasks which train the underlying network. Many of
these models include some form of reconstruction, where masked
representations, or elements of the original input waveform, are recreated.
The underlying network can then be used to extract representations which are
used as input to downstream tasks. This is instead of more conventional speech
features such as log Mel filter-bank features. These networks can either be
frozen, or further optimized to the downstream task by some form of transfer
learning.

wav2vec 2.0 is a leading model in this field [2]. However, the usual method of
fine-tuning wav2vec 2.0 for downstream ASR requires re-training most of the
model layers. In fact $95.6\%$ of the total model parameters are trained, and
must therefore be stored, per task. This is not an issue if monolingual ASR is
the only target task. However it is not ideal for multi-lingual speech
recognition, or indeed utilizing the wav2vec 2.0 encoder for multiple tasks
such as speech translation, as in [10].

Recently advances have been made to the wav2vec 2.0 architecture including
adding self-training [11]. While we use wav2vec 2.0 in this work, our approach
would scale well to any self-supervised transformer based speech model,
including HuBERT [3].

###  2.2 Adapters

Adapters were introduced in NLP as an efficient method of transfer learning
[4]. It was found that adapters approach the performance of full fine-tuning
with only a fraction of the parameters in NLP tasks using a pre-trained BERT
model. Adapters have also been successfully applied to NMT [12] and vision
tasks [13]. Recently dual adapters have been proposed for multi-lingual
transfer [5], where one adapter module is used to adapt for language, and
another captures task specific adaptations.

Adapters have also been applied to speech in streaming RNNT models [14] and
hybrid CTC-attention speech transformers [15] to solve the issue of multi-
lingual speech recognition. More recently adapters have been used in speech
translation to combine pre-trained modules [10]. Adapters have also been
employed for efficient SSL pre-training of new tasks in a continual learning
setting [16].

We hypothesize that we will see the same benefits of adapters in a speech
model as in an NLP model, namely parameter efficient transfer of the pre-
trained network to a downstream task with little performance degradation.

##  3 Method

In wav2vec 2.0 a convolutional encoder takes as input the raw waveform. Then a
transformer model is applied to the encoded output. The model is trained by
masking input frames and comparing predicted values with positive and negative
quantized versions of the masked frames. This is similar to contrastive
predictive coding [7].

When applying wav2vec 2.0 to a downstream model the feature encoder is frozen,
and a linear classifier is added on top of the transformer context model for
fine-tuning. Generally the transformer model is also frozen for the first N
updates. When applying wav2vec 2.0 to ASR, the linear classifier and model are
fine-tuned with a Connectionist Temporal Classification (CTC) loss using a
small amount of labeled speech data.

Adapters are small bottleneck modules consisting of a down projection, a non-
linearity, and an up-projection, with a skip connection (see Fig. 2). The
initial implementation [4] applies adapters after both the self-attention and
feedforward layers. However it is possible to apply adapters in different
positions throughout the transformer block [12]. The fully connected layers
are initialized as a near identity function. The identity initialization and
the skip connection allow the module to be ignored if not deemed necessary
during training.

In our experiments with adapters we apply adapters twice in each transformer
block. Our experimentation showed that this gave the best results. We apply a
linear classifier on top of the transformer network. A set of adapters, layer
normalization layers and a linear classifier are trained per task using a CTC
loss, the rest of the network remains frozen.

##  4 Experiments and Results

We use the standard wav2vec 2.0 BASE architecture which contains 12
transformer layers, and use the publicly released pre-trained checkpoint from
fairseq [17] which has been pre-trained using the LibriSpeech [18] dataset. We
use a standard size of 256 for all adapters, and initially apply adapters into
every transformer layer.

We also run our own fine-tuning experiments as a comparison and follow the
fine-tuning format of wav2vec 2.0. By tuning hyper-parameters we are able to
improve on the word error rate (WER) values reported in that paper [2]. We
trained for 20k steps, with the transformer layers frozen for just the first
4k updates, and a learning rate of 5e-5.

Table 1: A pre-trained BASE wav2vec2 model is transferred to the downstream
ASR task using the 10 hour LibriLight (LL) supervised set, or a random 10 hour
subset of the French Common Voice (CV) corpus. Word error rate (WER) is
reported on the dev-clean/dev-other sets of LibriSpeech for English, and the
CV test set for French. Results are all without a language model. % of
trainable parameters compared to the total model parameters are also reported

| wav2vec 2.0 FT[2] | Fine-tune | Adapter
---|---|---|---
10h LL dev-clean | $10.9\%$ | $8.98\%$ | $9.39\%$
10h LL dev-other | $17.4\%$ | $16.9\%$ | $17.0\%$
10h French CV test | N/A | $40.2\%$ | $39.4\%$
% trained params | $95.6\%$ | $95.6\%$ | $9.2\%$

We investigate performance of fine-tuning and adapters on the 10 hour
supervised subset of the LibriLight (LL) dataset [19], and evaluate on
standard LibriSpeech dev sets. We also experiment with French ASR to
demonstrate the multi-task scenario. We take a 10h subset of the Common Voice
(CV) [20] corpus and evaluate on the CV test set. We calculate WER in all
cases. We found that the optimal setup for adapter transfer was to run for 10k
steps with a learning rate of 5e-4. All experiments are run on 8 V100 GPUs,
and without language model fusion to enable distraction free comparison of
transfer method.

Fig. 3: A summary of number of trained parameters per task with fine-tuning
and adapters. Each fine-tuning task nearly doubles the number of trained
parameters which must be learnt and stored. Whereas adapters add only a small
number of additional parameters per task showing their scalability to multi-
task scenarios.

Adapters perform slightly worse than fine-tuning on English ASR (see table 1),
however the absolute WER increase is just $0.41\%$ and $0.17\%$ for dev-clean
and dev-other respectively. This performance decrease comes with a huge
decrease in the number of trained task-specific parameters. Fine-tuning
requires training $95.6\%$ of parameters, our adapter approach only trains
$9.2\%$ of parameters. For French ASR, there is in fact a slight performance
increase when using adapters compared to traditional fine-tuning.

However, the real benefit of this approach comes when scaling to multiple
languages or tasks, as can be seen in Fig. 3. Each fine-tuning task nearly
doubles the required number of parameters which must be learnt and stored,
however adapters add only a small number of additional parameters per task.
This makes adapter transfer much more scalable than fine-tuning while
attaining a similar performance.

Our results show that even when there is a mismatch between pre-trained and
downstream language, in the case of French transfer, both fine-tuning and
adapters are able to achieve some success in ASR. This is unsurprising for
fine-tuning, since the language specific features of the pre-trained network
are overwritten during training. However, with adapter transfer, the original
network remains unchanged. Therefore the adapters are able to compensate for
the language mismatch, and learn both language and task specific features.
This is in contrast to [5] which uses separate adapters for task and language.

It is also worth noting that adapter experiments run more quickly than fine-
tuning experiments, both due to the reduction in trainable parameters which
increases training speed, and the smaller optimal number of training steps.
Adapter experiments do not depend on a freeze steps hyper-parameter, and run
successfully with a simple bi-stage learning rate schedule, as in [4], rather
than the more complex tri-stage scheduler required for traditional fine-tuning
[2]. All these factors make this adapter approach much more experimentally
friendly for researchers.

Table 2: A pre-trained bilingual (French and English) model is transferred to English and French downstream ASR using fine-tuning and adapters. Results are reported in WER. | Fine-tune | Adapter
---|---|---
10h LL dev-clean | $12.16\%$ | $12.91\%$
10h French CV test | $27.75\%$ | $28.25\%$

We also pre-trained our own bi-lingual English and French wav2vec 2.0 model
using approximately 1000 hours of French CV data, as well as 960h of English
LibriSpeech data, and ran English and French ASR experiments (see table 2).
This demonstrates that our adapter method is also valid for multi-lingual pre-
trained models as adapters again get within close performance of fine-tuning.

Finally, we tested training just the layer normalization parameters as
suggested in [4]. However, as found in that work, this approach does not
perform well, WER on LibriSpeech dev-clean was $32.9\%$. This provides
evidence that the adapter modules themselves improve performance, rather than
the re-trained layer normalization parameters.

###  4.1 Ablations

Prior work suggests that lower layers of self-supervised models contain more
generic speech features, and higher layers contribute more to phone
discrimination [8]. Therefore we investigate training just the top N layers of
the network using both fine-tuning and adapters. All experiments are run with
the BASE English wav2vec 2.0 model, 10 hours LibriLight training data and
evaluated on LibriSpeech dev-clean.

Using adapters in just the top 4 layers, out of a total 12 layers, gives
nearly as good performance as adding adapters into every layer, and in fact
using just 6 adapters gives the best performance with $9.27\%$ WER. This
immediately halves the number of required parameters to just $4.85$% of total
parameters trained. It is also possible to optimize fine-tuning; training just
the top 8 transformer layers gives best performance. However fine-tuning
requires more layers to be trained than using adapters, and even 8 layers
equates to training $65.5\%$ of parameters (see Fig. 4). When just one layer
is trained, the method with adapters performs significantly better showing
that adapters are better able to utilize the pre-training knowledge of the
entire network.

The curves presented in Fig. 4 show that the top layers of the network are
more important for downstream performance than the bottom layers, supporting
the hypothesis that the higher layers of the network encode more phonemic
information, and lower layers more generic speech information.

Fig. 4: A comparison of trainable layers vs WER. For fine-tuning experiments,
transformer layers of the model are optionally trained or frozen in a top down
manner and evaluated for WER on LibriSpeech dev-clean. For adapter
experiments, adapters are only inserted into the top N layers. Models with
adapters perform better than models which are fine-tuned when the number of
trained layers is limited.

##  5 Discussion & Conclusion

This work provides an insight into how self-supervised speech models can be
utilized in a more parameter efficient manner without sacrificing performance.
When transferring to a downstream task, fine-tuning the majority of the model
is still generally the best performing method. However, using adapters instead
of fine-tuning achieves close to that performance with a fraction of
parameters and takes less training time. This allows quicker and cheaper
experimentation with models such as wav2vec 2.0 and HuBERT, and increases
scalability to multiple tasks, for example multi-lingual ASR. More work could
be done on utilizing adapters for additional tasks such as speaker recognition
and speech translation.

Ablations show that it is unnecessary to transfer all layers of the network to
the downstream task, only the top N layers, thereby supporting the theory that
the higher layers of these pre-trained networks encode more phonemic
information. Adapters are better able to utilize pre-trained information as we
achieve better performance with one layer of adapters, than with one layer of
fine-tuning, and best performance comes from adapting fewer layers compared to
fine-tuning.

While we are the first to utilize adapters in this way for speech, our
findings are similar to those in the NLP domain [4, 5]. It would be
interesting to perform the same layer ablations in the NLP domain. More
broadly, we show that the speech domain can benefit from future NLP work on
SSL.

## References

  * [1] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova,  “Bert: Pre-training of deep bidirectional transformers for language understanding,”  arXiv preprint arXiv:1810.04805, 2018.
  * [2] Alexei Baevski, Henry Zhou, Abdelrahman Mohamed, and Michael Auli,  “wav2vec 2.0: A framework for self-supervised learning of speech representations,”  arXiv preprint arXiv:2006.11477, 2020.
  * [3] Wei-Ning Hsu, Benjamin Bolte, Yao-Hung Hubert Tsai, Kushal Lakhotia, Ruslan Salakhutdinov, and Abdelrahman Mohamed,  “Hubert: Self-supervised speech representation learning by masked prediction of hidden units,”  arXiv preprint arXiv:2106.07447, 2021.
  * [4] Neil Houlsby, Andrei Giurgiu, Stanislaw Jastrzebski, Bruna Morrone, Quentin De Laroussilhe, Andrea Gesmundo, Mona Attariyan, and Sylvain Gelly,  “Parameter-efficient transfer learning for nlp,”  in International Conference on Machine Learning. PMLR, 2019, pp. 2790–2799.
  * [5] Jonas Pfeiffer, Ivan Vulić, Iryna Gurevych, and Sebastian Ruder,  “Mad-x: An adapter-based framework for multi-task cross-lingual transfer,”  arXiv preprint arXiv:2005.00052, 2020.
  * [6] Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov,  “Unsupervised cross-lingual representation learning at scale,”  arXiv preprint arXiv:1911.02116, 2019.
  * [7] Aaron van den Oord, Yazhe Li, and Oriol Vinyals,  “Representation learning with contrastive predictive coding,”  arXiv preprint arXiv:1807.03748, 2018.
  * [8] Yu-An Chung, Wei-Ning Hsu, Hao Tang, and James Glass,  “An unsupervised autoregressive model for speech representation learning,”  arXiv preprint arXiv:1904.03240, 2019.
  * [9] Steffen Schneider, Alexei Baevski, Ronan Collobert, and Michael Auli,  “wav2vec: Unsupervised pre-training for speech recognition,”  arXiv preprint arXiv:1904.05862, 2019.
  * [10] Hang Le, Juan Pino, Changhan Wang, Jiatao Gu, Didier Schwab, and Laurent Besacier,  “Lightweight adapter tuning for multilingual speech translation,”  arXiv preprint arXiv:2106.01463, 2021.
  * [11] Qiantong Xu, Alexei Baevski, Tatiana Likhomanenko, Paden Tomasello, Alexis Conneau, Ronan Collobert, Gabriel Synnaeve, and Michael Auli,  “Self-training and pre-training are complementary for speech recognition,”  in ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2021, pp. 3030–3034.
  * [12] Ankur Bapna, Naveen Arivazhagan, and Orhan Firat,  “Simple, scalable adaptation for neural machine translation,”  arXiv preprint arXiv:1909.08478, 2019.
  * [13] Sylvestre-Alvise Rebuffi, Hakan Bilen, and Andrea Vedaldi,  “Learning multiple visual domains with residual adapters,”  arXiv preprint arXiv:1705.08045, 2017.
  * [14] Anjuli Kannan, Arindrima Datta, Tara N Sainath, Eugene Weinstein, Bhuvana Ramabhadran, Yonghui Wu, Ankur Bapna, Zhifeng Chen, and Seungji Lee,  “Large-scale multilingual speech recognition with a streaming end-to-end model,”  arXiv preprint arXiv:1909.05330, 2019.
  * [15] Genta Indra Winata, Guangsen Wang, Caiming Xiong, and Steven Hoi,  “Adapt-and-adjust: Overcoming the long-tail problem of multilingual speech recognition,”  arXiv preprint arXiv:2012.01687, 2020.
  * [16] Samuel Kessler, Bethan Thomas, and Salah Karout,  “Continual-wav2vec2: an application of continual learning for self-supervised automatic speech recognition,”  arXiv preprint arXiv:2107.13530, 2021.
  * [17] Myle Ott, Sergey Edunov, Alexei Baevski, Angela Fan, Sam Gross, Nathan Ng, David Grangier, and Michael Auli,  “fairseq: A fast, extensible toolkit for sequence modeling,”  in Proceedings of NAACL-HLT 2019: Demonstrations, 2019.
  * [18] Vassil Panayotov, Guoguo Chen, Daniel Povey, and Sanjeev Khudanpur,  “Librispeech: an asr corpus based on public domain audio books,”  in 2015 IEEE international conference on acoustics, speech and signal processing (ICASSP). IEEE, 2015, pp. 5206–5210.
  * [19] Jacob Kahn, Morgane Rivière, Weiyi Zheng, Evgeny Kharitonov, Qiantong Xu, Pierre-Emmanuel Mazaré, Julien Karadayi, Vitaliy Liptchinsky, Ronan Collobert, Christian Fuegen, et al.,  “Libri-light: A benchmark for asr with limited or no supervision,”  in ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2020, pp. 7669–7673.
  * [20] Rosana Ardila, Megan Branson, Kelly Davis, Michael Henretty, Michael Kohler, Josh Meyer, Reuben Morais, Lindsay Saunders, Francis M Tyers, and Gregor Weber,  “Common voice: A massively-multilingual speech corpus,”  arXiv preprint arXiv:1912.06670, 2019.