# Calibrating VLBI Polarization Data Using GPCAL. II. Time-Dependent
Calibration

Jongho Park Department of Astronomy and Space Science, Kyung Hee University,
1732, Deogyeong-daero, Giheung-gu, Yongin-si, Gyeonggi-do 17104, Republic of
Korea  Korea Astronomy and Space Science Institute, Daedeok-daero 776,
Yuseong-gu, Daejeon 34055, Republic of Korea  Institute of Astronomy and
Astrophysics, Academia Sinica, P.O. Box 23-141, Taipei 10617, Taiwan  Keiichi
Asada Institute of Astronomy and Astrophysics, Academia Sinica, P.O. Box
23-141, Taipei 10617, Taiwan  Do-Young Byun Korea Astronomy and Space Science
Institute, Daedeok-daero 776, Yuseong-gu, Daejeon 34055, Republic of Korea
University of Science and Technology, Gajeong-ro 217, Yuseong-gu, Daejeon
34113, Republic of Korea

(Received –; Revised –; Accepted –)

###### Abstract

We present a new method of time-dependent instrumental polarization
calibration for Very Long Baseline Interferometry (VLBI). This method has been
implemented in the recently developed polarization calibration pipeline GPCAL.
Instrumental polarization, also known as polarimetric leakage, is a direction-
dependent effect, and it is not constant across the beam of a telescope.
Antenna pointing model accuracy is usually dependent on time, resulting in
off-axis polarimetric leakages that can vary with time. The method is designed
to correct for the off-axis leakages with large amplitudes that can severely
degrade linear polarization images. Using synthetic data generated based on
real Very Long Baseline Array (VLBA) data observed at 43 GHz, we evaluate the
performance of the method. The method was able to reproduce the off-axis
leakages assumed in the synthetic data, particularly those with large
amplitudes. The method has been applied to two sets of real VLBA data and the
derived off-axis leakages show very similar trends over time for pairs of
nearby sources. Furthermore, the amplitudes of the off-axis leakages are
strongly correlated with the antenna gain correction factors. The results
demonstrate that the method is capable of correcting for the off-axis leakages
present in VLBI data. By calibrating time-dependent instrumental polarization,
the rms-noise levels of the updated linear polarization images have been
significantly reduced. The method is expected to substantially enhance the
quality of linear polarization images obtained from existing and future VLBI
observations.

high angular resolution — techniques: interferometric — techniques:
polarimetric — methods: data analysis

††journal: ApJ††facilities: VLBA (NRAO)††software: AIPS (Greisen, 2003),
Difmap (Shepherd, 1997), GPCAL (Park et al., 2021b), ParselTongue (Kettenis et
al., 2006), Scipy (Virtanen et al., 2020)

##  1 Introduction

Polarization observations using Very Long Baseline Interferometry (VLBI) are
ideally suited for studying the processes of mass accretion and jet formation
(see, e.g., Park & Algaba 2022 and references therein). These processes occur
at small physical scales and can only be observed with instruments with high
angular resolution. Polarized emission from plasma in mass accretion flows and
jets is directly related to their magnetic fields. Recent observations of the
nearby elliptical galaxy M87 with the Event Horizon Telescope (EHT)
demonstrate the power of VLBI polarimetry (Event Horizon Telescope
Collaboration et al., 2019a, b, c, d, e, f). The EHT total intensity image of
M87 reveals the presence of a prominent ring-like structure. This structure is
interpreted as the result of gravitational light bending of synchrotron
emission from a hot plasma surrounding the black hole, along with photon
capture occurring at the event horizon. In the corresponding linear
polarization images of the ring, the polarization position angles are arranged
in a nearly azimuthal pattern (Event Horizon Telescope Collaboration et al.,
2021a). The comparison between these images and the General Relativistic
Magnetohydrodynamic (GRMHD) simulation results indicates that a strong
ordered, poloidal-dominated magnetic field may exist around the black hole,
which is capable of generating the powerful jets seen in this source. (Event
Horizon Telescope Collaboration et al., 2021b). In addition, polarization
observations of AGN jets using VLBI have provided important information about
the collimation, acceleration, and particle acceleration processes in AGN jets
(e.g., Asada et al., 2002; Jorstad et al., 2017; Gabuzda, 2018; Lister et al.,
2018; Park et al., 2019a, 2021a; Lisakov et al., 2021).

The linear polarization of AGN jets typically ranges from a few to a few tens
of percent (e.g., Jorstad et al., 2017; Lister et al., 2018). In the past, the
quality of VLBI linear polarization images of AGN jets was determined
primarily by their sensitivity. However, recent VLBI arrays have increased
their sensitivity considerably by enlarging the recording bandwidth and
including very high sensitivity stations such as the phased Atacama Large
Millimeter/submillimeter Array (ALMA; Matthews et al. 2018) and the phased
Very Large Array (VLA). Consequently, systematic errors are becoming a more
critical factor in determining the quality of VLBI linear polarization images.
The most dominant systematic error affecting VLBI polarization data is
“instrumental polarization”. The source of this polarization signal is due to
the antenna polarization not being exactly circular or linear, causing an
unpolarized source to appear polarized. This polarization is commonly referred
to as “polarization leakage” or “D-terms” and it must be corrected in the
visibility data before producing an image.

LPCAL is a task incorporated into Astronomical Image Processing System (AIPS,
Greisen, 2003) and based on the linearized leakage model (Leppanen et al.,
1995). It has been a standard program for calibrating instrumental
polarization in VLBI data for a long time. In spite of its success for many
studies utilizing various VLBI arrays (e.g., Casadio et al., 2017; Jorstad et
al., 2017; Lister et al., 2018; Park et al., 2018, 2019a; Gómez et al., 2022;
Zhao et al., 2022), there are some limitations which prevent accurate
calibration. An important limitation is that the “similarity
approximation”111As far as we are aware, the term “similarity approximation”
was initially introduced in Leppanen et al. (1995). The concept was originally
proposed by Cotton (1993), who also noted that the assumption of a direct
proportionality between total intensity and linear polarization emission may
not be valid., which assumes similar total intensity and linear polarization
structures for the calibrators (Cotton, 1993; Leppanen et al., 1995), is often
violated for VLBI. The problem becomes more severe at high frequencies, where
nearly all calibrators are resolved. Recently, a number of new calibration and
imaging pipelines have been developed to improve calibration accuracy and
linear polarization images, including the Generalized Polarization Calibration
pipeline (GPCAL; Park et al. 2021b), polsolve (Martí-Vidal et al., 2021), the
eht-imaging software library (Chael et al., 2016, 2018), D-term Modeling Code
(DMC, Pesce 2021), THEMIS (Broderick et al., 2020), and Comrade (Tiede, 2022).
Some of the pipelines were applied to the first linear polarization imaging of
the M87 black hole (Event Horizon Telescope Collaboration et al., 2021a).

Nevertheless, most of the existing pipelines rely on the fundamental
assumptions that instrumental polarization remains constant over a wide range
of frequencies of a receiver and over a period of time during observation.
VLBI arrays of recent generations can violate these assumptions, which results
in systematic errors in the visibility data. This series of papers presents
new methods for modeling the frequency- and time-dependent instrumental
polarization in VLBI data, which have been implemented in GPCAL. The method
for correcting frequency-dependent instrumental polarization is presented in a
companion paper (Park et al. 2023; henceforth Paper I). The purpose of this
paper is to introduce the method for correcting for time-dependent
instrumental polarization.

Instrumental polarization itself is believed to not change significantly
during the observation. However, it is a direction-dependent effect (e.g.,
Smirnov, 2011), so the instrumental polarization varies depending on the
direction of the antenna beam (due to the cross polarized sidelobes; see,
e.g., Napier 1989; Thum et al. 2008). More specifically, instrumental
polarization of an antenna consists of two main components (see Section 4.3 in
Napier 1989). There is one at the center of the antenna beam that can be
considered constant across the beam ($D_{\rm on-axis}$). The other component
varies across the beam ($D_{\rm off-axis}$). The accuracy of antenna pointing
during observation can be subject to variations over time, arising from
stochastic winds and the deformation of antennas caused by sunlight, among
other factors. Thus, effective D-terms can vary over time during an
observation due to the direction-dependent effect. Only if the antenna
pointing is sufficiently accurate throughout the observing period does
instrumental polarization remain constant. Time-dependent leakage amplitudes
are expected to be directly related to antenna pointing accuracy in this case.
In general, high-frequency observations and antennas with large diameters
usually have less accurate antenna pointing due to small beam sizes and
stronger dish deformation at low elevations. Moreover, the instrumental
polarization of phased arrays (such as the phased ALMA and VLA) is variable
(e.g., Event Horizon Telescope Collaboration et al., 2019c), as is the phasing
efficiency. As demonstrated in this paper, even the D-terms for the Very Long
Baseline Array (VLBA) may vary greatly due to inaccurate antenna pointing and
changing weather conditions. GPCAL can increase the dynamic range of the VLBA
linear polarization images by a factor of several by correcting the time-
dependent instrumental polarization.

The paper is organized as follows. In Section 2, we describe the radio
interferometer measurement equation, which forms the basis of the GPCAL
instrumental polarization model. We discuss in Section 3 the strategy for
calibrating time-dependent polarization leakages. We validate the performance
of the method with a synthetic data set in Section 4. In Section 5, we apply
the method to real VLBA data sets and verify that time-dependent leakage
correction using GPCAL can significantly enhance the quality of VLBI linear
polarization images. We summarize and conclude in Section 6.

##  2 The Radio Interferometer Measurement Equation

We use the radio interferometer measurement equation (RIME; Hamaker et al.
1996; Sault et al. 1996; Hamaker & Bregman 1996; Hamaker 2000; Smirnov 2011),
as described in Section 2 in Paper I. In this paper, we provide a brief
explanation of the equation for the convenience of readers.

The observed complex visibilities between two VLBI antennas, $m$ and $n$, are
expressed in a visibility matrix, ${\bm{V}}_{mn}$, of the form

${\bm{V}}_{mn}=\begin{pmatrix}r^{RR}_{mn}&r^{RL}_{mn}\\\
r^{LR}_{mn}&r^{LL}_{mn}\end{pmatrix},$ (1)

where $R$ and $L$ refer to the right- and left-handed circular polarizations
(RCP and LCP), respectively. The observed $V_{mn}$ is corrupted by antenna
gain and polarization leakage. It is convenient to arrange all the corruptions
into a single Jones matrix (Jones, 1941):

$\displaystyle{\bm{J}}_{m}$ $\displaystyle=$
$\displaystyle{\bm{G}}_{m}{\bm{D}}_{m}{\bm{P}}_{m}$ (2) $\displaystyle=$
$\displaystyle\begin{pmatrix}G^{R}_{m}&0\\\
0&G^{L}_{m}\end{pmatrix}\begin{pmatrix}1&D^{R}_{m}\\\
D^{L}_{m}&1\end{pmatrix}\begin{pmatrix}e^{j\phi_{m}}&0\\\
0&e^{-j\phi_{m}}\end{pmatrix},$

where $G$ is the complex antenna gain, $D$ is the leakage factor (D-term), and
$\phi$ is the antenna field rotation angle. Subscripts and superscripts denote
antenna numbers and polarization, respectively. The field rotation angle is a
combination of the source’s elevation ($\theta_{\rm el}$), parallactic angle
($\psi_{\rm par}$), and a constant offset for the rotation of antenna feed
with respect to the azimuth axis ($\phi_{\rm off}$) via:

$\phi=f_{\rm el}\theta_{\rm el}+f_{\rm par}\psi_{\rm par}+\phi_{\rm off},$ (3)

Cassegrain mounts have $f_{\rm par}=1$ and $f_{\rm el}=0$. Nasmyth-Right type
mounts have $f_{\rm par}=1$ and $f_{\rm el}=+1$ and Nasmyth-Left type mounts
have $f_{\rm par}=1$ and $f_{\rm el}=-1$.

The observed ${\bm{V}}_{mn}$ are modifications of the true $\bar{\bm{V}}_{mn}$
through:

${\bm{V}}_{mn}={\bm{J}}_{m}\bar{\bm{V}}_{mn}{\bm{J}}^{H}_{n},$ (4)

where H is the Hermitian operator.

For circular feeds, the $\bar{\bm{V}}$ are related to the Fourier transforms
of the Stokes parameters ($\tilde{I}$, $\tilde{Q}$, $\tilde{U}$, and
$\tilde{V}$) via

$\bar{\bm{V}}_{mn}\equiv\begin{pmatrix}\mathscr{RR}&\mathscr{RL}\\\
\mathscr{LR}&\mathscr{LL}\end{pmatrix}=\begin{pmatrix}\tilde{I}_{mn}+\tilde{V}_{mn}&\tilde{Q}_{mn}+j\tilde{U}_{mn}\\\
\tilde{Q}_{mn}-j\tilde{U}_{mn}&\tilde{I}_{mn}-\tilde{V}_{mn}\end{pmatrix}.$
(5)

As explained in Paper I, GPCAL assumes that the antenna field-rotation angles
are already corrected at an upstream calibration stage (before performing
global fringe fitting). Thus, Equation 4 becomes:

${\bm{V}}_{mn}={\bm{P}}^{-1}_{m}{\bm{G}}_{m}{\bm{D}}_{m}{\bm{P}}_{m}\bar{\bm{V}}_{mn}{\bm{P}}^{H}_{n}{\bm{D}}^{H}_{n}{\bm{G}}^{H}_{n}({\bm{P}}^{H}_{n})^{-1},$
(6)

##  3 Calibration Procedure

###  3.1 Model Equation

GPCAL assumes that antenna gains are already corrected during the upstream
calibration and imaging/self-calibration procedures. The original GPCAL
pipeline (Park et al., 2021b) fits Equation 6 with the assumption of
${\bm{G}}={\bm{I}}$ to the observed cross-hand visibilities averaged over the
frequency bandwidth. The pipeline assumes that polarimetric leakages are
constant during the observation. However, if this assumption is violated,
there are residual leakages in the data. Following the Appendix in Pesce
(2021), we can write the true leakage matrix ${\bm{D}}$ as:

${\bm{D}}=\begin{pmatrix}1&D^{R}\\\ D^{L}&1\end{pmatrix},$ (7)

and the estimated leakage matrix $\hat{\bm{D}}$ with the assumption of
constant leakage terms during the observation, which can be different from
${\bm{D}}$, as:

$\hat{{\bm{D}}}=\begin{pmatrix}1&D^{R}+\Delta^{R}\\\
D^{L}+\Delta^{L}&1\end{pmatrix}.$ (8)

The visibility matrix after running the GPCAL pipeline would become:

$\displaystyle{\bm{V}}_{mn}$ $\displaystyle=$
$\displaystyle[{\bm{P}}^{-1}_{m}\hat{\bm{D}}_{m}{\bm{P}}_{m}]^{-1}[{\bm{P}}^{-1}_{m}{\bm{D}}_{m}{\bm{P}}_{m}]\bar{\bm{V}}_{mn}$
(9)
$\displaystyle[{\bm{P}}^{H}_{n}{\bm{D}}^{H}_{n}({\bm{P}}^{H}_{n})^{-1}][{\bm{P}}^{H}_{n}\hat{\bm{D}}^{H}_{n}({\bm{P}}^{H}_{n})^{-1}]^{-1}$
$\displaystyle=$
$\displaystyle{\bm{P}}^{-1}_{m}{\bm{R}}_{m}{\bm{P}}_{m}\bar{\bm{V}}_{mn}{\bm{P}}^{H}_{n}{\bm{R}}^{H}_{n}({\bm{P}}^{H}_{n})^{-1},$

where ${\bm{R}}$ is a residual leakage matrix:

$\displaystyle{\bm{R}}$ $\displaystyle\equiv$
$\displaystyle\hat{\bm{D}}^{-1}{\bm{D}}$ (10) $\displaystyle\approx$
$\displaystyle\begin{pmatrix}1&-\Delta^{R}\\\ -\Delta^{L}&1\end{pmatrix}.$

The approximation holds when dropping out second-order terms. With this
approximation, the off-diagonal terms in ${\bm{R}}$ are the “residual” leakage
terms, i.e., $[{\bm{R}}_{i,j}={\bm{D}}_{i,j}-\hat{\bm{D}}_{i,j}]_{\in i\neq
j}$.

First, we use the GPCAL pipeline to remove polarimetric leakages that are
assumed to remain constant over time ($\hat{\bm{D}}$). We then fit Equation 9
to the corrected data in order to derive the “residual” time-dependent
leakages. According to Pesce (2021), it may be more appropriate to redo
calibration entirely with raw data (after data pre-processing) rather than
incrementally calibrating a partially calibrated data set. In spite of this,
we use this two-step procedure because, as we will demonstrate below, the
signal-to-noise ratio of data plays a crucial role in the accurate estimation
of time-dependent leakages. Therefore, before performing time-dependent
instrumental polarization calibration, it is preferable to average the data
over frequency. Nevertheless, as demonstrated in Paper I, modern VLBI arrays
provide a large recording bandwidth, which can lead to significant variations
in the D-terms over frequency. The frequency-dependent D-terms may introduce
non-negligible non-closing errors in the data if they are not removed before
averaging the data over frequency. It is possible to accurately correct
frequency-dependent D-terms using our method presented in Paper I or even the
GPCAL pipeline, which uses averaged data within each IF, which can remove the
gross variations of D-terms across the entire frequency band. Thus, we use the
data after correcting for $\hat{\bm{D}}$, where the gross variations of
D-terms over frequency have already been removed, and correct for the residual
time-dependent leakages.

The cross-hand visibilities in Equation 9 can be re-written as:

$\displaystyle r^{RL}_{mn}$ $\displaystyle\approx$
$\displaystyle(\tilde{Q}_{mn}+j\tilde{U}_{mn})+\Delta^{R}_{m}(t)e^{2j\phi_{m}}r^{LL}_{mn,{\rm
cal}}$ $\displaystyle+$
$\displaystyle\Delta^{L*}_{n}(t)e^{2j\phi_{n}}r^{RR}_{mn,{\rm cal}}$
$\displaystyle+$
$\displaystyle\Delta^{R}_{m}(t)\Delta^{L*}_{n}(t)e^{2j(\phi_{m}+\phi_{n})}(\tilde{Q}_{mn}-j\tilde{U}_{mn})$
$\displaystyle r^{LR}_{mn}$ $\displaystyle\approx$
$\displaystyle(\tilde{Q}_{mn}-j\tilde{U}_{mn})+\Delta^{L}_{m}(t)e^{-2j\phi_{m}}r^{RR}_{mn,{\rm
cal}}$ $\displaystyle+$
$\displaystyle\Delta^{R*}_{n}(t)e^{-2j\phi_{n}}r^{LL}_{mn,{\rm cal}}$
$\displaystyle+$
$\displaystyle\Delta^{L}_{m}(t)\Delta^{R*}_{n}(t)e^{-2j(\phi_{m}+\phi_{n})}(\tilde{Q}_{mn}+j\tilde{U}_{mn}),$

where we replaced $\mathscr{RL}$ and $\mathscr{LR}$ using Equation 5 and
$\mathscr{RR}$ and $\mathscr{LL}$ by the final calibrated parallel-hand
visibilities $r^{RR}_{mn,{\rm cal}}$ and $r^{LL}_{mn,{\rm cal}}$,
respectively. We assume that $\tilde{Q}$ and $\tilde{U}$ are constant over
time during the observation and depend only on the baseline coordinate
$(u,v)$. We fit Equation 3.1 to the data for each scan to derive
$\Delta^{R}(t)$ and $\Delta^{L}(t)$ using the Scipy curve_fit package.

###  3.2 Calibration Strategy

We fit Equation 3.1 to the data for each scan. There is a limited number of
data points used for fitting in this case. In a scan with $N$ antennas, there
are $4N$ free parameters (the real and imaginary parts of the D-terms for RCP
and LCP for each antenna). If one tries to fit a model with $4N$ degrees of
freedom to a limited number of data points within a scan as compared to the
entire dataset, there is likely to be a significant correlation among the
parameters. We will demonstrate below that the best-fit D-terms in this case
have large amplitudes because of the high correlation between parameters,
while the assumed D-terms have small amplitudes in the synthetic data (Section
4). This is similar to performing an amplitude self-calibration at a very
short solution interval on the total intensity data of weak sources, which
results in antenna gain solutions with very high amplitudes.

Due to the limited signal-to-noise ratio of the data (primarily as a result of
the limited number of total data points) compared to the number of free
parameters, it is challenging to accurately constrain the polarimetric
leakages of all stations for each scan. Additionally, antenna pointing should
be reasonably accurate for most stations and most scans. Nevertheless, some
stations may have more inaccurate antenna pointing than others due to poor
weather conditions or large diameters. For stations with inaccurate pointing
during particular scans, the linear polarization models produced after
correcting for on-axis instrumental polarization may differ substantially from
the cross-hand visibilities of all baselines associated with those stations.
If this is the case, one can assume that there are significant residual
leakages ($\Delta$ in Equation 3.1) only for those stations for the scans, and
fix the residual D-terms for the other stations to be zero for fitting. By
doing so, we will be able to avoid the strong correlation between fitting
parameters and improve the accuracy of the fitting222For fitting, the method
uses the visibility weights that are stored in the UVFITS files provided by
the users. However, the users have the option of scaling down the visibility
weights of particular antennas by a constant factor, as was done in the
original GPCAL pipeline (Park et al., 2021b). This feature proves advantageous
for arrays featuring a large variance in sensitivity among the antennas, as it
serves to prevent the fitting process from being dominated by the most
sensitive stations..

The first polarization imaging of M87 from the 2017 EHT observations (Event
Horizon Telescope Collaboration et al., 2021a) was also carried out using a
similar approach. In some calibration pipelines, short baselines are used to
calibrate the D-terms of the stations comprising the short baselines (ALMA and
the Atacama Pathfinder Experiment (APEX) telescope in Chile and the James
Clerk Maxwell Telescope (JCMT) and Submillimeter Array (SMA) on Maunakea in
Hawaii). Those pipelines then assume that the D-terms of those stations have
already been corrected and perform fitting for only the D-terms of the other
stations based on the long baseline data. Through this two-step calibration
strategy, which takes advantage of the fact that short baselines are less
sensitive to complex source structures and have a high signal-to-noise ratio,
it was possible to achieve good calibration accuracy and avoid strong
correlations between the D-terms of many stations.

Accordingly, we adopt the following calibration strategy for time-dependent
calibration of instrumental polarization.

  1. 1.

By using CLEAN with Difmap on the data corrected for the on-axis D-terms,
Stokes $Q$ and $U$ images are produced.

  2. 2.

For each scan and for each station $m$, we compute a norm

$\chi^{2}_{m}=\sum_{n,k}w_{mn,k}\left(|r^{RL}_{mn,k}-\bar{r}^{RL}_{mn,k}|^{2}+|r^{LR}_{mn,k}-\bar{r}^{LR}_{mn,k}|^{2}\right),$
(12)

where $w_{mn,k}$ is the weight of the $k$th visibility matrix $V^{k}$ in the
scan for the baseline between stations $m$ and $n$,
$\bar{r}^{RL}_{mn,k}\equiv\tilde{Q}_{mn,k}+j\tilde{U}_{mn,k}$ and
$\bar{r}^{LR}_{mn,k}\equiv\tilde{Q}_{mn,k}-j\tilde{U}_{mn,k}$ are the model
cross-hand visibilities corresponding to the $k$th visibility for the baseline
derived from the Fourier Transforms of the Stokes $Q$ and $U$ CLEAN images.

  3. 3.

Identify $l$ number of stations having the largest norm values (i.e., the
stations giving the largest and the second largest $\chi_{m}^{2}$ for $l=2$).
$l$ is a free parameter controlled by the variable timecal_freepar in GPCAL,
which is provided by users. The determination of the optimal value for
timecal_freepar may vary depending on the data and the SNR of the source, as
indicated in Table 1 and Section 4. In the case of VLBA data, it is
recommended to adopt a conservative approach and set a small value, such as
timecal_freepar = 1–2.

  4. 4.

We fit Equation 3.1 to the cross-hand visibilities in the scan. Let the
residual D-terms ($\Delta^{R}$ and $\Delta^{L}$) for the stations identified
in Step 3 be free parameters and those for the other stations be fixed to be
zero during the fitting.

  5. 5.

Repeat steps 1–4 for all scans.

  6. 6.

Correct for the derived residual D-terms by inverting the Jones matrices in
Equation 9.

  7. 7.

Iterate the above steps as many times as specified by users. Utilize the
results obtained in Step 6 to generate Stokes $Q$ and $U$ images during Step
1, which are subsequently employed in the remaining calibration steps. The
method yields the reduced $\chi^{2}$ of the fit at every iteration, which can
serve as a valuable tool for selecting the appropriate number of iterations.

While this procedure is similar to the instrumental polarization self-
calibration procedure implemented in the original GPCAL pipeline (Park et al.,
2021b), it differs in that fitting is performed only for leakages of stations
that exhibit large residuals between the model and cross-hand visibilities.
Using synthetic data sets, we will demonstrate the effectiveness of this
strategy in correcting for time-dependent leakages of large amplitudes
(Section 4). The method is implemented in GPCAL, which is publicly available
at https://github.com/jhparkastro/gpcal.

##  4 Validation using Synthetic Data

To validate the performance of the method, we used synthetic data derived from
real VLBA data observed on 2018 Dec 08 at 43 GHz as part of the VLBA-BU-BLAZAR
monitoring program333https://www.bu.edu/blazars/VLBAproject.html. Data
calibration and analysis methods are described in Paper I. Stokes $Q$ and $U$
images were produced with CLEAN in Difmap for 10 sources; 3C 279, 3C 273, OJ
287, 3C 454.3, 1633+382, 3C 345, CTA 102, 3C 84, MKN 501, and 0235+164. The
sources provide a wide range of total intensity and linear polarization
structures and cover a wide range of signal-to-noise ratios.

Figure 1: Comparison of the time-dependent, ground-truth D-term components
(real and imaginary parts) assumed in the synthetic data sets shown on x-axis
with the reconstructed D-term components derived by GPCAL shown on y-axis for
each scan. Each station’s result is presented using distinct colors. For each
source, a correlation is shown between the true and estimated D-terms in
percentage units. Over all antennas, the $L_{1}\equiv|D_{\rm Truth}-D_{\rm
Recon}|$ norm is averaged over the real and imaginary components of the
D-terms for RCP and LCP.  Figure 2: PDF of $D_{\rm Truth}-D_{\rm Recon}$ for
each source shown in each panel. The red solid lines represent Gaussian
distributions fitted to the PDFs. In each panel, the standard deviation
($\sigma$) of the best-fit Gaussian distribution is indicated in percentage
units. The average S/N of the parallel-hand ($\langle SNR\rangle_{\rm par}$)
and cross-hand visibilities ($\langle SNR\rangle_{\rm cro}$) for each source
is denoted in each panel.  Table 1: Goodness of reconstruction of time-
dependent leakages by changing the timecal_freepar parameter

timecal_freepar | 3C 279 | 3C 273 | OJ 287 | 3C 454.3 | 1633+382 | 3C 345 | CTA 102 | 3C 84 | MKN 501 | 0235+164
---|---|---|---|---|---|---|---|---|---|---
1 | 0.17 | 0.17 | 0.17 | 0.19 | 0.14 | 0.19 | 0.19 | 0.15 | 1.28 | 0.30
2 | 0.17 | 0.16 | 0.16 | 0.22 | 0.15 | 0.17 | 0.18 | 0.17 | 1.08 | 0.40
5 | 0.21 | 0.18 | 0.16 | 0.22 | 0.15 | 0.17 | 0.22 | 0.22 | 1.14 | 0.65
10 | 0.50 | 0.76 | 2.16 | 0.84 | 0.81 | 1.64 | 1.56 | 0.73 | 6.37 | 7.27

Note. — Standard deviations of Gaussian distributions fitted to PDFs of
$D_{\rm Truth}-D_{\rm Recon}$ are shown in units of % depending on the
timecal_freepar parameter. Different columns are used to display the results
from different sources.

The synthetic data for these 10 sources were generated using GPCAL as
described in Park et al. (2021a). GPCAL generates synthetic data sets using
Equation 6 and the $I$, $Q$, and $U$ CLEAN models as ground truth source
structures (i.e., $\tilde{I}$, $\tilde{Q}$, $\tilde{U}$ in Equation 5),
assuming $\tilde{V}=0$. Based on the uncertainties of each real data point, we
added thermal noise to the corresponding synthetic data point. Afterwards, the
data was averaged over the entire frequency band in order to increase the
signal-to-noise ratio. We assumed unity antenna gains, i.e.,
${\bm{G}}={\bm{I}}$. We introduced on-axis D-terms, which are randomly chosen
based on the D-term distribution estimated by GPCAL from the real data, and
are assumed to remain constant over time. We then introduced off-axis D-terms
that vary randomly from scan to scan and whose real and imaginary components
follow a Gaussian distribution with a zero mean and a certain standard
deviation. As explained in Section 3 and will be shown in Section 5, the off-
axis D-terms for most stations are anticipated to have small amplitudes in
most cases. However, some stations that experience large pointing errors may
have large amplitudes in their off-axis D-terms. In order to reflect this
realistic situation in the synthetic data set, we assumed a standard deviation
of 0.02 and 0.01 for the VLBA PT and NL station’s leakages, respectively, and
a standard deviation of 0.002 for the other station’s leakages444Obtaining the
precise standard deviation values from actual data is a non-trivial task due
to the dependence of off-axis D-term amplitudes on sources and scans. To
address this issue, we rely on assuming values that would result in off-axis
D-term amplitudes that are reasonably comparable to those obtained from the
actual data.. The method aims to identify and remove off-axis D-terms with
large amplitudes from the data.

First, we ran the GPCAL pipeline on the synthetic data set to correct for the
on-axis D-terms. Based on the assumption that the source 3C 84 is unpolarized,
the pipeline estimates the initial D-terms using the weakly polarized source
3C 84 (with the degree of linear polarization $\lesssim 1\%$ at 43 GHz; Kim et
al. 2019). As a next step, it implements 10 iterations of instrumental
polarization self-calibration using bright calibrators with moderate or high
linear polarization: 3C 279, 3C 273, 3C 454.3, and OJ 287. The data for 3C 84
is included as well in this procedure under the assumption that it is
unpolarized. For all sources, the derived on-axis D-terms are removed from the
data.

Using the method, we corrected for residual time-dependent polarimetric
leakages for all simulated sources. We used timecal_freepar = 2. This means
that during each iteration of calibration, the D-terms of the two stations
with the largest $\chi^{2}_{m}$ values for each scan are free parameters,
while the D-terms of the other stations are fixed to be zero (see Steps 2 and
3 in Section 3). The time-dependent leakage calibration procedure was repeated
ten times.

We present in Figure 1 a comparison of the real and imaginary components of
the time-dependent D-terms assumed in the synthetic data generation (x-axis)
with the reconstructed D-terms estimated by GPCAL (y-axis). The reconstructed
D-terms are the sum of the on-axis and off-axis terms. Each panel indicates
the average $L_{1}\equiv|D_{\rm Truth}-D_{\rm Recon}|$ norm. Our method can
reconstruct the true D-terms accurately, except for MKN 501 and 0235+164,
which have a low signal-to-noise ratio.

Figure 2 shows the probability density function (PDF) of the difference
between the reconstructed and truth D-term components for each source. As
indicated in each panel, we fitted a Gaussian function to each PDF and derived
its standard deviation ($\sigma$). For most sources, $\sigma$ is close to
0.2%, which is the standard deviation of a Gaussian distribution used to
generate random time-dependent D-term components, except for NL and PT.
Consequently, the method can correct for time-dependent D-terms with large
amplitudes from some stations, but cannot accurately constrain time-dependent
D-terms with small amplitudes. A linear polarization image’s quality is
primarily limited by the former D-terms, while the latter is naturally
expected in most realistic scenarios. Thus, the synthetic data test validates
the effectiveness of the method, which is primarily designed to correct for
time-dependent leakages with large amplitudes that usually appear at certain
stations and/or scans.

For the weak sources MKN 501 and 0235+164, $\sigma$ was substantially greater
than 0.2%. We denote the average S/N of the parallel-hand and cross-hand
visibilities for each source in each panel of Figure 2. It is evident that the
S/Ns of the cross-hand visibilities for these sources are notably
low555However, we note that the S/Ns of the parallel-hand visibilities are
also crucial for the precise reconstruction of time-dependent leakages. In our
method, we assume that $\mathscr{RR}$ and $\mathscr{LL}$ are equivalent to the
self-calibrated visibilities $r^{RR}_{\rm cal}$ and $r^{LL}_{\rm cal}$,
respectively, as expressed in Equation 3.1. This is the underlying reason why
the reconstruction for 0235+164 exhibits superior results compared to that of
MKN 501, despite their cross-hand SNRs being similar to each other. This
result demonstrates that we cannot derive accurate solutions due to the
limited signal to noise ratio for these sources despite the small number of
parameters available for each round of fitting. Thermal noise would likely
limit the quality of linear polarization images of these sources, and
systematic errors such as time-dependent leakages would not have a significant
impact. As a result, it is not recommended to perform time-dependent leakage
calibration on weak sources.

Probably the most critical parameter for the method would be timecal_freepar,
since high values of this parameter would result in large correlations between
the fitting parameters and overfitting. As a demonstration of this effect, we
applied the method by changing this parameter and presented the sigma values
for each source from each run in Table 1. When timecal_freepar = 2, the sigma
values are the lowest for most of the bright sources. For weak sources MKN 501
and 0235+164, the values tend to become smaller as timecal_freepar decreases.
With timecal_freepar = 10, i.e., when the D-terms of all 10 stations are free
parameters, the values become very high for all sources. Using small
timecal_freepar values is recommended for VLBA data at 43 GHz, though the
optimal parameter may differ depending on the quality and array of the data.

##  5 Application to Real Data

###  5.1 Data & Analysis

Using two real VLBA data sets, we evaluate the effectiveness of the method.
One is the same data analyzed in Paper I and in Section 4 (project code:
BM462M) observed on 2018 Dec 08. In the period between August 2018 and October
2019, the VLBA pointing model was not optimized due to problems that occurred
during the transition of the telescope control computers (Blanchard, 2021).
Due to this, the quality of the VLBA data observed during this period was
affected, particularly at high frequencies, where the antenna beam size is
small. The inaccurate antenna pointing model also affected the BM462M data we
analyzed. As a consequence, off-axis D-terms will be pronounced and will
significantly limit the quality of linear polarization images derived from
this data.

A second set of data is from the M87 monitoring program observed using the
VLBA simultaneously at 4.7, 5.1, 6.5, and 6.8 GHz in order to investigate the
kinematics of the M87 jet (Park et al. 2019b; Park et al. in prep.). We have
selected the data observed on 2021 July 02 at 4.7 GHz (project code: BP249H).
The VLBA Los Alamos (LA) station reported rain near the end of the
observation, and the Stokes $Q$ and $U$ CLEAN models of all sources were not
able to fit the data well after the rain began. Weather conditions may affect
the leakage characteristics of this station, which motivates us to test our
method.

Data postcorrelation was performed with AIPS and hybrid imaging with CLEAN and
self-calibration in Difmap. On-axis D-terms as well as frequency-dependent
leakages were corrected using GPCAL, as described in Paper I. Following this,
we corrected for time-dependent leakages based on our method, using the
parameter timecal_freepar = 2 and iterating the calibration procedure 10
times.

Figure 3: The derived off-axis D-terms for the BM462M data as a function of
time. In the diagram, the filled circles and open squares represent the
D-terms for RCP and LCP, respectively. The upper and lower panels in each
panel display the real and imaginary parts of the results for each antenna.
Different rows of results are presented for different pairs of adjacent
sources in the sky (3C 273 and 3C 279 in the top, 3C 454.3 and CTA 102 in the
middle, and 3C 345 and 1633+382 in the bottom). In each figure, the averaged
norm $L_{1}\equiv|D_{i}-D_{j}|$ between the neighboring scans (separated by
less than 10 minutes) for the adjacent sources $i$ and $j$ is shown for each
polarization. The symbol $\Delta d$ denotes the distance between adjacent
sources in the celestial sphere. While we provide results for only three
antennas as an illustration, the findings for other antennas were also
qualitatively comparable.  Figure 4: Amplitudes of gain correction factors
(upper) and off-axis D-terms for the VLBA Pie Town station for different pairs
of adjacent sources in the BM462M data set as a function of time. The filled
circles and open squares represent the results for RCP and LCP, respectively.
Figure 5: Off-axis D-term amplitudes as a function of gain correction factors
for the BM462M data set. RCP and LCP are represented by filled circles and
open squares, respectively. Different stations are represented by different
colors.

###  5.2 Results

Figure 3 shows the off-axis leakages as a function of time for the BM462M
data. Although we present the results for three antennas only as an example,
the results for other antennas were also qualitatively similar. The primary
purpose of this test is to compare the trend of the off-axis D-terms derived
from sources located near each other in the sky. The accuracy of the antenna
pointing model depends on the direction of the antenna and the time of day. As
a result, nearby scans of nearby sources may also experience similar levels of
antenna pointing errors. We can expect to observe similar trends in off-axis
D-terms for these scans if the direction-dependent leakages and inaccurate
antenna pointing model are the primary causes of time-dependent leakages. We
compare the derived off-axis D-terms between 3C 273 and 3C 279 (separated by
10.4 degrees in the sky), 3C 454.3 and CTA 102 (6.8 degrees), and 3C 345 and
1633+382 (2.2 degrees).

For all antennas, the trends of the derived off-axis D-terms from nearby
sources are very similar. In each panel, we present the average $L_{1}$ norm
of the D-term components between nearby scans (within ten minutes) of these
source pairs. We obtained smaller average $L_{1}$ norms for closer pairs of
sources, which is reasonable since antenna pointing offsets, and therefore
off-axis D-terms, are direction-dependent.

Figure 6: Same as Figure 3 but for the BP249H data. After UT 3, rain was
reported at the VLBA Los Alamos station (LA; middle panel), following which
the off-axis D-term amplitudes increased substantially.

The method did not attempt to derive time-dependent D-terms for the Mauna Kea
(MK) station for 3C 273 (the top middle panel in Figure 3). It is due to the
fact that the $\chi^{2}_{m}$ value for this station is always low during any
calibration round. The 3C 273 jet displays a complex linear polarization
structure with a moderate to high level of linear polarization (e.g., Attridge
et al., 2005; Hada et al., 2016; Park et al., 2021b). Thus, the cross-hand
visibility of the long baselines associated with MK station has low amplitudes
and low signal-to-noise ratios (at a $\lesssim 0.1$ Jy level for this data).
The effects of systematic errors, such as time-dependent leakages, are not
significant for these baselines, and therefore GPCAL does not attempt to model
them.

Figure 4 illustrates the amplitudes of the gain-correction factors (1/$|G|$)
and the amplitudes of the derived off-axis D-terms as a function of time for
the VLBA Pie Town station. As a result of the inaccurate antenna pointing
model, this station exhibits large gain correction factors and off-axis
D-terms. Both quantities for nearby source pairs exhibit very similar trends
as expected. In addition, the trends between the gain correction factors and
the off-axis D-terms are also very similar. Due to the direction-dependent
nature of antenna leakage, the effects of off-axis D-terms become more
prominent for large antenna pointing offsets.

Figure 5 shows the amplitudes of the derived off-axis D-terms as a function of
the gain correction factors for all scans and all sources. Regardless of the
antenna, there is a strong positive correlation between the two quantities. It
is demonstrated that time-dependent leakages are present in the data due to
the direction-dependent D-terms and imperfect pointing offsets. Their
amplitudes can be very large for scans that are affected by large pointing
offsets.

In Figure 6, we present the derived off-axis D-terms as a function of time for
3C 273 and M87, separated by 10.3 degrees in the sky, for the BP249H data. In
line with the results obtained for BM462M data, the off-axis D-terms for these
sources also exhibit similar trends. As the antenna pointing is accurate at
this low frequency (4.7 GHz), the amplitudes of the off-axis D-terms are
generally small. For LA station, however, large off-axis D-terms are observed
after around 3 UT, when rainfall began. Both 3C 273 and M87 exhibited this
trend, indicating that severe weather had an impact on the large off-axis
D-terms. Observations at other frequencies have also yielded similar results,
although they are not included in the present paper.

Figure 7: Linear polarization images of 3C 279 (upper) and 3C 345 (lower) from
the BM462M data set. The hot color indicates the total intensity distribution
and the colored ticks indicate the electric vector position angles (EVPAs),
with the color scales indicating Ricean de-biased linearly polarized
intensity. Presented in the left and right panels are images obtained before
and after time-dependent leakage correction. There is a note in each panel
indicating the linear polarization dynamic range, which is defined as the
ratio of the peak linear polarization intensity to the linear polarization
off-source rms noise. In the bottom right corner, the shape of the synthesized
beam is shown.

Figure 8: Same as Figure 7 but for the BP249H data for 3C 273 (upper) and M87
(lower).

###  5.3 Evaluation

Figure 7 illustrates the linear polarization images of 3C 279 and 3C 345
derived from the BM462M data before (left) and after (right) time-dependent
leakage correction. A Ricean de-biasing correction has been applied to the
images (i.e., $P_{\rm corr}=P_{\rm obs}\sqrt{1-(\sigma_{P}/P_{\rm obs})}$),
where $P_{\rm corr}$ and $P_{\rm obs}$ denote the corrected and original
linearly polarized intensities (Wardle & Kronberg, 1974; Lister et al., 2018),
respectively, and $\sigma_{P}\equiv(\sigma_{Q}+\sigma_{U})/2$ is the average
of the rms noise in the off-source regions in Stokes $Q$ and $U$ images
(Hovatta et al., 2012). In each image, we note the polarization dynamic range,
which is defined as the ratio of the peak linear polarization intensity to
$\sigma_{P}$. For each image, linear polarization vectors are drawn for
linearly polarized intensities exceeding three times the $\sigma_{P}$ value.

Through the use of the time-dependent leakage correction method, the noise
level of the polarization images has been significantly reduced, resulting in
an increase in polarization dynamic range by factors of two to three. After
the time-dependent leakage correction, linear polarization emission was
observed in the jet regions of weak total intensity that was not evident in
the original linear polarization images. In the case of 3C 279, which has a
high linearly polarized flux ($\sim 0.8$ Jy), the noise level in the
polarization images is more dominated by systematic errors than by thermal
noise. Using the method, residual systematic errors in the data are removed,
making it possible to detect the weak polarization emission in the jet.

A comparison of linear polarization images of 3C 273 and M87 based on the
BP249H data with and without time-dependent leakage correction is presented in
Figure 8. Because of the improvement in polarization dynamic range, weaker
linear polarization emission was detected in both sources, similar to the
results obtained from the BM462M data.

##  6 Conclusions

In this series of papers, we present new methods for calibrating frequency and
time-dependent leakages in VLBI data, which have been implemented in GPCAL. A
wide fractional bandwidth is provided by modern VLBI arrays. The instrumental
polarization is affected by chromatic effects since telescope systems are
usually designed to produce the most optimal polarization response at a given
nominal frequency (e.g., Martí-Vidal et al., 2021). In Paper I, we present a
method for correcting leakages that vary in frequency. It is also possible for
instrumental polarization to change over time. The reason for this is that
instrumental polarization is a direction-dependent effect (e.g., Smirnov,
2011). Therefore, instrumental polarization can vary depending on the accuracy
of the antenna pointing. Antenna pointing is subject to uncertainty, and can
be adversely affected by the weather, such as strong winds and the deformation
of antennas caused by sunlight, which may result in instrumental polarization
that is time-dependent. In poor weather conditions and with large dishes, this
effect is more pronounced. The purpose of this paper is to introduce a method
for correcting time-dependent leakages.

This method works on the data after correcting the “on-axis” D-terms, which
are calculated assuming that the D-terms are constant during observation. In
this manner, “residual” leakages are derived. It calculates how well the
Stokes $Q$ and $U$ images fit the cross-hand visibilities for each antenna and
for each scan based on the data. Next, it determines which antennas exhibit
the greatest difference between the data and the model. By assuming that the
D-terms of all other antennas are zero, the method determines the best-fit
D-terms for these antennas666We plan to implement another calibration strategy
that will correct for off-axis leakages only for scans and stations that are
affected by large antenna pointing offsets in the future. It is based on the
tight correlation between the antenna gain correction factors and the
amplitudes of the derived off-axis D-terms (Figures 4 and 5). This strategy
may reduce the degrees of freedom in the fitting process, but it requires good
images of the source’s total intensity.. The procedure is iterated until the
solutions are converged, and the number of iterations can be controlled by the
user.

The method was tested using a synthetic data set containing time-dependent
leakages generated from real VLBA data observed at 43 GHz. Leakages consist of
two components: the on-axis, stable D-terms drawn from the distribution of the
D-terms from the real data, as well as the off-axis, variable D-terms that are
randomly selected from Gaussian distributions for each scan. It is assumed
that two antennas have large standard deviations (0.02 and 0.01 for PT and NL
stations, respectively) for the Gaussian distributions and other antennas have
small standard deviations (0.002), in order to mimic some recently observed
VLBA data sets at 43 GHz.

Based on the method, we derived time-dependent leakages for the synthetic data
set. We found that the method was able to successfully reconstruct the time-
dependent leakages that were assumed in the generation of synthetic data. In
most cases, the distribution of the difference between the estimated and
ground-truth D-term components (the real and imaginary parts) could be
described by Gaussian distributions with standard deviations less than 0.002.
This result indicates that the method is capable of capturing the time-
dependent leakages of large amplitudes, which are the dominant factors in
limiting the quality of linear polarization images. The purpose of developing
the method is actually to correct for the largest systematic errors in the
data by using prior information that the off-axis D-term amplitudes of most
antennas in realistic VLBI arrays are usually small.

Our method was applied to two sets of real data obtained with the VLBA at 4.7
and 43 GHz. In the VLBA 4.7 GHz data, significant deviations are observed
between the source polarization models and cross-hand visibilities for all
baselines associated with the VLBA Los Alamos station after the time the
station reported rain. An antenna control computer upgrade caused inaccurate
antenna pointing models in the VLBA 43 GHz data. We derived the off-axis
D-terms for several pairs of sources located close to one another in the sky
as a function of time. The trends of the derived off-axis D-terms between
nearby sources are very similar. Off-axis D-terms for nearby scans of less
separated sources displayed a higher degree of consistency. Additionally, the
amplitudes of the off-axis D-terms and the gain-correction factors show very
similar trends. Based on these results, it can be concluded that the method is
capable of capturing time-dependent leakages caused by direction-dependent
leakages and imperfect antenna pointing.

For all sources analyzed in this paper, the dynamic ranges of linear
polarization images have been improved by factors of $\gtrsim 2$ after
correcting for time-dependent leakages using the method. According to these
results, the method can significantly enhance the polarization image fidelity
of VLBI data. Additionally, the method will be useful for improving global
VLBI arrays operating at millimeter wavelengths, such as the GMVA and EHT,
which have very small antenna beams because of the high observing frequencies
and large dishes.

We express our gratitude to the anonymous referee for conducting a
comprehensive review of our manuscript, which greatly enhanced the quality of
the paper. J.P. acknowledges financial support through the EACOA Fellowship
awarded by the East Asia Core Observatories Association, which consists of the
Academia Sinica Institute of Astronomy and Astrophysics, the National
Astronomical Observatory of Japan, Center for Astronomical Mega-Science,
Chinese Academy of Sciences, and the Korea Astronomy and Space Science
Institute. This work is supported by the Ministry of Science and Technology of
Taiwan grant MOST 109-2112-M-001-025 and 108-2112-M-001-051 (K.A). The VLBA is
an instrument of the National Radio Astronomy Observatory. The National Radio
Astronomy Observatory is a facility of the National Science Foundation
operated by Associated Universities, Inc.

## References

  * Asada et al. (2002) Asada, K., Inoue, M., Uchida, Y., et al. 2002, PASJ, 54, L39, doi: 10.1093/pasj/54.3.L39
  * Attridge et al. (2005) Attridge, J. M., Wardle, J. F. C., & Homan, D. C. 2005, ApJ, 633, L85, doi: 10.1086/498392
  * Blanchard (2021) Blanchard, J. 2021, VLBA Test Memo 73 (Socorro: NRAO)
  * Broderick et al. (2020) Broderick, A. E., Gold, R., Karami, M., et al. 2020, ApJ, 897, 139, doi: 10.3847/1538-4357/ab91a4
  * Casadio et al. (2017) Casadio, C., Krichbaum, T., Marscher, A., et al. 2017, Galaxies, 5, 67, doi: 10.3390/galaxies5040067
  * Chael et al. (2018) Chael, A. A., Johnson, M. D., Bouman, K. L., et al. 2018, ApJ, 857, 23, doi: 10.3847/1538-4357/aab6a8
  * Chael et al. (2016) Chael, A. A., Johnson, M. D., Narayan, R., et al. 2016, ApJ, 829, 11, doi: 10.3847/0004-637X/829/1/11
  * Cotton (1993) Cotton, W. D. 1993, AJ, 106, 1241, doi: 10.1086/116723
  * Event Horizon Telescope Collaboration et al. (2019a) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019a, ApJ, 875, L1, doi: 10.3847/2041-8213/ab0ec7
  * Event Horizon Telescope Collaboration et al. (2019b) —. 2019b, ApJ, 875, L2, doi: 10.3847/2041-8213/ab0c96
  * Event Horizon Telescope Collaboration et al. (2019c) —. 2019c, ApJ, 875, L3, doi: 10.3847/2041-8213/ab0c57
  * Event Horizon Telescope Collaboration et al. (2019d) —. 2019d, ApJ, 875, L4, doi: 10.3847/2041-8213/ab0e85
  * Event Horizon Telescope Collaboration et al. (2019e) —. 2019e, ApJ, 875, L5, doi: 10.3847/2041-8213/ab0f43
  * Event Horizon Telescope Collaboration et al. (2019f) —. 2019f, ApJ, 875, L6, doi: 10.3847/2041-8213/ab1141
  * Event Horizon Telescope Collaboration et al. (2021a) Event Horizon Telescope Collaboration, Akiyama, K., Algaba, J. C., et al. 2021a, ApJ, 910, L12, doi: 10.3847/2041-8213/abe71d
  * Event Horizon Telescope Collaboration et al. (2021b) —. 2021b, ApJ, 910, L13, doi: 10.3847/2041-8213/abe4de
  * Gabuzda (2018) Gabuzda, D. 2018, Galaxies, 7, 5, doi: 10.3390/galaxies7010005
  * Gómez et al. (2022) Gómez, J. L., Traianou, E., Krichbaum, T. P., et al. 2022, ApJ, 924, 122, doi: 10.3847/1538-4357/ac3bcc
  * Greisen (2003) Greisen, E. W. 2003, Astrophysics and Space Science Library, Vol. 285, AIPS, the VLA, and the VLBA, ed. A. Heck, 109, doi: 10.1007/0-306-48080-8_7
  * Hada et al. (2016) Hada, K., Kino, M., Doi, A., et al. 2016, ApJ, 817, 131, doi: 10.3847/0004-637X/817/2/131
  * Hamaker (2000) Hamaker, J. P. 2000, A&AS, 143, 515, doi: 10.1051/aas:2000337
  * Hamaker & Bregman (1996) Hamaker, J. P., & Bregman, J. D. 1996, A&AS, 117, 161
  * Hamaker et al. (1996) Hamaker, J. P., Bregman, J. D., & Sault, R. J. 1996, A&AS, 117, 137
  * Hovatta et al. (2012) Hovatta, T., Lister, M. L., Aller, M. F., et al. 2012, AJ, 144, 105, doi: 10.1088/0004-6256/144/4/105
  * Jones (1941) Jones, R. C. 1941, Journal of the Optical Society of America (1917-1983), 31, 488
  * Jorstad et al. (2017) Jorstad, S. G., Marscher, A. P., Morozova, D. A., et al. 2017, ApJ, 846, 98, doi: 10.3847/1538-4357/aa8407
  * Kettenis et al. (2006) Kettenis, M., van Langevelde, H. J., Reynolds, C., & Cotton, B. 2006, in Astronomical Society of the Pacific Conference Series, Vol. 351, Astronomical Data Analysis Software and Systems XV, ed. C. Gabriel, C. Arviset, D. Ponz, & S. Enrique, 497
  * Kim et al. (2019) Kim, J. Y., Krichbaum, T. P., Marscher, A. P., et al. 2019, A&A, 622, A196, doi: 10.1051/0004-6361/201832920
  * Leppanen et al. (1995) Leppanen, K. J., Zensus, J. A., & Diamond, P. J. 1995, AJ, 110, 2479, doi: 10.1086/117706
  * Lisakov et al. (2021) Lisakov, M. M., Kravchenko, E. V., Pushkarev, A. B., et al. 2021, ApJ, 910, 35, doi: 10.3847/1538-4357/abe1bd
  * Lister et al. (2018) Lister, M. L., Aller, M. F., Aller, H. D., et al. 2018, ApJS, 234, 12, doi: 10.3847/1538-4365/aa9c44
  * Martí-Vidal et al. (2021) Martí-Vidal, I., Mus, A., Janssen, M., de Vicente, P., & González, J. 2021, A&A, 646, A52, doi: 10.1051/0004-6361/202039527
  * Matthews et al. (2018) Matthews, L. D., Crew, G. B., Doeleman, S. S., et al. 2018, PASP, 130, 015002, doi: 10.1088/1538-3873/aa9c3d
  * Napier (1989) Napier, P. J. 1989, in Astronomical Society of the Pacific Conference Series, Vol. 6, Synthesis Imaging in Radio Astronomy, ed. R. A. Perley, F. R. Schwab, & A. H. Bridle, 39
  * Park & Algaba (2022) Park, J., & Algaba, J. C. 2022, Galaxies, 10, doi: 10.3390/galaxies10050102
  * Park et al. (2021a) Park, J., Asada, K., Nakamura, M., et al. 2021a, ApJ, 922, 180, doi: 10.3847/1538-4357/ac26bf
  * Park et al. (2021b) Park, J., Byun, D.-Y., Asada, K., & Yun, Y. 2021b, ApJ, 906, 85, doi: 10.3847/1538-4357/abcc6e
  * Park et al. (2019a) Park, J., Hada, K., Kino, M., et al. 2019a, ApJ, 871, 257, doi: 10.3847/1538-4357/aaf9a9
  * Park et al. (2018) Park, J., Kam, M., Trippe, S., et al. 2018, ApJ, 860, 112, doi: 10.3847/1538-4357/aac490
  * Park et al. (2019b) Park, J., Hada, K., Kino, M., et al. 2019b, ApJ, 887, 147, doi: 10.3847/1538-4357/ab5584
  * Pesce (2021) Pesce, D. W. 2021, AJ, 161, 178, doi: 10.3847/1538-3881/abe3f8
  * Sault et al. (1996) Sault, R. J., Hamaker, J. P., & Bregman, J. D. 1996, A&AS, 117, 149
  * Shepherd (1997) Shepherd, M. C. 1997, in Astronomical Society of the Pacific Conference Series, Vol. 125, Astronomical Data Analysis Software and Systems VI, ed. G. Hunt & H. Payne, 77
  * Smirnov (2011) Smirnov, O. M. 2011, A&A, 527, A106, doi: 10.1051/0004-6361/201016082
  * Thum et al. (2008) Thum, C., Wiesemeyer, H., Paubert, G., Navarro, S., & Morris, D. 2008, PASP, 120, 777, doi: 10.1086/590190
  * Tiede (2022) Tiede, P. 2022, The Journal of Open Source Software, 7, 4457, doi: 10.21105/joss.04457
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: https://doi.org/10.1038/s41592-019-0686-2
  * Wardle & Kronberg (1974) Wardle, J. F. C., & Kronberg, P. P. 1974, ApJ, 194, 249, doi: 10.1086/153240
  * Zhao et al. (2022) Zhao, G.-Y., Gomez, J. L., Fuentes, A., et al. 2022, arXiv e-prints, arXiv:2205.00554.  https://arxiv.org/abs/2205.00554


# Flat or crumpled: states of active symmetric membranes

Sudip Mukherjee  sudip.bat<EMAIL_ADDRESS>Barasat Government
College, 10, KNC Road, Gupta Colony, Barasat, Kolkata 700124, West Bengal,
India  Abhik Basu  abhik.123<EMAIL_ADDRESS>Theory Division,
Saha Institute of Nuclear Physics, 1/AF Bidhannagar, Calcutta 700064, West
Bengal, India

###### Abstract

We set up and study the hydrodynamic theory for active fluid and tethered
membranes. We focus on those membranes which are inversion-symmetric. We show
that such membranes are either described by appropriate linear hydrodynamic
equations, which are exact in the asymptotic long wavelength limit, giving
stable flat phases with positional quasi long range orders, or be linearly
unstable, implying crumpling of the membranes. We argue that in a such an
active membrane thermal noises dominate over any active noises, and use it to
calculate the correlation functions of the undulation and in-plane
displacements of the membrane in the stable case, and the associated
correlation functions of the embedding bulk flow velocities.

Biomembranes or cell membranes at small scales have complex structures alberts
. Nonetheless, equilibrium properties of single-component membranes (i.e.,
those which are made of only one type of lipid molecules) at scales much
larger than the typical molecular scales can be described by continuum
theories, parametrised by a set of elastic constants nelson ; safran . For
instance, an arbitrarily large tensionless fluid membrane that has vanishing
shear modulii, cannot remain orientationally ordered or statistically flat at
non-zero temperature $T$, no matter how low $T$ is nelson ; luca ; tirtha-pre
. Instead, it is always crumpled at any finite $T$ nelson ; luca ; tirtha-pre
, a phenomenon known as thermal crumpling of fluid membranes. This can be
explained by invoking the Mermin-Wagner-Hohenberg (MWHT) theorem mwt ; mwt1
.In contrast, a tethered membrane, which has finite shear modulii, can display
orientational LRO at sufficiently low but finite $T$ nelson ; kantor ;
lubensky . This apparent breakdown of MWHT is attributed to the phonon-
mediated effective long range interactions in a tethered membrane, that take
it out of the scope of MWHT.

Due to the underlying microscopic active processes, e.g., nonequilibrium
ﬂuctuations of cell cytoskeletons active-cyto  and active proteins in lipid
membranes active-protein , live cell membranes display unusual dynamical
fluctuation properties fluc1 ; fluc2 ; fluc3 . Similarly, experimental
comparisons of fluctuations in healthy, live red blood cell membranes and
those in ATP (adenosine triphosphate) depleted red blood cells suggest that
active processes should be responsible in live red blood cell membrane
fluctuations rbc-exp1 ; rbc-exp2 ; rbc-exp3 . More recently, by comparing
membrane response and fluctuations on single red blood cells, a violation of
the Fluctuation–Dissipation-Theorem chaikin  is reported in Ref. joanny-
natphys , that clearly demonstrates the nonequilibrium process underlying the
flickering of live red blood cell joanny-natphys . The generality of
experimental observations on the manifestations of nonequilibrium nature of
cell membrane dynamics has prompted scientists to look for minimal physical
descriptions. Such theoretical attempts have led to theories of active
membranes for a various settings joanny-natphys ; nir ; sriram . We here
consider an open, inversion-symmetric, isotropic single-component fluid or
tethered membrane, which is tensionless in the equilibrium limit and embedded
in a three-dimensional (3D) isotopic active fluid, and construct hydrodynamic
equations for the such systems. Our results complement the existing theories.
Our theoretical predictions are general and can be tested in artificial in-
vitro settings.

Our principal results are as follows: (i) The nonequilibrium dynamics of
nearly flat symmetric active membranes, fluid or tethered, emerge from their
couplings with an embedding active isotropic bath, provided by a 3D active but
isotropic fluid. (ii) The linearised hydrodynamic equations predict that such
a nearly flat patch can be stable, corresponding to orientational LRO for a
fluid membrane, or orientational LRO and in-plane translational quasi long
range order (QLRO). The orientational LRO is due to a positive active tension,
generated due to the interplay between active processes and 3D embedding fluid
flows. These results are exact in the asymptotic long wave length limit, with
all nonlinear effects being irrelevant in that limit. (iii) Although these
results are reminiscent of those for an equilibrium fluid or tethered membrane
with a finite tension at some effective temperature, in this asymptotically
exact linear theory the equal-time undulation mode correlation function
depends on a dynamic (active, see below) coefficient, whereas the in-plane
displacement correlators depend only on thermodynamic or equilibrium
coefficients. This indicates the inherent nonequilibrium nature of the system.
(iv) Even though the linearised hydrodynamic equations are exact at the
largest scales, the undulation modes can be linearly unstable at the
intermediate scales, if the membrane is permeable and the associated
permeation flow is destabilising. (v) The active tension can also be negative,
in which case the undulation modes are linearly unstable, signalling membrane
crumpling. Interestingly, the in-plane displacement modes are still linearly
stable, but should get nonlinearly unstable, possibly due to their nonlinear
couplings with the undulation modes, which are irrelevant in the linearly
stable case.

We now construct the hydrodynamic equations for a permeable active symmetric
tethered membrane embedded in a 3D bulk isotopic active fluid containing
active orientable particles. The hydrodynamic equation for an active symmetric
fluid membrane can be obtained by simply suppressing the in-plane displacement
modes. The conformation ${\bf R}({\bf s})$ of any membrane in 3D is
parametrised by a two-dimensional (2D) position vector ${\bf s}$ that
designates points in the membrane. Imposing that all the points in the
membrane labelled by $\bf s$ retain their coordinates, i.e., a convected
coordinate system is used convec1 ; convec2 , the local membrane velocity
${\bf V}_{m}({\bf s})\equiv\partial_{t}{\bf R}({\bf s})$, which formally
defines the equation of motion of the membrane. For a nearly flat patch of a
membrane without any overhangs, it is convenient to use the Monge gauge nelson
; chaikin , in which ${\bf R}=({\bf r}+{\bf u},h({\bf r},t))$, where $h({\bf
r},t)$ is the height above points ${\bf r}=(x,y)$ on a reference plane; see
Fig. 1 for a schematic diagram of the membrane. Here, ${\bf u}$ is the
relative in-plane displacement of a point whose location in the undistorted
membrane is at ${\bf R}=({\bf r},h=0)$. Notice that this reference plane is
completely arbitrary; no substrate is present to provide a fixed preferred
reference horizontal plane in this problem.

Figure 1: Schematic model diagram of a membrane depicted in the Monge gauge
that is embedded in an active isotropic fluid. Here, ${\bf R}=({\bf r}+{\bf
u},h({\bf r},t))$, where $h({\bf r},t)$ is the height above points ${\bf
r}=(x,y)$ on a reference plane. Further, $\hat{n}$ is the local normal to the
membrane.

For the ease of notations, let us denote the 3D bulk hydrodynamic flow
velocity by $v^{\prime}_{\alpha},\,\alpha=x,y,z$ and the 2D flow velocity
along the membrane by $v_{i},\,i=x,y$. Ignoring any number conservation, $h$
and ${\bf u}\equiv u_{i}$ are the two hydrodynamic variables in the problem
along with the total momentum of the system. In the presence of a permeation
flow velocity $v_{\text{perm}}$ locally normal to the membrane, considered as
a permeable fluid film cai , the equation for $h$ reads

$\frac{\partial h}{\partial t}-v_{\text{perm}}=v_{z}^{\prime}|_{z=h}.$ (1)

In general, $v_{\text{perm}}$ is a function of $h$ and $u_{i}$ that is
inversion-symmetric, tilt-invariant due to the rotational symmetry nelson  and
also invariant under $u_{i}\rightarrow u_{i}+const.$. For small fluctuations,
expanding in $h$ and $u_{i}$, we can phenomenologically write

$v_{\text{perm}}=\tilde{\mu}\nabla^{2}h+\mu_{1}\nabla^{2}h({\bm{\nabla}}h)^{2}+\mu_{2}{\bm{\nabla}}\cdot{\bf
u}\nabla^{2}h+\mu_{p}\frac{\delta{\cal F}}{\delta h},$ (2)

where $\mu_{p}>0$ is an equilibrium permeation coefficient,
$\tilde{\mu},\,\mu_{1},\,\mu_{2}$ are all active permeation coefficients, each
of which can be positive or negative. Here, $\cal F$ is the Landau-Ginzburg
free energy functional of a tensionless tethered membrane. General symmetry
considerations dictate the form of $\cal F$:

${\cal F}=\frac{1}{2}\int d^{2}r\,[K(\nabla^{2}h)^{2}+\lambda u_{ii}+2\mu
u_{ij}u_{ij}].$ (3)

Here,

$u_{ij}=\frac{1}{2}(\partial_{i}u_{j}+\partial_{j}u_{i}+\partial_{i}h\,\partial_{j}h)$
(4)

is the strain tensor nelson ; chaikin , $K>0$ is the bend modulus and
$\lambda,\mu$ are the Lamé coefficients. The equation of $u_{i}$ is

$\frac{\partial u_{i}}{\partial t}=v_{i}.$ (5)

In the extreme Stokesian limit, dropping the fluid inertia,
$v_{\alpha}^{\prime}$ can be obtained from the force balance condition
$\nabla_{\beta}\sigma^{\prime}_{\alpha\beta}=0$, where stress tensor
$\sigma_{\alpha\beta}^{\prime}$ is given by

$\sigma^{\prime}_{\alpha\beta}=\frac{\eta^{\prime}}{2}(\partial_{\alpha}v^{\prime}_{\beta}+\partial_{\beta}v^{\prime}_{\alpha})-P^{\prime}\delta_{\alpha\beta}$
(6)

for an incompressible bulk fluid landau-fluid , where $\eta^{\prime}$ is the
shear viscosity and $P^{\prime}$ is the 3D pressure. The 2D in-plane velocity
$v_{i}$ satisfies the Stokes equation

$\nabla_{j}\sigma^{\text{tot}}_{ij}+f_{i}=0.$ (7)

Here, 2D stress tensor
$\sigma^{\text{tot}}_{ij}=\eta(\partial_{i}v_{j}+\partial_{j}v_{i})/2+\eta_{b}{\bm{\nabla}}\cdot{\bf
v}\delta_{ij}-P\delta_{ij}+\sigma_{ij}^{a}$, $\eta$ and $\eta_{b}$ are the 2D
shear and bulk visocities respectively and $P$ is a 2D partial pressure. The
last contribution to $\sigma^{\text{tot}}_{ij}$ comes from the active stress
$\sigma_{\alpha\beta}^{a}=An_{\alpha}n_{\beta}$, where $\hat{\bf
n}=(-\hat{z},\nabla_{i}h)$ is the local normal to the leading order in
$h$-fluctuations with $\hat{z}$ being the unit vector along the $z$-direction.
This active stress can arise, e.g., in a membrane with actin filaments grafted
normally to it. The actin polymersation or depolymerisation can force the
membrane. Now define a polarisation vector ${\bf p}$ that is non-vanishing
only at the location of the membrane. Elsewhere in the bulk active fluid being
in its isotropic phase, we set ${\bf p}=0$. For a nearly flat membrane, the
condition of normal anchoring gives $p_{\alpha}=-\hat{n}_{\alpha}$ niladri-
epje3 . A non-zero $p_{i}$ on the membrane allows us to define an active
stress $\sim Ap_{\alpha}p_{\beta}$ sriram-RMP , which with normal anchoring
gives $\sigma^{a}_{ij}$ as defined above. Lastly, $f_{i}$ is the force on the
2D membrane that originates from the shear stress of the bulk flow:

$f_{i}=\eta^{\prime}\big{(}\frac{\partial v_{i}^{\prime}}{\partial
z}+\frac{\partial v_{z}^{\prime}}{\partial
x_{i}}\big{)}|_{z=h_{+}}-\eta^{\prime}\big{(}\frac{\partial
v_{i}^{\prime}}{\partial z}+\frac{\partial v_{z}^{\prime}}{\partial
x_{i}}\big{)}|_{z=h_{-}}.$ (8)

Furthermore, the condition of normal stress balance at the membrane yields

$(2\eta^{\prime}\frac{\partial v_{z}^{\prime}}{\partial
z}-P^{\prime})|_{z=h_{+}}-(2\eta^{\prime}\frac{\partial
v_{z}^{\prime}}{\partial z}-P^{\prime})|_{z=h_{-}}=\frac{\delta{\cal
F}}{\delta h}+\nabla_{j}\sigma^{a}_{zj}.$ (9)

We impose $v_{\alpha}^{\prime}(z=\pm\infty)=0$ as the boundary conditions on
$v_{\alpha}^{\prime}$. Lastly, the no-slip boundary condition on the membrane
implies $v^{\prime}_{\alpha=i}(z=h_{+})=v^{\prime}_{\alpha=i}(z=h_{-})=v_{i}$
and the continuity condition
$v^{\prime}_{z}(z=h_{+})=v_{z}^{\prime}(z=h_{-})$. We now solve for
$v^{\prime}_{\alpha}({\bf r},z)$ subject to the above conditions. It is
convenient to first define $v_{\alpha}^{\prime}({\bf q},z,t)$ and $v_{i}({\bf
q},t)$ as the in-plane Fourier transforms of $v_{\alpha}({\bf r},z,t)$ and
$v_{i}({\bf r},t)$ [see, e.g., Refs. brochard ; kramer ; niladri-epje4 ]. We
obtain

$\frac{\partial h({\bf q},t)}{\partial t}=\frac{Aq}{4\eta^{\prime}}h({\bf
q},t)+v_{\text{perm}}-\frac{1}{4\eta^{\prime}q}\frac{\delta{\cal F}}{\delta
h}.$ (10)

Using (2) and (3) and ignoring nonlinearities, we then get

$\frac{\partial h({\bf q},t)}{\partial t}=\frac{Aq}{4\eta^{\prime}}h({\bf
q},t)-\tilde{\mu}q^{2}h({\bf q},t)-\frac{Kq^{3}}{4\eta^{\prime}}h({\bf q},t),$
(11)

which is the linearised hydrodynamic equation for the undulation modes.The
intermediate calculational details, similar to Refs. sriram ; niladri-epje4 ,
are straight forward, and are given in Supplemental Material (SM) sm .

To obtain the hydrodynamic equation for $\bf u$, we decompose ${\bf
u}=(u_{\parallel},\,u_{\perp})$ and ${\bf v}=(v_{\parallel},\,v_{\perp})$,
where $\parallel$ and $\perp$ refer to in-plane directions, which are parallel
and perpendicular to $\hat{\bf q}$. To calculate $v_{\perp}$ and
$v_{\parallel}$, we use (7) and resolve it perpendicular and parallel to
$\hat{\bf q}$ giving $v_{\perp}({\bf q},t)=-\mu
q/(2\eta^{\prime})\,u_{\perp}({\bf q},t)$ to the linear order in (small)
fluctuations. Then

$\frac{\partial u_{\perp}({\bf q},t)}{\partial t}=-\frac{\mu
q}{2\eta^{\prime}}\,u_{\perp}({\bf q},t)$ (12)

is the linearised hydrodynamic equation for $u_{\perp}$. In the long
wavelength limit, the equation for $v_{\parallel}({\bf q},t)$ reads

$4\eta^{\prime}qv_{\parallel}({\bf
q},t)=-iqP-(\lambda+2\mu)q^{2}u_{\parallel}({\bf q},t).$ (13)

Here, $P$, the “2D partial pressure” is given by $P=\chi(\rho)$, where $\rho$
is the concentration of the active elements in the 2D layer. Now assume the
total number of these active elements $\int d^{2}r\,\rho$ is constant only on
average due to possible non-number conserving events, e.g., birth-death of the
active elements, $\rho$ is not a conserved density and hence is not a
hydrodynamic variable. In fact, fluctuations $\delta\rho$ about the mean
concentration $\rho_{0}$ relaxes fast. In the spirt of coarse-grained
hydrodynamic appraoches, we then replace $\chi(\rho)$ by $\chi(\rho_{0})$, a
constant. This ensures that $P$ drops out of (13). Solving for
$v_{\parallel}({\bf q},t)$ we get

$\frac{\partial u_{\parallel}({\bf q},t)}{\partial
t}=-\frac{\lambda+2\mu}{4\eta^{\prime}}qu_{\parallel}({\bf q},t)$ (14)

to the linear order in (small) fluctuations.

Linear stability:- We now study the linear stability of the hydrodynamic
equations (11), (12) and (14). Clearly, Eqs. (12) and (14) are always linearly
stable, since $\mu>0$ and $\lambda+2\mu>0$ for the thermodynamic stability of
${\cal F}$. In contrast, Eq. (11) is linearly stable or unstable in the long
wavelength limit, when the active parameter $A$ is positive or negative. Even
when $A>0$ corresponding to linear stablity, if $\tilde{\mu}<0$ and there is a
sufficient scale separation between $A$ and $\tilde{\mu}$ on one hand, and
$\tilde{\mu}$ and $K$ on the other, $h$-dynamics displays linear instability
for an intermediate range of wavevectors. This instability is due to the
destabilising nature of the active permeation fow. See Fig. 2 for a schematic
plot of $\lambda_{h}\equiv
Aq/(4\eta^{\prime})-\tilde{\mu}q^{2}+Kq^{3}/(4\eta^{\prime})$ versus $q$,
where $\lambda_{h}$ is the effective damping in (11).

Figure 2: Schematic plot of $\lambda_{h}(q)$ versus $q$ showing an unstable
intermediate band of wavevectors. Within the range $q_{q}<q,q_{2}$,
$\lambda_{h}(q)>0$, implying linear instability.

Correlation functions in the linearised hydrodynamics:- To calculate the
correlation functions, we now must include noises in the governing dynamical
equations. Let $g_{h},\,g_{\perp}$ and $g_{\parallel}$ be the noises added in
Eqs. (11), (12) and (14), all of which are assumed to be Gaussian-distributed
with zero mean. The noises have two sources - thermal and active noises. The
thermal noises that survive in the equilibrium limit of the dynamics is
controlled by the FDT chaikin , which relates the noise variances to the
damping coefficients. This gives brochard

$\displaystyle\langle g_{h}({\bf q},t)g_{h}({\bf-q},0)\rangle$
$\displaystyle=$ $\displaystyle\frac{2D_{h}\delta(t)}{4\eta^{\prime}q},$ (15)
$\displaystyle\langle g_{a}({\bf q},t)g_{\perp}({\bf-q},0)\rangle$
$\displaystyle=$ $\displaystyle\frac{2D_{a}\delta(t)}{2\eta^{\prime}q},$ (16)

where $a=\parallel,\,\perp$; $D_{h},\,D_{\perp}$ and $D_{\parallel}$ are all
positive and scale with $T$. Unsurprisingly, all the noises are long range due
to the long range hydrodynamic damping present even in the equilibrium limit
of the dynamics. The active noises are expected to be short range frank-motor
; casimir ; madan-jacques-PRL , and hence subdominant to the thermal noises in
the long wavelength limit. We thus conclude that although the model is active,
the noises are essentially thermal. This forms a major conclusion of this
work. Do we then have the FDT being valid in the linearised hydrodynamics? To
know this, we must calculate the correlation function of the fields, which can
be done by using (15)-(16). We consider the linearly stable case with no
intermediate wavevector instabilities. In the long wavelength limit, we get

$\displaystyle\langle|h({\bf q},t)|^{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{D_{h}}{Aq^{2}},\langle|u_{a}({\bf
q},t)|^{2}\rangle=\frac{D_{a}}{\tilde{a}q^{2}}.$ (17)

Here, $a=\parallel,\,\perp$ and correspondingly
$\tilde{a}=\mu,\,\lambda+2\mu$, respectively. Now if we compare (17) with the
known equilibrium results for a fluid membrane with a finite tension to be
identified with $A$, and a 2D crystal, we identify
$D_{h}=D_{\perp}=D_{\parallel}=k_{B}T$, where $k_{B}$ is the Boltzmann
constant. Nonetheless, we must notice that $A$ is a dynamical coefficient (an
active coefficient actually) and not a thermodynamic coefficient, unlike $\mu$
and $\lambda$. The fact that $A$ appears in a static or time-independent
quantity $\langle|h({\bf q},t)|^{2}\rangle$ [see Eq. (17) above], is a
reminder to us that in spite of the apparent equilibrium-like behaviour the
underlying dynamics is actually a nonequilibrium dynamics. Lastly, the forms
of (17) reveal that all the fields have dynamic exponent $z=1$ and roughness
or spatial scaling exponent $\chi=0$ in 2D, the physically relevant dimension.
Here, $z=1$ means the fluctuations relax much faster than that in diffusive
relaxations.

In the linearly stable case, all of $\langle h^{2}({\bf r},t)\rangle,\,\langle
u_{\perp}^{2}({\bf r},t)\rangle,\,\langle u_{\parallel}^{2}({\bf r},t)\rangle$
scale as $\ln(L/a_{0})$, giving positional QLRO, where $L$ is the linear
system size and $a_{0}$ is a microscopic cut-off. Likewise, equal-time
correlation functions $\langle h({\bf r},t)h(0,t)\rangle,\,\langle
u_{\perp}({\bf r},t)u_{\perp}(0,t)\rangle,\,\langle u_{\parallel}({\bf
r},t)u_{\parallel}(0,t)$ all behave as $\ln(r/a_{o})$, growing logarithmically
with $r\equiv|{\bf r}|$, another hallmark of QLRO. Furthermore,
$\langle({\bm{\nabla}}h)^{2}\rangle$ is independent of $L$ for large $L$,
giving orientational LRO.

In the unstable case ($A<0$), clearly both $\langle h^{2}({\bf r},t)\rangle$
and $\langle({\bm{\nabla}}h)^{2}\rangle$ diverge as soon as $L$ exceeds a
threshold determined by the zero of $\lambda_{h}$ for a finite $q$. This
implies crumpling of the membrane. This is different from the well-known
thermal crumpling phenomenon of equilibrium fluid membranes at any finite $T$.
Thermal crumpling occurs in tensionless fluid membranes, whereas in the
present case it is due to a negative “active tension” $A<0$. This brings us to
another question: what happens if $A=0$? Since $A$ can be varied by varying
the strength of the active stress, let us assume it has been tuned to zero.
Even in this case, there would be no thermal crumpling, assuming
$\tilde{\mu}>0$ giving

$\langle h^{2}({\bf r},t)\rangle\sim
L,\;\langle({\bm{\nabla}}h)^{2}\rangle\sim\text{independent of}\,L,$ (18)

corresponding to positional short range order and orientational long range
order. While with $A=0$, there is no active tension and the membrane is
tensionless, still $\tilde{\mu}>0$ means the effective wavevector-dependent
bend modulus $K_{e}\sim 1/q$, diverging strongly for small $q$. This renders
thermal crumpling ineffective.

Nonlinear effects:- We now discuss the nonlinear effects so far neglected in
our studies above. There are two distinct sources of nonlinearities in the
present model - those that exist in equilibrium and appear from $\cal F$ and
another class which appear in the nonequilibrium dynamics and have their
origin in the underlying active processes. Some of these nonlinear terms have
a geometric origin arising from the definition of the strain tensor (4), while
some others have purely dynamical (active) origin. These nonlinear terms
appear both in the dynamics of $h$ and ${\bf u}$; see SM for details. However,
scaling and renormalisation arguments can be used to show that all these
nonlinear terms are irrelevant in the linearly stable states. This means the
scaling of the correlation functions in the linearly stable case as reported
above are in fact exact in the asymptotic long wavelength limit.

The nonlinearities are however expected to be important in two other
situations: (i) With $A<0$, the linearised $h$-dynamics given by (11) is
linearly unstable. However, the dynamics of $u_{\perp},\,u_{\parallel}$ remain
linealy stable, as both $\mu$ and $\lambda+2\mu$ are assumed to be always
positive. Thus the undulation modes get linearly unstable, but the in-plane
displacements remain linearly stable, a situation that is physically
unacceptable. After all, with the undulation modes getting bigger in the
linearly unstable regime, the strain (4) must also get bigger, eventually
affect the in-plane displacements. This apparent contradiction between this
physically acceptable expectation and the linear stability of (12) and (14)
even when $A<0$ can be reconciled by the nonlinear effects. As shown in SM,
these nonlinearities in (12) and (14) involve $\nabla_{i}h$, which also gets
bigger with $A<0$, which in turn makes fluctuations in ${\bf u}$ growing
leading to nonlinear instabilities. (ii) The exponential growth of $h$ in the
unstable case with $A<0$ ultimately must saturate due to the nonlinear
effects. Studies of this is outside the scope of the present work.

Bulk velocity solutions:- As the membrane fluctuates, it stirs the embedding
bulk fluid, giving bulk flows. This can be measured, which calls for
evaluation of $v_{\alpha}$. Solving the 3D Stokes equation for $z>0,\,z<0$, we
find sm

$\displaystyle\langle|v^{\prime}_{\perp}({\bf q},z,t)|^{2}\rangle$
$\displaystyle=$ $\displaystyle\frac{\mu^{2}}{4\eta^{\prime
2}}q^{2}\langle|u_{\perp}({\bf q},t)|^{2}\rangle\exp(-2qz),$ (19)
$\displaystyle\langle|v^{\prime}_{z}({\bf q},z,t)|^{2}\rangle$
$\displaystyle=$ $\displaystyle\frac{A^{2}q^{2}}{16\eta^{\prime
2}}(1-qz)^{2}\langle|h({\bf q},t)|^{2}\rangle$ (20)
$\displaystyle+q^{2}z^{2}\frac{(\lambda+2\mu)^{2}}{16\eta^{\prime 2}}$

for $z>0$. Correlation of $v^{\prime}_{\parallel}({\bf q},z,t)$ can be
calculated from those of $v^{\prime}_{z}({\bf q},z,t)$.Thus the correlations
of $v_{\alpha}$ are expressed in terms of the correlations of $h$ and $\bf u$.
Since the latter correlations diverge in the unstable case, the bulk flow
should also accordingly stirred strongly. This can be measured by tracking
passive tracer particles, which in turn should give valuable information about
the membrane itself.

Summary and outlook:- We have presented a comprehensive hydrodynamic theory of
a nearly flat active tethered membrane immersed in a bulk isotopic active
fluid. The activity can make the system linearly stable or unstable. We make
quantitative, experimentally testable predictions on the nature of order and
scaling of the in-plane and out of plane fluctuations by using the linearised
hydrodynamic theory, which we argue to be exact in the asymptotic long
wavelength limit. In the unstable case, the membrane is shown to be crumpled.
Our theory can be tested in in-vitro set up involving, say, a graphene sheet
coated with identical lipid layers on both sides and immersed in an active
fluid bath (e.g., a solution of live orientable bacteria or actin filaments)
in its isotropic phase.

Acknowledgement:- S.M. and A.B. thank the SERB, DST (India) for partial
financial support through the TARE scheme [file no.: TAR/2021/000170].

## References

  * (1) B Alberts, D Bray, J Lewis, M Raff, K Roberts and J D Watson, Molecular Biology of the Cell, 3rd edn (New York: Garland, 1994).
  * (2) Statistical Mechanics of Membranes and Surfaces, edited by D. Nelson, T. Piran, and S. Weinberg (World Scientific, Singapore, 1989).
  * (3) S. A. Safran, Statistical Thermodynamics of Surfaces, Interfaces, and Membranes (Westview, Boulder, 2003).
  * (4) L. Peliti and S. Leibler, Effects of Thermal Fluctuations on Systems with Small Surface Tension, Phys. Rev. Lett. 54, 1690 (1985).
  * (5) T. Banerjee and A. Basu, Thermal fluctuations and stiffening of symmetric heterogeneous fluid membranes, Phys. Rev. E 91, 012119 (2015).
  * (6) N. D. Mermin and H. Wagner, Absence of Ferromagnetism or Antiferromagnetism in One-or Two-Dimensional Isotropic Heisenberg Models, Phys. Rev. Lett. 17, 1133 (1966).
  * (7) P. C. Hohenberg, Existence of long-range order in one and two dimensions, Phys. Rev. 158, 383 (1967).
  * (8) Y. Kantor, M. Kardar, and D. R. Nelson, Tethered surfaces: Statics and dynamics, Phys. Rev. A 35, 3056 (1987).
  * (9) J. A. Aronovitz and T. C. Lubensky, Fluctuations of Solid Membranes, Phys. Rev. Lett. 60, 2634 (1988).
  * (10) E. Frey and D. R. Nelson, Dynamics of flat membranes and flickering of red blood cells, J. Phys. I France 1, 1715 (1991).
  * (11) D. Mizuno, C. Tardin, C., C. F. Schmidt and F. C. MacKintosh, Nonequilibrium mechanics of active cytoskeletal networks, Science 315 370 (2007); Stuhrmann B et al, Phys. Rev. E 86, 020901 (2012).
  * (12) M-J Huang, H-Y Chen and A S Mikhailov 2012 Eur. Phys. J. E 35 119 (2012); M-J Huang, R. Kapral, A. S. Mikhailov and H-Y Chen, J. Chem. Phys. 138 195101 (2013).
  * (13) H. G. Döbereiner, B. J. Dubin-Thaler, J. M. Hofman, H. S. Xenias, T. N. Sims, G. Giannone, M. L. Dustin, C. H. Wiggins, and M. P. Sheetz, Lateral Membrane Waves Constitute a Universal Dynamic Pattern of Motile Cells, Phys. Rev. Lett. 97, 038102 (2006).
  * (14) C. H. Chen, F. C. Tsai, C. C. Wang, and C. H. Lee, Three-Dimensional Characterization of Active Membrane Waves on Living Cells, Phys. Rev. Lett. 103, 238101 (2009).
  * (15) J. Allard and A. Mogilner, Traveling waves in actin dynamics and cell motility, Curr. Opin. Cell Biol. 25, 107 (2013).
  * (16) T. Betz, M. Lenz, J.-F. Joanny and C. Sykes, ATP-dependent mechanics of red blood cells. Proc. Natl Acad. Sci. USA 106, 15320 (2009).
  * (17) Y. Park et al, Metabolic remodeling of the human red blood cell membrane, Proc. Natl Acad. Sci. USA 107, 1289 (2010).
  * (18) R. Rodríguez-García et al, Direct cytoskeleton forces cause membrane softening in red blood cells Biophys. J. 108, 2794 (2015).
  * (19) H. Turlier, D. A. Fedosov, B. Audoly, T. Auth, N. S. Gov, C. Sykes, J.-F. Joanny, G. Gompper and T. Betz, Equilibrium physics breakdown reveals the active nature of red blood cell flickering, Nature Phys. 12, 513 (2016).
  * (20) R. Shlomovitz and N. S. Gov, Membrane Waves Driven by Actin and Myosin, Phys. Rev. Lett. 98, 168103 (2007).
  * (21) A. Maitra, P. Srivastava, M. Rao and S. Ramawasmy, Activating Membranes, Phys. Rev. Lett. 112, 258101 (2014).
  * (22) T. banerjee and A. Basu, Active processes make mixed lipid membranes either flat or crumpled, New J. Phys. 20, 013028 (2018).
  * (23) R. Aris, Vectors, Tensors, and the Basic Equations of Fluid Mechanics (Dover, USA, 1989).
  * (24) D. Hu, P. Zhang and Weinan E, Continuum theory of a moving membrane, Phys. Rev. E 75, 041605 (2007).
  * (25) P. M. Chaikin and Tom C Lubensky, Principles of condensed matter physics, (Cambridge university press Cambridge, 2000).
  * (26) W. Cai and T. C. Lubensky, Hydrodynamics and dynamic fluctuations of fluid membranes, Phys. Rev. E 52, 4251 (1995).
  * (27) L. D. Landau and E. M. Lifshitz, Fluid Mechanics (Butterworth and Heinemann, Oxford, 1998).
  * (28) N. Sarkar and A. Basu, Generic instabilities in a ﬂuid membrane coupled to a thin layer of ordered active polar ﬂuid, Eur. Phys. J E 36, 86 (2013).
  * (29) M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Hydrodynamics of soft active matter, Rev. Mod. Phys. 85, 1143 (2013).
  * (30) F. Brochard and J. F. Lennon, Frequency spectrum of the flicker phenomenon in erythrocytes, J. Phys. (Paris) 36, 1035 (1975).
  * (31) L. Kramer, Theory of Light Scattering from Fluctuations of Membranes and Monolayers, J. Chem. Phys. 55, 2097 (1971).
  * (32) N. Sarkar and A. Basu, Role of interfacial friction for flow instabilities in a thin polar-ordered active fluid layer, Phys. Rev. E 92, 052306 (2015).
  * (33) Supplemental Material.
  * (34) F. Jülicher, A. Ajdari and J. Prost, Modeling molecular motors, Rev. Mod. Phys. 69, 1269 (1997).
  * (35) A. Basu, J.-F. Joanny, F. Jülicher and J. Prost, Casimir stresses in active nematic ﬁlms, New J. Phys. 21 123046 (2019).
  * (36) J.-F. Rupprecht, A. Singh Vishen, G. V. Shivashankar, M. Rao, and J. Prost, Phys. Rev. Lett. 120, 098001 (2018).

##  Appendix A Derivation of the hydrodynamic equations

Assume a nearly flat membrane to be spread out along the XY plane and for
concreteness $z=0$ be its average position. We consider small fluctuations.
The three-dimensional (3D) flow field $v_{\alpha}^{\prime}$ satisfies the 3D
Stokes equation $\nabla_{\beta}\sigma^{\prime}_{\alpha\beta}=0$ giving

$\eta^{\prime}\nabla^{2}v_{\alpha}^{\prime}-\nabla_{\alpha}P^{\prime}=0,$ (21)

where $\nabla_{\beta}\sigma_{\alpha\beta}^{\prime}$ is given in the main text.
Similarly, in the Stokesian limit, the 2D flow field $v_{i}$ satisfies
$\nabla_{j}\sigma_{ij}^{\text{tot}}+f_{i}=0$, where $f_{i}$ is the shear force
due to the bulk flow on the 2D flow in the membrane.

We solve $v_{\alpha}^{\prime}$ in the Fourier space. It is convenient to
resolve $v_{\alpha}^{\prime}$ into three components
$v_{z}^{\prime},\,v_{\perp}^{\prime}=(\hat{z}\times\hat{\bf q})\cdot{\bf
v}^{\prime}$ and $v_{\parallel}^{\prime}=\hat{\bf q}\cdot{\bf v}^{\prime}$,
where ${\bf q}$ is an in-plane Fourier wavevector. These three components then
satisfy

$\displaystyle\eta^{\prime}(-q^{2}+\partial_{z}^{2})v_{\perp}^{\prime}$
$\displaystyle=$ $\displaystyle 0,$ (22)
$\displaystyle\eta^{\prime}(-q^{2}+\partial_{z}^{2})v_{\parallel}^{\prime}$
$\displaystyle=$ $\displaystyle iqP^{\prime},$ (23)
$\displaystyle\eta^{\prime}(-q^{2}+\partial_{z}^{2})v_{z}^{\prime}$
$\displaystyle=$ $\displaystyle\frac{\partial P^{\prime}}{\partial z}.$ (24)

Here $\bf q$ is an in-plane Fourier wavevector. Now impose incompressibility
condition on the 3D bulk velocity: ${\bm{\nabla}}\cdot{\bf v}^{\prime}=0$.
This in turn gives $\nabla^{2}P^{\prime}=0$ as the equation for pressure
$P^{\prime}$.

We write the general solutions for $v_{\alpha}^{\prime}$ and $P^{\prime}$ as

$\displaystyle v_{\perp}^{\prime}({\bf q},z,t)$ $\displaystyle=$
$\displaystyle C_{1}({\bf q},t)\exp(-q\,z),\;z>0,$ (25) $\displaystyle=$
$\displaystyle C_{2}({\bf q},t)\exp(q\,z),\;z<0,$ $\displaystyle
v_{z}^{\prime}({\bf q},z,t)$ $\displaystyle=$ $\displaystyle[A_{1}({\bf
q},t)+z\,B_{1}({\bf q},t)]\exp(-q\,z),\,z>0,$ (26) $\displaystyle=$
$\displaystyle[A_{2}({\bf q},t)+z\,B_{2}({\bf q},t)]\exp(q\,z),\,z<0,$
$\displaystyle P^{\prime}({\bf q},z,t)$ $\displaystyle=$ $\displaystyle
D_{1}({\bf q},t)\exp(-q\,z),\;z>0,$ (27) $\displaystyle=$ $\displaystyle
D_{2}({\bf q},t)\exp(q\,z),\;z<0.$

Lastly, $v_{\parallel}^{\prime}$ can be solved by using the 3D
incompressibility condition giving

$v_{\parallel}^{\prime}=\frac{i}{q}\frac{\partial v_{z}^{\prime}}{\partial
z}.$ (28)

Using the continuity of the 3D bulk flow velocity at $z=h\approx 0$, we get

$C_{1}=C_{2},\;A_{1}=A_{2},\;\text{and}\;-qA_{1}+B_{1}=qA_{2}+B_{2}.$ (29)

Then using he $z$-component of the 3D Stokes equation, we find

$D_{1}=-2\eta^{\prime}B_{1},\;D_{2}=-2\eta^{\prime}B_{2}.$ (30)

Next, we find

$(2\eta^{\prime}\frac{\partial v_{z}^{\prime}}{\partial
z}-P^{\prime})|_{z=0_{+}}-(2\eta^{\prime}\frac{\partial
v_{z}^{\prime}}{\partial
z}-P^{\prime})|_{z=h_{-}}-4\eta^{\prime}qA_{1}=-4\eta^{\prime}qv_{z}(z=h).$
(31)

Then, $\nabla_{j}\sigma_{zj}^{a}=A\nabla_{j}n_{z}n_{j}=-A\nabla^{2}h$ to the
lowest order in ${\bm{\nabla}}h$. Combining everything, the normal stress
balance condition can be used to obtain

$\frac{\partial h}{\partial t}=v_{z}(z=h)=\frac{Aq}{4\eta^{\prime}}h({\bf
q},t)-\frac{1}{4\eta^{\prime}q}\frac{\delta F}{\delta h},$ (32)

which is the linearised hydrodynamic equation for the undulation mode to the
lowest order in $h$, as already obtained in the main text.

We now derive the hydrodynamic equations for the in-plane displacements
$u_{i}$ to the linear order in ${\bm{\nabla}}h$. It is convenient to resolve
$\bf u$ in components parallel and perpendicular to $\bf q$ in the XY plane:
${\bf u}\equiv(u_{\parallel},\,u_{\perp})$. Resolving the 2D Stokes equation
parallel and perpendicular to $\hat{\bf q}$, we get to to the linear order in
${\bm{\nabla}}h$

$\displaystyle-(\eta+\eta_{b})q^{2}\,v_{\parallel}({\bf
q},t)-iqP-(\lambda+2\mu)q^{2}u_{\parallel}({\bf q},t)=-f_{\parallel}({\bf
q},z=0,t),$ (33) $\displaystyle-\eta q^{2}v_{\perp}({\bf q},t)-\mu
q^{2}u_{\perp}({\bf q},t)=-f_{\perp}({\bf q},z=h,t).$ (34)

We find

$f_{\perp}({\bf q},z=h,t)=\eta^{\prime}[\frac{\partial v_{\perp}^{\prime}({\bf
q},t)}{\partial z}|_{z=h_{+}}-\frac{\partial v_{\perp}^{\prime}({\bf
q},t)}{\partial
z}|_{z=h_{-}}=-\eta^{\prime}(C_{1}+C_{2})q=-2\eta^{\prime}qv^{\prime}_{\perp}=-2\eta^{\prime}qv_{\perp}.$
(35)

Combining, we find in the long wavelength limit

$\frac{\partial u_{\perp}({\bf q},t)}{\partial t}=v_{\perp}({\bf
q},t)=-\frac{\mu q_{\perp}}{2\eta^{\prime}}u_{\perp}({\bf q},t)$ (36)

to the linear order in (assumed small) fluctuations, as obtained in the main
text.

We proceed in the same way to obtain the hydrodynamic equation for
$v_{\parallel}({\bf q},t)$. We find

$f_{\parallel}({\bf q},z=0,t)=\eta^{\prime}[\frac{\partial
v_{\parallel}^{\prime}}{\partial z}|_{z=0_{+}}-\frac{\partial
v_{\parallel}^{\prime}}{\partial
z}|_{z=0_{-}}]=-4\eta^{\prime}qv_{\parallel}({\bf q},z=0,t).$ (37)

This gives in the long wavelength limit

$4\eta^{\prime}qv_{\parallel}({\bf
q},t)+iqP=-(\lambda+2\mu)q^{2}u_{\parallel}({\bf q},t).$ (38)

Here, the 2D partial pressure of the active particles $P=\chi(\rho)$, where
$\rho$ is the active particle concentration. If there are birth-death
processes, then $\rho$ is nont conserved, but there is a mean density
$\rho_{0}$. Number density fluctuations $\delta\rho$ about $\rho_{0}$ are fast
variables, i.e., its relaxation time remains finite even in the limit
$q\rightarrow 0$, and hence is non-hydrodynamic. Therefore, in the spirit of
hydrodynamics, we drop $\delta\rho$ from our theory in the large time, long
wavelength limit and get

$\frac{\partial u_{\parallel}({\bf q},t)}{\partial t}\equiv v_{\parallel}({\bf
q},t)=-\frac{\lambda+2\mu}{4\eta^{\prime}}qu_{\parallel}({\bf q},t),$ (39)

as obtained in the main text.

##  Appendix B Solutions of the bulk flow fields

We use the general solutions of the 3D Stokes equation to express
$v^{\prime}_{\alpha}({\bf q},z,t)$ in terms of $h({\bf
q},t),u_{\parallel}({\bf q},t)$ and $u_{\perp}({\bf q},t)$. We separately
solve in the upper ($z>0$) and lower ($z<0$) half planes. First consider
$z>0$. We find

$v^{\prime}_{\perp}({\bf q},z,t)=C_{1}({\bf q},t)\exp(-qz)=v_{\perp}({\bf
q},t)\exp(-qz)-\frac{\mu}{2\eta^{\prime}}qu_{\perp}({\bf q},t)\exp(-qz),$ (40)

Similarly,

$\displaystyle v^{\prime}_{z}(q,z,t)$ $\displaystyle=$
$\displaystyle[v^{\prime}_{z}(q,z=0,t)+z(-qv^{\prime}_{z}(q,z=0,t)-iqv_{\parallel}(q,t)]\exp(-qz)$
(41) $\displaystyle=$ $\displaystyle\frac{A}{4\eta^{\prime}}h({\bf
q},t)(1+qz)\exp(-qz)+iqz\frac{\lambda+2\mu}{4\eta^{\prime}}qu_{\parallel}({\bf
q},t)\exp(-qz).$

Finally, $v_{\parallel}^{\prime}({\bf q},z,t)=(i/q)\partial
v_{z}^{\prime}({\bf q},z,t)/\partial z$ due to the incompressibility of the
embedding bulk fluid. This gives

$v_{\parallel}^{\prime}(q,z,t)=-iqz\frac{A}{4\eta^{\prime}}h({\bf
q},t)-(q-q^{2}z)\frac{\lambda+2\mu}{4\eta^{\prime}}u_{\parallel}({\bf
q},t)\exp(-qz).$ (42)

To get the solutions of $v_{\alpha}$ in the lower half plane $z<0$, replace
$z$ by $-z$ in the above solutions.

##  Appendix C Nonlinear terms

We now systematically derive the nonlinear terms which are to be added in the
dynamical equations for $h,u_{\parallel}$ and $u_{\perp}$. There are several
sources of possible nonlinearities. We first consider those which are
spatially local. These can come from (i) $v_{\text{perm}}$ as given in the
main text, which has both active and equilibrium contributions, and (ii)
equilibrium relaxation term $1/(4\eta^{\prime}q)\,(\delta{\cal F}/\delta h)$.
To proceed further, we start from the form of the nonlinear stress nelson :

$u_{ij}=\frac{1}{2}(\partial_{i}u_{j}+\partial_{j}u_{i}+\partial_{i}h\partial_{j}h),$
(43)

still neglecting the the contributions quadratic in $\bf u$ nelson . This
gives

$\displaystyle u_{mm}^{2}=({\bm{\nabla}}\cdot{\bf
u})^{2}+{\bm{\nabla}}\cdot{\bf u}({\bm{\nabla}}\cdot{\bf
u})^{2}+\frac{1}{4}({\bm{\nabla}}h)^{4},$ (44) $\displaystyle
u_{ij}u_{ij}=(\partial_{i}u_{j})^{2}+\frac{1}{4}({\bm{\nabla}}h)^{4}+\frac{1}{2}(\partial_{i}u_{j}+\partial_{j}u_{i})\,\partial_{i}h\partial_{j}h,$
(45)

where we have ignored total derivative terms. The terms listed in (45) appear
in the free energy $\cal F$. Then

$\displaystyle\frac{\delta{\cal F}}{\delta
h}|_{\text{nonlinear}}=-(\mu+\frac{\lambda}{2})\nabla_{j}[({\bm{\nabla}}h)^{2}\nabla_{j}h]-\mu\nabla_{j}[(\partial_{i}u_{j}+\partial_{j}u_{i})\partial_{i}h]-\lambda\nabla[({\bm{\nabla}}\cdot{\bf
u})\,{\bm{\nabla}}h].$ (46)

The terms in (46) contribute to the hydrodynamic equation of $h$ through
$v_{\text{perm}}$ as given in the main text and the relaxational term
$1/(4\eta^{\prime}q)\,(\delta{\cal F}/\delta h)$.

Next, we consider the nonlocal nonlinear terms. These will originate from the
solutions of the 3D bulk flow fields beyond the linear order in $h$. See,
e.g., Eq. (31), which we have evaluated by setting $z=h\approx 0$ that
produces the correct results to the linear order in $h$. Going beyond the
linear order results would produce factors of $qh$ and its higher powers
appearing multiplicatively with the various terms in the linear theory.

We now consider the possible nonlinear terms in the $\bf u$-equation. As for
the undulation modes, the nonlinearities can be local or nonlocal. We first
consider the local nonlinear terms. These originate from anharmonic
contributions in the free energy $\cal F$. We find

$\frac{\delta{\cal F}}{\delta
u_{i}}|_{\text{nonlinear}}=-\mu\nabla_{j}[\nabla_{i}h\,\nabla_{j}h]-\frac{\lambda}{2}\nabla_{i}({\bm{\nabla}}h)^{2}.$
(47)

Projecting (47) parallel and perpendicular to $\hat{\bf q}$ in the $XY$ plane
produces the local nonlinear contributions to the hydrodynamic equations for
$u_{\parallel}$ and $u_{\perp}$. Finally, similar to the undulation modes,
there are nonlocal nonlinear contributions to the equations for
$u_{\parallel}$ and $u_{\perp}$, which originate from the solutions of the 3D
ambient flow fields. See, e.g., Eq. (35) above. Going beyond the linear order
in fluctuations produce factors of $qh$ and its higher powers appearing
multiplicatively with the terms in the linear hydrodynamic equations for
$u_{\parallel}$ and $u_{\perp}$. We now find out whether or not these
nonlinear terms can change the predictions from the linear hydrodynamic
equations in the long wavelength limit.

To asses the renormalisation group (RG) relevance of these above nonlinear
terms, we rescale space, time and the fields as

$\displaystyle{\bf r^{\prime}}=b{\bf r},\,t^{\prime}=b^{z}t,\,h({\bf
r},t)=b^{\chi_{h}}h({\bf r^{\prime}},t^{\prime}),\,u_{\perp}({\bf
r},t)=b^{\chi_{\perp}}u_{\perp}({\bf
r^{\prime}},t^{\prime}),\,u_{\parallel}({\bf
r},t)=b^{\chi_{\parallel}}u_{\parallel}({\bf r},t).$ (48)

In the linear theory, $z=1,\,\chi_{h}=\chi_{\perp}=\chi_{\parallel}=0$ in 2D,
such that with these choices for the scaling exponents, the parameters in the
linearised hydrodynamic equations for $h,u_{\perp}$ and $u_{\parallel}$ and
the corresponding noise strengths remain unchanged under rescaling. Simple
power counting however shows that under such a rescaling procedure, all the
nonlinear terms acquire scale factor of $b$ raised to negative powers,
rendering them irrelevant in the long wavelength limit.