# Voice-preserving Zero-shot Multiple Accent Conversion

###### Abstract

Most people who have tried to learn a foreign language would have experienced
difficulties understanding or speaking with a native speaker’s accent. For
native speakers, understanding or speaking a new accent is likewise a
difficult task. An accent conversion system that changes a speaker’s accent
but preserves that speaker’s voice identity, such as timbre and pitch, has the
potential for a range of applications, such as communication, language
learning, and entertainment. Existing accent conversion models tend to change
the speaker identity and accent at the same time. Here, we use adversarial
learning to disentangle accent dependent features while retaining other
acoustic characteristics. What sets our work apart from existing accent
conversion models is the capability to convert an unseen speaker’s utterance
to multiple accents while preserving its original voice identity. Subjective
evaluations show that our model generates audio that sound closer to the
target accent and like the original speaker.

Index Terms—  accent conversion, adversarial learning, voice conversion,
speech synthesis

##  1 Introduction

For many people, their own or others’ accents present a severe obstacle to
communication. For some other people, watching a movie set in America with
British accented speakers creates a dissonance that takes them out of the
immersive experience. Therefore, a system that converts a speaker’s accent yet
still preserves the original voice identity can have a great impact in a wide
range of situations including communication, language learning, and
entertainment.

Conventional accent conversion approaches require the availability of parallel
speech, i.e. utterances reading the same text with target accents, during
training and/or inference [1, 2, 3, 4]. Practically, we often do not have
native and non-native utterances with the same linguistic contents. Usually,
each speaker in a dataset only speaks with one accent, and there is relative
scarcity in the available number of speakers and audio recordings for non-
native accents. Therefore, the main challenge in voice-preserving accent
conversion is the need to disentangle features related to a speaker’s voice,
accent, and linguistic contents.

Previous works have used adversarial learning to disentangle features[5, 6],
where both apply a discriminator to wipe-out speaker dependent information
from content embeddings. Other works, such as [7, 8], achieve disentanglement
via quantization of different features to obscure undesired information.

Our proposed model is most similar to the accent conversion models proposed in
[9] and [6]. Both models use a recognition-synthesis framework, where they use
an automatic speech recognition (ASR) model to extract speaker and accent-
independent linguistic features and a text-to-speech (TTS) model trained on
the speech with target accent to synthesize audio. The model from [9] only
converts Hindi accented speech to American accented speech, and the converted
audios are perceived as different from the voice of original speaker. On the
other hand, our model converts arbitrary accented input speech to speech with
multiple target accents. Listeners from our perceptual tests agree that our
model is good at preserving the original speaker’s voice identity. Our work
also differs significantly from [6], which only works for a few seen speakers,
in that our model performs well for unseen speakers.

In our model, we do not train an ASR model ourselves but utilize a pre-trained
model checkpoint to extract speaker and accent independent text predictions
prior to training. We further disentangle the accent-dependent features from
other features with an accent discriminator. Finally, the processed,
disentangled features are re-combined and fed to a HiFiGAN decoder to
reconstruct the audio waveform. Our work make contributions in three major
ways: 1) To the best of our knowledge, our model is the first to convert
arbitrary accented, unseen speech to multiple target accents while preserving
non-accent related voice characteristics. 2) We do not require text labels
associated with accented speech or speaker ID labels during training, although
we use an existing ASR model checkpoint trained on native accented English
speech to extract linguistic features. 3) We convert accent without modifying
the input audios’ duration so that the input and output audios can be
synchronized, allowing for applications such as dubbing a video with different
accents.

##  2 Proposed System

Our proposed system during training and inference is shown in Fig. 1. Unlike
[6] and [9], we do not train an accented ASR model with text labels
corresponding to accented speech. Instead, we take advantage of an available
wav2vec2.0 checkpoint111Specifically, we downloaded the checkpoint
corresponding to ”Wav2Vec 2.0 Large (LV-60 + CV + SWBD + FSH)” from the link
https://github.com/facebookresearch/fairseq/tree/main/examples/wav2vec. that
has been pre-trained using self-supervised learning and fine-tuned on ASR task
using $960$ hours of LibriSpeech data[10].

Let $\mathbf{x}$ denote the input audio waveform, $\mathbf{\hat{x}}$ the
output audio, and $M$ the function that transforms a waveform into the
corresponding 80-dimensional mel-spectrogram. We train our model to minimize
the reconstruction loss in (1)

$L_{mel}=\mathbb{E}_{\mathbf{x}\in\mathcal{X}}[||M(\mathbf{x})-M(\hat{\mathbf{x}})||_{1}]$
(1)

Fig. 1: (a) Training (b) During inference, we convert accent by feeding in the
target accent ID to the pronunciation encoder.

###  2.1 Pronunciation Encoder

Inspired by many early works that have used phonetic posteriograms (PPG) to
capture accent-dependent features, we use a pronunciation encoder, as shown in
Fig. 2, to synthesize accent-dependent pronunciation sequence given text
predictions and accent ID [2, 11]. For each accent ID, the pronunciation
encoder learns a unique embedding, which is concatenated with every frame of
character-level wav2vec2.0 prediction. The four transformer layers with 8-head
attention mechanism accounts for how context affects pronunciations. Applying
a dropout of $0.3$ promotes the decoder to rely more on the acoustic encoder
for non-accent related voice features.

Fig. 2: The pronunciation encoder.

###  2.2 Acoustic Encoder

The acoustic encoder maps the mel-frequency cepstrum coefficients (MFCC) and
periodicity features to a single 256-dimensional vector. The acoustic encoder
consists of four convolution layers with kernel sizes (5, 3, 3, 1) and
dilations (1, 2, 1, 1), a self-attention layer, and finally an average pooling
layer over time to output a single vector. Instead of using an absolute
positional embedding, we use the convolution layers to act as relative
positional embedding for the self-attention layer as in [10].

We use adversarial training to remove accent information from the output of
the acoustic encoder. All the accents are labeled as either native or foreign.
The accent discriminator, consisting of two fully connected layers, learns to
predict 1 if the source audio accent is native and 0 if foreign, while the
acoustic encoder tries to force the accent discriminator to predict 1 all of
the time. Mathematically, $\mathcal{N}$ the set of audios with native English
accent, $\mathcal{F}$ the set with foreign accents, and
$\mathcal{X}=\mathcal{N}\cup\mathcal{F}$ the entire set of training data. Let
$\mathbf{z}$ be the output of the acoustic encoder. The accent discriminator
(AD) attempts to minimize $L_{AD}$ from (2) while the acoustic encoder tries
to minimize $L_{AD,adv}$ from (3).

$L_{AD}=-\mathbb{E}_{\mathcal{N}}[\log(AD(\mathbf{z}))]-\mathbb{E}_{\mathcal{F}}[\log(1-AD(\mathbf{z}))]$
(2) $L_{AD,adv}=-\mathbb{E}_{\mathcal{F}}[\log(AD(\mathbf{z}))]$ (3)

###  2.3 HiFiGAN-based Voice Decoder

Finally, we re-combine the accent-dependent pronunciation encodings, the
acoustic features, and F0 sequence. We use a modified HiFiGAN to invert the
processed features back to audio waveform [12]. Our modifications to the
HiFiGAN architecture include adding an additional convolution layer with
kernel size 11, modifying the number of input channels to HiFiGAN, and
modifying the upsampling rates so that the output length matches the source
audio length.

We use the same multi-scale discriminator (MSD) and multi-period discriminator
(MPD) as in [12] to encourage the synthesis of natural sounding audio.
Mathematically, if we consider the HiFiGAN discriminators, MSD and MPD, as one
discriminator $HD$, the HiFiGAN discriminators try to minimize $L_{HD}$ in
(4).

$L_{HD}=\mathbb{E}_{\mathcal{X}}[(HD(\mathbf{x})-1)^{2}+(HD(\hat{\mathbf{x}}))^{2}]$
(4)

As in [12], the adversarial loss $L_{HD,adv}$, and the feature mapping loss
functions, $L_{FM}$ applied to the rest of the model in (5) and (6).

$L_{HD,adv}=\mathbb{E}_{\mathcal{X}}(HD(\hat{\mathbf{x}})-1)^{2}]$ (5)
$L_{FM}=\mathbb{E}_{\mathcal{X}}[\sum_{i=1}^{T}\frac{1}{N_{i}}||HD^{i}(\mathbf{x})-HD^{i}(\hat{\mathbf{x}})||]$
(6)

where $HD^{i}$, $N_{i}$ denote the features and the number of features in the
$i^{th}$ layer of the HiFiGAN discriminators.

##  3 Experiments

###  3.1 Data Sets

Our training data is summarized in Table 1. Our training data set includes $8$
different accents: American (AM), Arabic (AR), British (BR), Hindi (HI),
Korean (KO), Mandarin (MA), Spanish (SP), and Vietnamese (VI). We consider AM
and BR accents as native and the rest foreign. Since our data set is highly
unbalanced across different accents, we assign different weights to each
subset. A weight of $n$ assigned to a subset means that audio clips from that
subset will appear around $n$ times during one epoch of training. The weights
can also be found in Table 1.

Table 1: Training Data Descriptions [13, 14, 15, 16, 17] Data Set |  Accents |  Duration(hrs) |  Speakers |  Prompts |  Weight
---|---|---|---|---|---
LibriTTS |  AM |  $585$ |  $2456$ |  LibriSpeech text |  1
VCTK |  BR |  $42$ |  $109$ |  Newspaper clippings |  6
SAA |  BR |  $3.7$ |  $579$ |  “Please call Stella…” |  10
L2-Arctic |  AR, HI, KO, MA, SP, VI |  $24$ |  $24$ ($2$ male and $2$ female per accent) |  ARCTIC prompts |  15
Indic TTS |  HI |  $20.06$ |  $2$ ($1$ male and $1$ female) |  ARCTIC prompts, Fairy tales |  2

###  3.2 Feature Processing

All audios are re-sampled at 16kHz and divided into $1.12$s segments for
training. We use the YAAPT algorithm to extract the F0 sequence at a frame
shift of $5ms$ and window length of $20ms$ [18]. We up-sampled the outputs of
the pronunciation encoder by repeating each time frame $4$ times since text
predictions from wav2vec2.0 are extracted at a $20ms$ frame shift.

###  3.3 Training Configurations

We trained our model with learning rate = $0.0002$, and we decayed our
learning rate by factor of $0.999$ every $1000$ iterations. We used the AdamW
optimizer with $\beta_{1}=0.8$ and $\beta_{2}=0.99$. For adversarial training,
we optimized the accent discriminator for the first $50000$ iterations before
applying the adversarial loss to the acoustic encoder. The model was trained
for a total of $3$ million iterations using batch size $16$.

###  3.4 Evaluation Set-up

We conducted mean opinion score (MOS) tests to evaluate the performance of the
proposed model on audio naturalness and speaker similarity and an XAB test to
evaluate the results of accent conversion. A different set of $100$ listeners
participated in each listening test, and we recruited an additional $24$
“experts” to evaluate accent conversion results. As in [9], all audios were
shuffled before being presented to the listeners.

We trained an ablation model (AB) by disabling the accent discriminator. We
used the same 20 prompt examples used by [9] for evaluation. In our case,
instead of using all examples from the Hindi accented female speaker labeled
p248 in the VCTK dataset, we replaced the first 10 examples with the
corresponding ones from a male British accented speaker, labeled p226, from
the VCTK dataset. We excluded audios from p248 and p226 from our training set
for both proposed and ablation models so that they are unseen speakers. Using
our proposed model (P) and the ablation model (AB), we converted the $20$
audio clips to American (AM), Hindi (HI), and Korean (KO) accents.

##  4 Results

###  4.1 Audio Quality

The listeners were asked to rate the audio quality for the original (O) and
the converted (AB, P) clips on a five-point scale (1-bad, 5-excellent). As
seen in Table 2, the proposed model maintains a comparable quality to the
original audio. We hypothesize it does better than the ablation model because
accent invariance in the acoustic encoder improves generalization.

###  4.2 Speaker Similarity

The raters first listen to a reference audio from each speaker with different
text content. For each audio clip, the listeners were asked to rate how close
did the voice sound to that of the reference audio’s on a five-point scale.
The listeners were instructed to disregard the accent and recording
conditions. Table 2 demonstrates that most listeners believed that the voice
in the audios converted using the models sound almost as similar to the
reference audio as original audios from the speaker.

Table 2: MOS Study Results (with $95\%$ confidence interval) Samples | Audio Quality | Speaker Similarity
---|---|---
Original | $3.87\pm 0.07$ | $4.28\pm 0.06$
Ablation | $3.27\pm 0.12$ | $3.89\pm 0.10$
Proposed | $3.62\pm 0.09$ | $4.05\pm 0.09$

###  4.3 Accent Conversion

The listeners first listened to some reference audio clips with the target
accent. Then, for each pair containing an original audio clip and its
converted audio clip, the listeners were asked to choose the one that has a
closer accent to the reference audio. As seen from the results in Fig. 3, for
every target accent, the proposed model performs better than the ablation
model. Significantly more listeners preferred the converted audios as sounding
more American than the original audios. However, the listeners seemed unsure
about whether the converted audios sounded more Korean or Hindi, likely
because it is difficult to identify an accent after listening to only a few
audio clips. The goal of conversion to other foreign accents is to help a non-
native listener understand someone with a different foreign accent or to
create a better entertainment value for the people knowledgeable of the
desired accent. Therefore, we recruited $24$ additional “expert” participants
who were born and raised in USA, India, and Korea, to judge the corresponding
accents. The results shown in Fig. 4 indicate that the listeners who are very
familiar with the target accent were confident that the proposed model
converted audios to sound like the target accent.

Fig. 3: Preference results from $100$ random listeners. Inconclusiveness on
accents other than American suggests that recognizing attributes of a specific
accent may need more exposure than obtained from listening to a couple
reference clips. Yet, converted samples from the proposed system are preferred
more than the ablation.

Fig. 4: Preference results from raters having substantial exposure to
respective accents (4 raters per accent). With the exception of HI-to-HI
conversion for p248, the listeners strongly preferred the audios converted by
our proposed model. Preference is more evident when converting from native to
non-native accents or vice versa.

##  5 Conclusion

In this paper, we have presented a novel accent conversion model that converts
an unseen speaker’s utterance with an arbitrary accent to utterances with
multiple different target accents222 Demo samples at https://accent-
conversion.github.io/oct-2022-demo/. The model is able to noticeably convert
the accent of input audios while preserving speaker identity and audio
quality. Future work may improve the perceptive accent accuracy of the
converted audio by modifying prosody, but the fact that the outputs of our
model are synchronized with the inputs makes the model useful for applications
such as dubbing a video with different accents in the original speaker’s
voice.

## References

  * [1] Shaojin Ding, Guanlong Zhao, and Ricardo Gutierrez-Osuna,  “Accentron: Foreign accent conversion to arbitrary non-native speakers using zero-shot learning,”  Computer Speech & Language, vol. 72, pp. 101302, 2022.
  * [2] Guanlong Zhao, Sinem Sonsaat, John Levis, Evgeny Chukharev-Hudilainen, and Ricardo Gutierrez-Osuna,  “Accent conversion using phonetic posteriorgrams,”  in 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2018, pp. 5314–5318.
  * [3] Guanlong Zhao and Ricardo Gutierrez-Osuna,  “Using phonetic posteriorgram based frame pairing for segmental accent conversion,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 27, no. 10, pp. 1649–1660, 2019.
  * [4] Sandesh Aryal and Ricardo Gutierrez-Osuna,  “Can voice conversion be used to reduce non-native accents?,”  in 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2014, pp. 7879–7883.
  * [5] Yu Zhang, Ron J. Weiss, Heiga Zen, Yonghui Wu, Zhifeng Chen, R. J. Skerry-Ryan, Ye Jia, Andrew Rosenberg, and Bhuvana Ramabhadran,  “Learning to speak fluently in a foreign language: Multilingual speech synthesis and cross-language voice cloning,”  Proceedings of the Annual Conference of the International Speech Communication Association, INTERSPEECH, vol. 2019-Septe, pp. 2080–2084, 2019\.
  * [6] Zhichao Wang, Wenshuo Ge, Xiong Wang, Shan Yang, Wendong Gan, Haitao Chen, Hai Li, Lei Xie, and Xiulin Li,  “Accent and speaker disentanglement in many-to-many voice conversion,”  2021 12th International Symposium on Chinese Spoken Language Processing, ISCSLP 2021, pp. 2–6, 2021.
  * [7] Andros Tjandra, Ruoming Pang, Yu Zhang, and Shigeki Karita,  “Unsupervised learning of disentangled speech content and style representation,”  Proceedings of the Annual Conference of the International Speech Communication Association, INTERSPEECH, vol. 4, pp. 3191–3195, 2021.
  * [8] Adam Polyak, Yossi Adi, Jade Copet, Eugene Kharitonov, Kushal Lakhotia, Wei Ning Hsu, Abdelrahman Mohamed, and Emmanuel Dupoux,  “Speech resynthesis from discrete disentangled self-supervised representations,”  Proceedings of the Annual Conference of the International Speech Communication Association, INTERSPEECH, vol. 5, pp. 3531–3535, 2021.
  * [9] Songxiang Liu, Disong Wang, Yuewen Cao, Lifa Sun, Xixin Wu, Shiyin Kang, Zhiyong Wu, Xunying Liu, Dan Su, Dong Yu, et al.,  “End-to-end accent conversion without using native utterances,”  in ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2020, pp. 6289–6293.
  * [10] Alexei Baevski, Henry Zhou, Abdelrahman Mohamed, and Michael Auli,  “wav2vec 2.0: A framework for self-supervised learning of speech representations,”  Advances in Neural Information Processing Systems, vol. 2020-Decem, pp. 1–12, 2020.
  * [11] Guanlong Zhao, Shaojin Ding, and Ricardo Gutierrez-Osuna,  “Foreign accent conversion by synthesizing speech from phonetic posteriorgrams.,”  in INTERSPEECH, 2019, pp. 2843–2847.
  * [12] Jungil Kong, Jaehyeon Kim, and Jaekyoung Bae,  “Hifi-gan: Generative adversarial networks for efficient and high fidelity speech synthesis,”  Advances in Neural Information Processing Systems, vol. 2020-Decem, 2020.
  * [13] Heiga Zen, Viet Dang, Rob Clark, Yu Zhang, Ron J. Weiss, Ye Jia, Zhifeng Chen, and Yonghui Wu,  “Libritts: A corpus derived from librispeech for text-to-speech,”  CoRR, vol. abs/1904.02882, 2019.
  * [14] Christophe Veaux, Junichi Yamagishi, Kirsten MacDonald, et al.,  “Cstr vctk corpus: English multi-speaker corpus for cstr voice cloning toolkit,”  University of Edinburgh. The Centre for Speech Technology Research (CSTR), 2017.
  * [15] Nishanthi N L Arun Baby, Anju Leela Thomas and Hema A Murthy,  “The speech accent archive: towards a typology of english accents,”  in Community-based Building of Language Resources (International Conference on Text Speech and Dialogue), pp. 37–43. 2016.
  * [16] Guanlong Zhao, Sinem Sonsaat, Alif Silpachai, Ivana Lucic, Evgeny Chukharev-Hudilainen, John Levis, and Ricardo Gutierrez-Osuna,  “L2-arctic: A non-native english speech corpus,”  in Proc. Interspeech, 2018, p. 2783–2787.
  * [17] Steven H A. Baby and Stephen A Kunath,  “Resources for indian languages,”  in Corpus-based studies in language use, language learning, and language documentation, pp. 265–281. Brill, 2011.
  * [18] Kavita Kasi and Stephen A Zahorian,  “Yet another algorithm for pitch tracking,”  in 2002 ieee international conference on acoustics, speech, and signal processing. IEEE, 2002, vol. 1, pp. I–361.


# Transformation Hydrodynamic Metamaterials: Rigorous Arguments on Form
Invariance and Structural Design with Spatial Variance

Gaole Dai<EMAIL_ADDRESS>School of Sciences, Nantong University, Nantong
226019, China  Jun Wang<EMAIL_ADDRESS>School of Physics, East China
University of Science and Technology, Shanghai 200237, China  Wenzhou
Institute, University of Chinese Academy of Sciences, Wenzhou 325001, China

###### Abstract

The method of transformation optics has been a powerful tool to manipulate
physical fields if governing equations are formally invariant under coordinate
transformations. However, regulation of hydrodynamics is still far from
satisfactory due to the lack of rigorous arguments on the validation of
transformation theory for various categories of fluids. In this paper, we
systematically investigate the applicability of transformation optics to fluid
mechanics. We find that the Stokes equation and the Navier-Stokes equations,
respectively describing the Stokes flow and general flow, will alter their
forms under curvilinear transformations. On the contrary, the Hele-Shaw flow
characterized with shallow geometries rigidly retain the form of its governing
equation under arbitrary transformations. Based on the derived transformation
rules, we propose the design of multilayered structures with spatially varying
cell depth, instead of engineering the rank-2 shear viscosity tensor, to
realize the required anisotropy of transformation Hele-Shaw hydrodynamic
metamaterials. The theoretical certify and fabrication method revealed in this
work may pave an avenue for precisely controlling flow distribution with the
concept of artificial structure design.

##  I Introduction

Transformation optics (TO), based on the form invariance of Maxwell’s
equations under coordinate transformations, has triggered the discovery of
many exotic phenomena and rapid progress of metamaterials in electromagnetism
Science06-01 ; Science06-02 ; siam ; NM2010 ; laser . This geometric approach
to mimic physics in curved space has also been successfully applied to other
wave or diffusion systems rpp ; AM2021 . It is notable that not all physical
phenomena satisfy the form invariance required by TO. For example, Milton et
al. found that the form of the conventional elastodynamic equation with a
rank-4 elasticity tensor actually changes unless certain modifications are
added NJP2006 . However, such a rigorous treatment on the validation of TO is
still absent in another fundamental topic of continuum theories, the fluid
mechanics, whose governing equations include the Navier-Stokes (NS) equations
known for its complexity landau-6 . Early works on transforming fluid
mechanics took different ways to simplify the NS equations and realized
functions like invisibility in underwater acoustic waves NJP2007 ; APL2007 ;
PRB2017 , liquid surface waves epl09 ; prl2022hy ; prl08  and Darcy flows in
porous media PRL2011 ; PRE2018 . Recently, a series of transformation
hydrodynamic metamaterials PRL2019 ; EML2020-2 ; PRAP2019  were experimentally
demonstrated in Stokes or creeping flows with anisotropic inhomogeneous
viscosities, built upon the form invariance of the Stokes equation (neglecting
the inertial term in the NS equations). Later, this approach was extended to
general laminar flows EML2020-1 , which assumes the NS equations are also
form-invariant. Although the recent works appear to greatly expand the
applicability of TO for fluid manipulation, we find that there are still some
problems with the strictness of the theory that cannot be ignored.

First, the general viscosity should be a rank-4 tensor, while the recent works
PRL2019 ; EML2020-2 ; PRAP2019 ; EML2020-1  rely on the assumption of a rank-2
tensor as an anisotropic version of scalar shear viscosity. We have not yet
checked whether it can be reconciled to realistic viscous constitutive
relation. Second, even if the viscosity can be reduced to a rank-2 tensor, the
Stokes flows still might not meet the requirements of TO, let alone general
(laminar) flows pd . In fact, we prove that neither the Stokes equation nor
the NS equations can satisfy the form invariance in this work with both
theoretical and numerical demonstrations. However, recent works PRL2019 ;
EML2020-2 ; PRAP2019 ; EML2020-1  do give valid evidences that hydrodynamic
metamaterials with an effective rank-2 viscosity tensor could work under
certain conditions, both numerically and experimentally. Therefore, we also
aim to find the conditions under which the TO theory can hold in fluid
mechanics. In fact, some studies arxiv19 ; PR2021 ; arxiv20 ; PRL2021  have
noticed that Refs. PRL2019 ; PRAP2019 ; EML2020-2 ; EML2020-1  actually used
shallow geometries, i.e., the Hele-Shaw cells nature1898 , which restrict the
flows inside two plates (usually parallel) with a small gap. Hele-Shaw cells
or Hele-Shaw flows are widely used in many fields like microfluidics ARFM2004
; AC2002 . We further prove that the Stokes equation in Hele-Shaw cells can be
reduced to a form that satisfies TO when considering anisotropy and
inhomogenity. In addition, to realize hydrodynamic metamaterials in Hele-Shaw
cells, we propose an approach of transforming the depth of the cells instead
of the viscosities. We achieve anisotropic depths required for transformation
devices with engineered cell configurations in which layers with different
depths are alternately arranged.

##  II Necessity of a shallow geometry for TO

###  II.1 Theoretical modeling

In TO only involving spatial coordinate transformations, we consider a map
from the virtual space to the physical space. Usually, the virtual space is
isotropic and homogeneous. TO requires the equation’s forms in the two spaces
are the same or the form in the physical space can be an appropriate
generalization for anisotropic inhomogeneous fluids. In this work, we consider
steady incompressible isothermal Newtonian fluid. The governing equations in
the virtual space include the law of continuity

$\nabla\cdot\left(\rho\mathbf{v}(\mathbf{r})\right)=0,$ (1)

and the NS equations

$\rho(\mathbf{v}\cdot\nabla)\mathbf{v}=-\nabla p+\nabla\cdot{{\bm{\tau}}},$
(2)

where $\rho$ accounts for the density, $\mathbf{v}$ is the velocity,
$\mathbf{r}$ is the position vector, $p$ is the pressure, and ${{\bm{\tau}}}$
denotes the viscous stress tensor which can be expressed as

${{\bm{\tau}}}=\mu\left(\nabla\mathbf{v}+\nabla\mathbf{v}^{\top}\right).$ (3)

Here $\mu$ is the shear viscosity. For Stokes flows, the inertial term
$\rho(\mathbf{v}\cdot\nabla)\mathbf{v}$ can be neglected, and Eq. (2) can be
reduced to

$-\nabla p+\nabla\cdot{{\bm{\tau}}}=0.$ (4)

It is notable that the viscosity is a scalar in Eq. (3), which corresponds to
isotopic incompressible fluids with some common symmetries landau-6 . It is
easy to check that Eq. (1) is form-invariant, and its form in the physical
space is

$\nabla^{\prime}\cdot\left(\rho^{\prime}\mathbf{v}^{\prime}(\mathbf{r}^{\prime})\right)=0,$
(5)

where all the symbols are superscripted to represent quantities in the
physical space. Let J denote the Jacobian of the map
$\mathbf{r}\mapsto\mathbf{r}^{\prime}$. The transformation rule to make Eq.
(5) valid is PR2021

$\mathbf{v}^{\prime}(\mathbf{r}^{\prime})=\text{J}(\mathbf{r}^{\prime})\mathbf{v}(\mathbf{r}(\mathbf{r}^{\prime}))\operatorname{det}\text{J}^{-1}(\mathbf{r}^{\prime})$
(6)

since we can take density as a constant in incompressible isothermal flows.

Now, the question is whether Eqs. (2) and (4) are also form-invariant as
declared in Refs. PRL2019 ; PRAP2019 ; EML2020-2 ; EML2020-1 . Their
conclusions are based on a rank-2 viscosity tensor

$\bm{\mu}^{\prime}=\text{J}^{-\top}\text{J}^{-1}\mu\det\text{J}$ (7)

in the physical space, and the corresponding stress tensor is PRL2019 ;
PRAP2019 ; EML2020-2 ; EML2020-1

$\nabla^{\prime}\cdot\left({\bm{\mu}^{\prime}}\cdot\left(\nabla^{\prime}\mathbf{v}^{\prime}+(\nabla^{\prime}\mathbf{v}^{\prime})^{\top}\right)\right).$
(8)

However, the general form of the viscosity is a rank-4 tensor denoted by
$\mathbb{A}$ or $A_{ijkl}$ and the stress tensor should also be more
complicated than Eq. (8) landau-6 . When general anisotropy is introduced,
$\mathbb{A}$ cannot be reduced to a rank-2 tensor; see our discussion in Part
A of Appendix. Even if such a rank-2 viscosity tensor does exist
phenomenologically, whether Eqs. (3) and (8) can be transformed into each
other should also be reconsidered. Anyway, the experimental results in
Refs.PRL2019 ; PRAP2019 ; EML2020-2  do give evidences that TO can work in
some special flows, so the difference between the observed systems and general
Stokes flows or laminar flows should not be overlooked, i.e., the shallow
geometries in experiments.

Figure 1: (a) Geometry of a 3D cell with a depth denoted as $h$. The top view
of the cell is a rectangle with a length $L$ and a width $W$. The two $y$-$z$
boundaries (with regard to $x=\pm L/2$) are the flow inlet and outlet,
respectively. Other four boundaries (with regard to $y=\pm L/2$ and $z=\pm
h/2$) are nonslip walls. (b) Geometry of a 2D Hele-Shaw cell in the $x$-$y$
plane. It also represents the virtual space for coordinate transformations to
design hydrodynamic metamaterials. The region outside the red dashed circle
(with a radius $R_{2}$) undergoes the identity transformation. The region
inside the circle can be transformed into devices with different manipulation
functions including: (c) cloak; (d) rotator; (e) concentrator. The
transformation effects are also shown by the curved meshes.

To see the effect of shallow geometries, we consider flows in a cell or a
cuboid-shaped pipe [Fig. 1(a)]. Its length and width are $L$ and $W$,
respectively. For simplicity, we take $L=W$ in the following discussion. $h$
denotes its height or the fluid depth if the cell is fully occupied. The inlet
and outlet are on different open borders parallel to the $y$-$z$ plane. The
other four sides are solid walls with nonslip conditions. Below we focus on
the Stokes flows in the cell. If the Stokes equation changes its form under
coordinate transformations, so do the NS equations unless the inertial term
can cancel out the form change of the stress tensor. A shallow geometry means
$h\ll L$ (e.g., $h/L=5\times 10^{-3}$ in Ref. PRL2019 ) and then the
lubrication approximation can be used to obtain the velocity distribution for
isotropic fluids:
$\mathbf{v}=-\frac{h^{2}}{8\mu}\left(1-\left(\frac{2z}{h}\right)^{2}\right)\nabla
p$ InFlow . The average velocity in the $x$-$y$ plane (integrating along the
$z$ axis with a symmetry to the $z=0$ plane) is InFlow

$\overline{\mathbf{v}}=\frac{2}{h}\int_{0}^{\frac{h}{2}}\mathbf{v}dz=-\frac{h^{2}}{12\mu}\nabla
p.$ (9)

If we take the average velocity as an approximation of $\mathbf{v}$, the
three-dimensional (3D) model can be viewed as a two-dimensional (2D) one,
i.e., the Hele-Shaw flow. Considering the law of continuity [Eq.(1)] at the
same time, the governing equation of Hele-Shaw flows (Hele-Shaw equation) can
be written as InFlow

$\nabla\cdot\left(\frac{h^{2}}{12\mu}\nabla p(\mathbf{r})\right)=0,$ (10)

which reduces to the Laplace equation if $\frac{h^{2}}{12\mu}$ is a constant.
In addition, the Hele-Shaw flow has the same form as the Darcy’ law since
$\frac{h^{2}}{12\mu}$ behaves as the effective permeability. Similar to the
work on transforming Darcy’s law PRL2011 , Eq. (10) is also mathematically
form-invariant under coordinate transformations (The detailed derivation is
presented in Part B of Appendix). Its form in the physical space is

$\nabla^{\prime}\cdot\left(\text{J}\text{J}^{\top}\det\text{J}^{-1}\frac{h^{2}}{12\mu}\nabla^{\prime}p^{\prime}(\mathbf{r}^{\prime}(\mathbf{r}))\right)=0.$
(11)

If the depth keeps unchanged, we can obtain the same transformation rule for
the viscosity as Eq. (7). Further, we also need to verify whether the 2D model
described by Eq. (11) can be derived from anisotropic flows with a rank-4
viscosity tensor in 3D shallow cells. Based on the lubrication approximation,
we can give the affirmative conclusion. Similarly, the form of Eq. (9) in the
physical space is
$\mathbf{v}^{\prime}(\mathbf{r}^{\prime})=-\text{J}\text{J}^{\top}\det\text{J}^{-1}\frac{h^{2}}{12\mu}\nabla^{\prime}p^{\prime},$
and Eq. (6) still holds. The detailed derivation on form invariance is offered
in Part C of Appendix. On the other hand, in Part D of Appendix, we give a
rigorous analysis on how the form of the Stokes equation and the NS equations,
especially the viscous stress term in them, change under curvilinear
transformations.

A more intuitive way is to use numerical simulations to check the performance
of transformation devices under different flow models. For that, we consider a
modification to Eqs. (10) and (11). If we assume the Stokes equation is form
invariant with a rank-2 viscosity tensor, the combination of Eqs. (4) and (9),
i.e.,

$\nabla
p-\nabla\cdot\left(\mu\left(\nabla\mathbf{v}+\nabla\mathbf{v}^{\top}\right)\right)+\frac{12\mu}{h^{2}}\mathbf{v}=0,$
(12)

should also be form invariant since their transformation rules are all the
same. When the depth is small ($h/L\to 0$), Eq. (12) reduces to Eq. (9). In
contrast, if $h$ is large enough, Eq. (12) will reduce to Eq. (4). The
anisotropic version of Eq. (12) is

$\nabla^{\prime}p^{\prime}-\nabla^{\prime}\cdot\left({\bm{\mu}^{\prime}}\cdot\left(\nabla^{\prime}\mathbf{v}^{\prime}+(\nabla^{\prime}\mathbf{v}^{\prime})^{\top}\right)\right)+\frac{12}{h^{2}}{\bm{\mu}^{\prime}}\cdot\mathbf{v}^{\prime}=0.$
(13)

Therefore, we can deal with Hele-Shaw flows and more general Stokes flows in a
unified 2D model by considering Eqs. (12) and (13) with different values of
$h/L$. This also helps alleviate the huge demands on computing resources for
3D simulations. If we can construct an example that changes the form of Eq.
(12) under coordinate transformations when $h/L$ is not very small, we can
prove that TO is not valid in general Stokes flows.

###  II.2 Numerical verification

Numerical simulations are performed to verify our assertions about the form
invariance of fluid mechanisms equations. We use the Coefficient Form PDE
module in the commercial finite-element software COMSOL Multiphysics
(https://www.comsol.com/) to establish the anisotropic Hele-Shaw model. The 2D
model has the same geometry and boundary conditions as Fig. 1(b). In
simulations, we set $L=10$ mm and apply a pressure bias $\Delta p$ on the $x$
direction. The properties of the background material are referenced to water,
including the viscosity $\mu=10^{-3}$ Pa s and the density $\rho=1000$ kg m-3.
We use three commonly used transformation devices, i.e., the cloak [Fig.
1(c)], the rotator [Fig. 1(d)], and the concentrator [Fig. 1(e)].

Figure 2: Simulation results for cells with different depths. (a1)–(a3) show
the pressure distributions of the cloak when $h/L$ ($\Delta p$) is $10^{-3}$
(10 Pa), $10^{-2}$ ($10^{-1}$ Pa) and $10^{-1}$ ($10^{-3}$ Pa), respectively.
The white lines are isobars. (b1)–(b3)/(c1)–(c3)/(d1)–(d3) are the
corresponding pressure distributions of the reference/rotator/concentrator.

First, we verify the performance of hydrodynamic cloaks. Let $(r,\theta)$ and
$(r^{\prime},\theta^{\prime})$ denote the radius and azimuth of the polar
coordinates in the virtual and physical spaces, respectively. The
corresponding transformation Science06-01  is
$r^{\prime}=R_{1}+\frac{R_{2}-R_{1}}{R_{2}}r$ $(0<r<R_{1})$ [Figs. 1(b) and
1(c)]. The cloak is the shell $S=\\{\mathbf{r}^{\prime}:R_{1}\leqslant
r^{\prime}\leqslant R_{2}\\}$ which prevents the obstacle in it from
disturbing the outside flows. Figs. 2(a1)–2(a3) present the pressure
distributions when the cloak is put in cells with different depths. The region
inside the cloak is a solid obstacle so it can be removed from the simulated
domain with nonslip condition on $\\{\mathbf{r}^{\prime}:r^{\prime}=R_{1}\\}$.
For comparison, the corresponding reference pressure distributions in the
virtual space are also illustrated in Figs. 2(b1)–2(b3). Here, we take
$R_{1}=2$ mm and $R_{2}=3$ mm, and let $\Delta p$ vary with $h$ or $h/L$ so
that the Reynolds numbers (Re) can be close in all the cases. In fact, we have
$\text{Re}\approx 8.3\times 10^{-5}$ for the reference in the following
simulations if we use $h$ as the characteristic length. This small value of Re
ensures that the flow is still a creeping one even if Re is increased in cases
with a cloak. The pressure patterns in Figs. 2(b1)–2(b3) are visually
indistinguishable while those in Figs. 2(a1)–2(a3) are quite different. The
cloak with regard to $h/L=10^{-3}$ [Fig. 2(a1)] has a background pressure
pattern looking no different from its counterparts in Fig. 2(b1). In Fig.
2(a2) with a larger depth, we can see the isobars near the cloak are slightly
curved. When $h/L$ increases to $10^{-1}$ in Fig. 2(a3), the isobars are
sharply bent, and the obstacle can be easily detected. We can see the shallow
geometry is necessary to maintain invisibility. When the stress tensor term is
neglected, i.e., using the real Hele-Shaw equation, we can obtain a nearly
perfect cloak for any value of $h/L$ in simulations, although the
corresponding 3D models may no longer be called shallow. In this cases, all
the pressure distributions should have the same pattern close to Fig. 2(a1).
Here, we say nearly perfect because the nonslip condition requires zero
velocity at $\\{\mathbf{r}^{\prime}:r^{\prime}=R_{1}\\}$. This conflicts with
Eq. (6), which does not prohibit tangential flow around the obstacle. In Fig.
2(a3), both the changes of equation forms and the nonslip boundary contribute
to the failure of cloaking.

Further, we check the performance of another two invisible devices without
introducing extra nonslip boundaries. Figs. 2(c1)–2(c3) present the simulated
pressure distributions of the (invisible) rotator. Its transformation [Figs.
1(b) and 1(d)] is $\theta^{\prime}=\theta+\theta_{0}$ $(r<R_{1})$ and
$\theta^{{}^{\prime}}=ar+b$ $(R_{1}\leqslant r\leqslant R_{2})$, where
$a=\frac{\theta_{0}}{R_{1}-R_{2}}$,
$b=\theta+\frac{R_{2}}{R_{2}-R_{1}}\theta_{0}$ ChenAPL07 . The rotator itself
is still the shell $S$ and the physical fields inside it are expected to be
rotated by $\theta_{0}$ (counterclockwise for a positive value), which is
$-60\degree$ in Figs. 2(c1)–2(c3). Also, the region outside it undergoes an
identity transformation so the pressure distribution should keep invariant
(invisibility). Similar to the cloak, the first two rotators in Figs. 2(c1)
and 2(c2) with shallower geometries work better than that in Fig. 2(c3).
Another typical device in TO is the (invisible) concentrator in which the
physical fields can be amplified without changing the outside flows. Its
transformation [Figs. 1(b) and 1(e)] is $r^{\prime}=\frac{R_{1}}{R_{3}}r$
$(r<R_{3})$ and
$r^{\prime}=\frac{R_{1}-R_{3}}{R_{2}-R_{3}}R_{2}+\frac{R_{2}-R_{1}}{R_{2}-R_{3}}r$
$(R_{3}\leqslant r\leqslant R_{2})$, where $R_{3}$ is an auxiliary parameter
that determines the concentrating effect RahmPNFA08 . Here we take $R_{3}=2.5$
mm. The corresponding simulation results are shown in Figs. 2(d1)–2(d3),
demonstrating denser isobars in the center. It is interesting that when $h/L$
reaches $10^{-1}$, the uniform pressure distribution outside the concentrator
hasn’t been distorted much. In fact, it needs a quite large $h/L$ to see the
obvious failure of invisible concentrator. Anyway, we can conclude that Eq.
(12) cannot be transformed to Eq. (13). At the same time, we can confirm the
form invariance of the Hele-Shaw flows. Thus, the only reasonable explanation
for Fig. 2 is that the Stokes flows are not form invariant under coordinate
transformations. If we add the inertial term back, the pressure distributions
will be almost the same as those in Fig. 2 due to the creeping nature. So the
general laminar flows governed by the NS equations are not form invariant
either.

Figure 3: Simulated velocity in cells with different depths. (a1)–(a3)
Pressure distributions of the cloak when $h/L$ is $10^{-3}$, $10^{-2}$ and
$10^{-1}$, respectively. Different colors represent the magnitude of velocity.
The black lines with arrows are streamlines. (b1)–(b3)/(c1)–(c3)/(d1)–(d3) The
corresponding velocity distributions of the reference/rotator/concentrator.
(e1)–(e3) Speed data read from the line $x=0$ for the devices or
$x^{\prime}=0$ for the reference from (a1)–(d3): (e1) $h/L=10^{-3}$; (e2)
$h/L=10^{-2}$; (e3) $h/L=10^{-1}$. The different shades behind the plotted
lines represent the background area (white), the device (lighter gray), and
the inside of the device (darker gray).

Velocity distribution [Fig. 3] provides another perspective for judging and
analyzing the validation of TO in different cells. By comparing the background
velocity (invisibility) [Figs. 3(a1)–3(d3)] and checking the uniformness
inside the rotators [Figs. 3(c1)–3(c3)] and concentrators [Figs. 3(d1)–3(d3)],
we can obtain the same conclusion as that from comparing pressure in Fig. 2.
In particular, Fig. 3(d3) shows the deviation of the concentrator from the
ideal effect more clearly than Fig. 2(d3). Also, we can see the viscous drag
plays an important role in the failure of invisible devices. Besides the inner
edge of the cloak, non-slip also exists on the side walls of the cell [Fig.
1(b)]. In Figs. 3(b1)–3(b3), the boundary layer becomes thicker as $h/L$
increases in the virtual space. The flows in Figs. 3(b1) and 3(b2) are
approximately uniform while that in Fig. 3(b3) has an obvious gradient along
the direction of $y$-axis near the nonslip boundaries. The latter case
contradicts the lubrication approximation, which requires the velocity
gradient in the $x$-$y$ plane to be negligible compared to that along the
direction of $z$-axis. Such boundary layer features are also present in the
three devices composed of hydrodynamic metamaterials. More specifically, we
plot the corresponding speed on the line $x=0$ for the reference or
$x^{\prime}=0$ for the devices in Figs. 3(e1)–3(e3), which intuitively show
the size of the boundary layers.

Based on our analysis in Part D of Appendix, some constraints that can make
the form of Stokes flows transformation-invariant include (1) an isotropic
homogeneous virtual space that is drag-free and (2) a pressure distribution
satisfying Eq. (A36) (differing from Eq. (11) by a coefficient before the
pressure gradient) in the physical space at the same time. It is easy to see
that both of these conditions no longer hold when the lubrication
approximation fails. In addition, the relatively better performance of the
concentrator in deep cells can be explained by the less space distortion
caused by its transformation. In Part E of Appendix, we do another set of
simulations with different parameters for the three devices and quantify their
invisibility effects together with the results in Figs. 2 and 3. By adjusting
$R_{3}$ for a larger concentrating effect, we observe a more obvious failure
of the concentrator at $h/L=10^{-1}$. In conclusion, the rigorous
establishment of TO needs lubrication approximation although some devices can
still have a well performance in non-shallow flows if the geometric
transformations only cause small space distortions.

##  III Multilayered structure with spatial variance

The devices designed by TO need anisotropic inhomogeneous viscosities, which
are quite difficult to fabricate. A common technique to simplify the
fabrication in metamaterials is discretizing the space and then achieving a
specific isotropic (or diagonally anisotropic) and uniform material property
in each unit oe07 ; PRB2008 ; PRL2012con . For example, the cloak can be
realized by multilayered structures composed of two bulk materials: one has a
higher and the other has a lower viscosity than the background. Further, it is
still not easy to flexibly increase or decrease the scalar shear viscosity.
External factors such as magnetic field pre2011  and temperature prap2021sb
can change the viscosity in a wide range, but this will make the system more
complicated. Previous studies PRL2019 ; EML2020-2 ; PRAP2019  put pillars into
the fluid, which is an effective method to increase the viscosity. However,
this technique can’t reduce the viscosity so pillars were also put into the
background region to match the parameters. On the other hand, putting in
pillars with the same height as the cell actually reduces the fluid depth to
zero. So, we can use a cell with a spatially varying depth to avoid additional
adjustments to the background flow. This also means that we can go back to the
realistic 3D model. Actually, in metamaterials controlling water waves,
engineering depth has been an important technique hs ; PRL2021zj  and even
anisotropic depth can be realized epl09 ; prl2022hy . Due to the variable flow
cross section, the law of continuity in the physical space needs to be
modified as PF2013 ; pre2013 ; JFM2014

$\nabla^{\prime}\cdot(h^{\prime}(x^{\prime},y^{\prime})\mathbf{v^{\prime}})=0,$
(14)

and then Eq. (11) should be generalized to

$\nabla^{\prime}\cdot\left(\frac{(h^{\prime})^{3}}{12\mu}\nabla^{\prime}p^{\prime}\right)=0.$
(15)

It’s easy to check Eqs. (14) and (15) are mathematically form-invariant,
although how the transformation matrix
$\text{J}\text{J}^{\top}\det\text{J}^{-1}$ acts on $h$ is a bit tricky.
Nevertheless, we can still try to use multilayered structures to approximately
mimic anisotropy.

Figure 4: Structures of multilayered metamaterials with variable depth and
their simulation results. (a1)–(a3) Top view (or 2D model) of the devices. The
regions illustrated by different colors have different depth $h^{\prime}$ in
the $z^{\prime}$ direction. (a1) The cloak has 10 annular layers with the same
difference between inner and outer radii. (a2) The rotator has 40 layers whose
helical boundaries are determined by $x^{\prime}=R_{1}\exp(s)\cos(ks+n\pi/20)$
and $y^{\prime}=R_{1}\exp(s)\sin(ks+n\pi/20)$, where
$s\in\left[0,\,\ln\left(R_{2}/R_{1}\right)\right]$,
$k=-\theta_{0}/(\ln(R_{2}/R_{1}))$ and $n$ is an integer from 1 to 40 PRB2008
. (a3) The concentrator has 40 layers with symmetry of a $9\degree$ rotation.
(b1)–(b3) Side views of the 3D multilayered devices. They all have a symmetry
to the $z^{\prime}=0$ plane, which can be generated by extending the 2D model
(plane $z^{\prime}=0$) by $h^{\prime}/2$ in both the positive and negative
directions of the $z^{\prime}$-axis. (c1)–(c3) are the pressure distributions
on the surface (top view, i.e., the upper boundary $z^{\prime}=h^{\prime}/2$).
(d1)–(d3) are the corresponding pressure distributions on the central plane
$z^{\prime}=0$. (e1)–(e3) illustrate the horizontal components of the velocity
vector on $z^{\prime}=0$. The red arrows show the direction of horizontal
velocity. Their length and the colors show the magnitude of horizontal
velocity. We set $L=D=10$ mm, $h/L=2\times 10^{-3}$ and $a=2$, and apply a
pressure bias equal to 1 Pa along the $x^{\prime}$-axis.

Fig. 4 shows the schematic diagram of the multilayered structures. For the
cloak, the annular shell in the plane $z^{\prime}=0$ is divided into ten
concentric annuli [Fig. 4(a1)] and the 3D structure is generated by extruding
the plane in both directions of the $z^{\prime}$-axis by
$h^{\prime}(x^{\prime},y^{\prime})/2$ [Fig. 4(a2)]. In each layer (annulus),
$h^{\prime}$ is a constant. The product of $h^{\prime}$ of two adjacent layers
is equal to $h^{2}$. Similarly, the 3D rotator [Figs. 4(a2) and 4(b2)] or
concentrator [Figs. 4(a3) and 4(b3)] can also be designed by dividing the
shell into helical or rotationally symmetrical multilayered structures.
Figure. 4 also confirms our design by 3D numerical modeling. For simplicity,
we show the (flattened) pressure distributions on the upper surface of the
cell in Figs. 4(c1)–4(c3) and those on the central plane $z^{\prime}=0$ in
Figs. 4(d1)–4(d3). The difference of the pressure distributions on the surface
and the central horizontal section of the same device is actually very small.
Figures. 4(e1)–4(e3) illustrate the horizontal velocity distributions on plane
$z^{\prime}=0$. We can see the demonstrated pressure and velocity
distributions indicate the good invisibility performance of multilayered
structures, although the isobars for the cloak in Figs. 4(c1) and 4(d1) are
slightly bent–due to the extra boundary layer and approximation of the 3D
model itself. From Fig. 4(e2), the velocity is turned about $38\degree$
clockwise inside the rotator, not fully reaching the theoretical value
$60\degree$ clockwise. For the concentrator [Fig. 4(e3)], the velocity inside
it is amplified by a factor of 1.31 (calculated based on the spatial average),
exceeding $R_{m}/R_{1}=1.25$. In fact, the expected concentrating effect
should be $R_{2}/R_{1}$ for multilayered structures since $R_{m}$ doesn’t
affect the design PRL2012con . The performance of our designs can be enhanced
by using more layers, larger depth difference between two adjacent layers,
and, of course, a shallower cell.

##  IV Conclusion

In summary, we give a careful discussion on choosing a proper flow model to
design hydrodynamic transformation metamaterials. Based on rigorous
mathematical proof and numerical verification, we show that the form
invariance required by TO can hold in Hele-Shaw flows, but not in general
Stokes flows or laminar flows. When the viscous drag cannot be neglected, the
flows in cells will deviate from the lubrication approximation. Then, the
obstacles inside the cloaks become detectable, and the rotators and
concentrators can no longer achieve the desired regulation effects.
Furthermore, we propose a method to equivalently achieve the desired
anisotropic viscosity of metamaterials by spatially varying cell depth. Our
results not only clarify the formal basis for the application of TO in fluid
mechanics, but also reveal a freedom of manipulating fluids in cells through
spatially-varying configuration.

## Acknowledgments

We acknowledge financial support by the National Natural Science Foundation of
China under Grants No. 12147169 and 12205101.

## Appendix

###  IV.1 Can the rank-4 viscosity tensor be reduced to a rank-2 one?

Here we give a brief discussion on the existence of a rank-2 viscosity tensor.
In general, the viscosity is a rank-4 tensor $\mathbb{A}$ and the constitutive
relation of stress tensor is

${\bm{\tau}}=\mathbb{A}:\nabla\mathbf{v}=A^{ijkl}\epsilon_{kl}\mathbf{g}_{i}\otimes\mathbf{g}_{j}.$
(A1)

Since $\nabla\mathbf{v}$ is a rank-2 tensor, we can remark it as
$\bm{\epsilon}=\epsilon^{ij}\mathbf{g}_{i}\otimes\mathbf{g}_{j}$. Although a
phenomenological rank-2 tensor or a matrix of shear viscosity has been
revealed in surface water waves prl08  and the recent works on hydrodynamic
metamaterials PRL2019 ; EML2020-2 ; PRAP2019 ; EML2020-1 , it’s still worth
verifying the reconciliation between the two forms of viscosity. For
simplicity, we use the Cartesian coordinate system in this note and Eq. (A1)
can be expressed as

$\tau_{ij}=A_{ijkl}\frac{\partial v_{k}}{\partial x_{l}}.$ (A2)

For 3D flows, $A_{ijkl}$ has 81 components and the independent ones can be
reduced under some common assumptions. For example, with the Onsager
reciprocal relations ($A_{ijkl}=A_{klij}$), the rotational symmetry
($A_{ijkl}=A_{jikl}$) and no-internal-friction condition under a pure uniform
rotation ($A_{ijkl}=A_{ijlk}$) landau-6 , only 21 independent components are
left and the number is reduced from 16 to 6 for 2D flows. More details on the
symmetry of viscosity can be found in Ref. arxiv21 . It’s obvious that 21 and
6 exceed the component numbers of a rank-2 tensor or matrix in 3D and 2D
spaces, respectively.

Since the rank-2 viscosity tensors in Refs. PRL2019 ; EML2020-2 ; PRAP2019 ;
EML2020-1  are only anisotropic in the horizontal plane, we consider the 2D
anisotropic case in the $x$–$y$ plane. With the three symmetries mentioned
above, the viscous stress tensor can be expressed as OM2003

$\begin{bmatrix}\tau_{xx}\\\ \tau_{yy}\\\
\tau_{xy}\end{bmatrix}=\begin{bmatrix}A_{1}&A_{2}&A_{3}\\\
A_{2}&A_{4}&A_{5}\\\
A_{3}&A_{5}&A_{6}\end{bmatrix}\begin{bmatrix}\dot{e}_{xx}\\\ \dot{e}_{yy}\\\
\dot{e}_{xy}\end{bmatrix}.$ (A3)

$\\{A_{1},A_{2},...,A_{6}\\}$ are the 6 independent components of $A_{ijkl}$.
Further, we only consider incompressible fluids, which gives an extra
restriction for the components in strain rate: $\dot{e}_{xx}+\dot{e}_{yy}=0$.
However, the number of independent components in $A_{ijkl}$ can only be
reduced to 5, not 4 or less, so

$\displaystyle\begin{bmatrix}\tau_{xx}&\tau_{xy}\\\
\tau_{yx}&\tau_{yy}\end{bmatrix}\neq\begin{bmatrix}\mu_{xx}&\mu_{xy}\\\
\mu_{yx}&\mu_{yy}\end{bmatrix}\begin{bmatrix}\dot{e}_{xx}&\dot{e}_{xy}\\\
\dot{e}_{xy}&\dot{e}_{yy}\end{bmatrix}$ (A4)

for any matrix $\begin{bmatrix}\mu_{xy}\end{bmatrix}$ as real valued functions
of $(x,y)$. In conclusion, for general anisotropic fluids, we can’t reduce the
rank-4 viscosity tensor to a rank-2 one.

###  IV.2 Form invariance of Hele-Shaw flows

In this note, we give a detailed proof for the form invariance of Hele-Shaw
flows. As the starting point, we consider a map $f$ in 3D Euclidean space
$\mathsf{E}^{3}$: $U\to V$ ($U,V\subseteq\mathsf{E}^{3}$) or
$\mathbf{r}\mapsto\mathbf{r}^{\prime}$ ($\mathbf{r}\in
U,\mathbf{r}^{\prime}\in V$) [Fig. A1]. $U$ (the virtual space) and $V$ (the
physical space) can seen as two sub-manifolds of the simple manifold
$\mathsf{E}^{3}$. The position vectors can be expressed by the global
Cartesian coordinates in $\mathsf{E}^{3}$, i.e.,
$\mathbf{r}=x\mathbf{e}_{x}+y\mathbf{e}_{y}+z\mathbf{e}_{z}$ and
$\mathbf{r}^{\prime}=f(\mathbf{r})=x^{\prime}\mathbf{e}_{x}+y^{\prime}\mathbf{e}_{y}+z^{\prime}\mathbf{e}_{z}.$
Here $\\{\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z}\\}$ is the orthonormal
basis. Then we can see $f$ has the same effect as another map $\hat{f}$
between two subsets $X,Y$ in real 3-space $\mathsf{R}^{3}$:

${\displaystyle{\begin{cases}\hat{f}:X\subseteq\mathsf{R}^{3}\to
Y\subseteq\mathsf{R}^{3},\\\
(x,y,z)\mapsto(x^{\prime},y^{\prime},z^{\prime}).\end{cases}}}$ (A5)

Here $X=\chi(U)$, $Y=\psi_{\alpha}(V)$ and
$\hat{f}=\psi_{\alpha}\circ(f\circ\chi^{-1})$. $\chi$ and $\psi_{\alpha}$ are
both bijections (or more exactly, homeomorphisms) that use only one coordinate
map to give the coordinate of every point in each manifold, i.e., the atlas
jost . The pairs of manifold and coordinate map, i.e., (coordinate) charts,
can be denoted as $(U,\chi)$ [Figs. A1(a) and A1(c)] and $(V,\psi_{\alpha})$
[Figs. A1(f) and A1(d)] . For TO, $f$ must be a bijection and thus its inverse
$f^{-1}$ exists. Further, both $f$ and $f^{-1}$ should be smooth enough so $f$
and $\hat{f}$ are diffeomorphisms. In this way, Jacobian
$\text{J}_{x^{\prime}x}=\partial x/\partial x^{\prime}$ and its inverse are
both well-defined.

Figure A1: Illustration of the ways to understand the form invariance of TO.
(a) and (b) both represent manifold $U$. The different meshes corresponding to
different coordinate maps. (c) and (d) are value ranges in $\mathsf{R}^{3}$,
linked by the map $\hat{f}$. (e) and (f) both represent manifold $V$. $f$ is
the map (diffeomorphism) from $U$ to $V$. We actually take the cloak as an
example. (c) gives the value range of chart $(U,\chi)$ and $(V,\psi_{\beta})$.
(d) gives the value range of chart $(U,\chi_{\gamma})$ and
$(V,\psi_{\alpha})$. Form invariance means the governing equations written by
the coordinates in (a) and (e), or (b) and (f), are the same.

First, the law of continuity [Eq.(1)] in $U$ under a curvilinear coordinate
system of $\mathsf{E}^{3}$ can be written as

$\nabla\cdot(\rho\mathbf{v})=\frac{1}{\sqrt{g}}\partial_{i}\left(\sqrt{g}\rho
v^{i}\right)=0.$ (A6)

Here $v^{i}$ is the contravariant component of velocity using contravariant
coordinates $\\{x^{i},x^{j},x^{k}\\}$ and the corresponding basis
$\\{\mathbf{g}_{i},\mathbf{g}_{j},\mathbf{g}_{k}\\}$. $g$ is the determinant
of the metric matrix $g_{ij}=\mathbf{g}_{i}\cdot\mathbf{g}_{j}$ which
satisfies $\sqrt{g}=\frac{1}{\operatorname{det}\operatorname{J}}$. The
Jacobian matrix is $\operatorname{J}_{ix}=\partial x_{i}/\partial x$, with
regard to the coordinate transformation from $\\{x,y,z\\}$ to
$\\{x^{i},x^{j},x^{k}\\}$. The form of Eq. (A6) in arbitrary curvilinear
coordinate system differs by only the value of $g$. This is the basic
requirement for building a TO theory. In particular, we can choose a
coordinate system which makes $\partial x_{i}/\partial x=\partial
x^{\prime}/\partial x$ since $f$ is a diffeomorphism. In fact, using
coordinates $\\{x^{i},x^{j},x^{k}\\}$ means constructing another atlas with a
single chart for $U$ (denoted by $(U,\chi_{\gamma})$), and there must exist a
homeomorphism satisfying $\chi_{\gamma}=\hat{f}\circ\chi=\psi_{\alpha}\circ
f$. Also, it’s easy to see $\chi_{\gamma}(U)=Y$ [Figs. A1(b) and A1(d)]. Then
we can obtain

$\partial_{x^{\prime}}\left(\dfrac{1}{\operatorname{det}\operatorname{J}}\rho
v^{x^{\prime}}(x^{\prime},y^{\prime},z^{\prime})\right)=0$ (A7)

from Eq. (A6). With
$v^{x^{\prime}}(x^{\prime},y^{\prime},z^{\prime})\text{J}^{-1}_{xx^{\prime}}=v^{x}(x,y,z)$,
Eq. (A7) can be rearranged as the law of continuity in $V$ [Eq. (5)] and the
transformation rule for velocity [Eq. (6)]. The key point here is that Eq.
(A6) can describe the law of continuity in different manifolds using the
charts $(U,\chi_{\gamma})$ and $(V,\psi_{\alpha})$ at the same time, which are
linked by the diffeomorphism $f$. Actually, the form invariance required by TO
in previous studies Science06-01  was usually explained by another two charts
including $(U,\chi)$. The left chart we haven’t mentioned is
$(V,\psi_{\beta})$, where $\psi_{\beta}=\hat{f}^{-1}\circ\psi_{\alpha}$ [Figs.
A1(e) and A1(c)]. We have $\psi_{\beta}(V)=X=\chi(U)$ and
$\chi=\psi_{\beta}\circ f$. If the law of continuity written by
$(V,\psi_{\beta})$ is exactly the same as its counterpart by $(U,\chi)$, we
can say it is form-invariant. By doing a coordinate transformation from
$(V,\psi_{\beta})$ to the partial Cartesian coordinates $(V,\psi_{\alpha})$,
the transformation rules for velocity can be obtained again. In fact, Figs.
1(b)–1(e) are drawn based on the latter way of understanding.

Second, we observe the form of Hele-Shaw equation [Eq. (9)] in $U$ using
contravariant coordinates $\\{x^{i},x^{j},x^{k}\\}$:

$\left(\dfrac{h^{2}}{12}(\partial_{i}p)g^{ij}+\mu
v^{j}\right)\mathbf{g}_{j}=0.$ (A8)

Do the same procedure as we have done to the law of continuity and we can
obtain

$\dfrac{1}{\operatorname{det}\text{J}}\left(\dfrac{h^{2}}{12}(\partial_{y^{\prime}}p^{\prime}(x^{\prime},y^{\prime},z^{\prime}))\text{J}_{y^{\prime}x}\text{J}_{xx^{\prime}}^{\top}+\mu
v^{x^{\prime}}\right)\mathbf{e}_{x}=0$ (A9)

for $(U,\chi_{\gamma})$ and $(V,\psi_{\alpha})$, which indicates the form
invariance. The pressure in $V$ should satisfy

$p^{\prime}(\mathbf{r}^{\prime})=p(f^{-1}(\mathbf{r}^{\prime})),$ (A10)

and the viscosity in $V$ should be transformed to
$\mu^{\prime}_{x^{\prime}y^{\prime}}(\mathbf{r}^{\prime})=\mu\det\text{J}\text{J}^{-\top}_{x^{\prime}x}\text{J}^{-1}_{xy^{\prime}}$
[Eq. (7)]. Another approach is directly using Eq. (10):

$\nabla\cdot\left(\dfrac{h^{2}}{12\mu}\nabla
p\right)=\frac{1}{\sqrt{g}}\partial_{j}\left(\sqrt{g}\dfrac{h^{2}}{12\mu}g^{ij}\partial_{i}p\right).$
(A11)

With $g^{ij}=\text{J}_{ix}\text{J}_{xj}^{\top}$, we can obtain Eq. (7) again.

If we directly consider the form invariance of the 3D model described by
$\mathbf{v}=-\frac{h^{2}}{8\mu}\left(1-\left(\frac{2z}{h}\right)^{2}\right)\nabla
p$ or
$\nabla\cdot\left(\frac{h^{2}}{8\mu}\left(1-\left(\frac{2z}{h}\right)^{2}\right)\nabla
p\right)=0$, we can get the same conclusion as Eq. (10) since they differ only
by a factor along the $z$-axis. Of course, this factor contains $z$. If the
coordinate transformation only happens in the horizontal plane, the 3D and 2D
models can share the same streamlines in the horizontal plane and the devices
we design, i.e., the cloak, the rotator, and the concentrator, can still work
in the 3D model.

###  IV.3 Anisotropic and inhomogeneous Hele-Shaw flows.

We have known the Hele-Shaw equation is form-invariant under coordinate
transformations and the transformation rule requires a rank-2 viscosity
tensor. In this note , we show that Eq. (11) is indeed the governing equation
for anisotropic inhomogeneous Stokes flows in shallow cells. Under lubrication
approximation, the strain rate with the Cartesian coordinates is

$\left(\nabla\mathbf{v}+\nabla\mathbf{v}^{\top}\right)\simeq\begin{bmatrix}{0}&{0}&{\frac{\partial
v_{x}}{\partial z}}\\\ {0}&{0}&{\frac{\partial v_{y}}{\partial z}}\\\
{\frac{\partial v_{x}}{\partial z}}&{\frac{\partial v_{y}}{\partial
z}}&{0}\end{bmatrix},$ (A12)

which comes from the fact that $\frac{\partial}{\partial
x}\ll\frac{\partial}{\partial z}$ and $\frac{\partial}{\partial
y}\ll\frac{\partial}{\partial z}$ InFlow  since the scale in the direction of
$z$-axis is much smaller than those in the $x$–$y$ plane. Similarly, the
Stokes equation [Eq. (4)] and the law of continuity [Eq. (1)] can respectively
be reduced to InFlow

$\frac{\partial p}{\partial x}=\frac{\partial\tau_{xz}}{\partial
z},\quad\frac{\partial p}{\partial y}=\frac{\partial\tau_{yz}}{\partial
z},\quad\frac{\partial p}{\partial z}=\frac{\partial\tau_{zz}}{\partial z}=0,$
(A13)

and

$\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{y}}{\partial y}=0.$ (A14)

Now we consider the case when the shear viscosity is anisotropic and
inhomogeneous in the $x$–$y$ plane, which can be expressed as

$\displaystyle\begin{bmatrix}\mu_{xy}\end{bmatrix}=\begin{bmatrix}\mu_{xx}(x,y)&\mu_{xy}(x,y)&0\\\
\mu_{yx}(x,y)&\mu_{yy}(x,y)&0\\\ 0&0&\mu_{zz}\end{bmatrix}.$ (A15)

Then the stress is

$\displaystyle\begin{bmatrix}\tau_{xy}\end{bmatrix}=\begin{bmatrix}0&0&\tau_{xz}\\\
0&0&\tau_{yz}\\\
\tau_{zx}&\tau_{zy}&0\end{bmatrix}=\begin{bmatrix}0&0&\mu_{xx}\frac{\partial
v_{x}}{\partial z}+\mu_{xy}\frac{\partial v_{y}}{\partial z}\\\
0&0&\mu_{yx}\frac{\partial v_{x}}{\partial z}+\mu_{yy}\frac{\partial
v_{y}}{\partial z}\\\ \mu_{zz}\frac{\partial v_{x}}{\partial
z}&\mu_{zz}\frac{\partial v_{y}}{\partial z}&0\end{bmatrix}.$ (A16)

We can notice that the pressure is only the function of $x$ and $y$ (recall
$\tau_{zz}$=0). Do an integral in Eq. (A13) and we have

$\displaystyle\begin{bmatrix}\mu_{xx}(x,y)&\mu_{xy}(x,y)\\\
\mu_{yx}(x,y)&\mu_{yy}(x,y)\end{bmatrix}\begin{bmatrix}\frac{\mathrm{d}v_{x}}{\mathrm{d}z}\\\
\frac{\mathrm{d}v_{y}}{\mathrm{d}z}\end{bmatrix}=\begin{bmatrix}\frac{\partial
p}{\partial x}\\\ \frac{\partial p}{\partial y}\end{bmatrix}z.$ (A17)

Finally we obtain the velocity distribution in the $x$–$y$ plane:

$\displaystyle\begin{bmatrix}v_{x}\\\
v_{y}\end{bmatrix}=-\frac{h^{2}}{8}\left(1-\left(\frac{2z}{h}\right)^{2}\right)\begin{bmatrix}\mu_{xx}(x,y)&\mu_{xy}(x,y)\\\
\mu_{yx}(x,y)&\mu_{yy}(x,y)\end{bmatrix}^{-1}\begin{bmatrix}\frac{\partial
p}{\partial x}\\\ \frac{\partial p}{\partial y}\end{bmatrix},$ (A18)

and the average values along the $z$ axis:

$\displaystyle\begin{bmatrix}\overline{v}_{x}\\\
\overline{v}_{y}\end{bmatrix}=-\frac{h^{2}}{12}\begin{bmatrix}\mu_{xx}(x,y)&\mu_{xy}(x,y)\\\
\mu_{yx}(x,y)&\mu_{yy}(x,y)\end{bmatrix}^{-1}\begin{bmatrix}\frac{\partial
p}{\partial x}\\\ \frac{\partial p}{\partial y}\end{bmatrix}.$ (A19)

In particular, when the shear viscosity is a scalar, we obtain the familiar
results standing for isotropic Hele-Shaw flows [Eq. (9)].

Here we use the anisotropic rank-2 shear viscosity in the derivations, which
can be proven consistent with the anisotropic rank-4 viscosity tensor
$A_{ijkl}$. Under the lubrication approximation, a simple form of $\tau_{xz}$
and $\tau_{yz}$ and can be obtained, writing

$\tau_{xz}=A_{xzxz}\frac{\partial v_{x}}{\partial z}+A_{xzyz}\frac{\partial
v_{y}}{\partial z},\quad\tau_{yz}=A_{yzxz}\frac{\partial v_{x}}{\partial
z}+A_{yzyz}\frac{\partial v_{y}}{\partial z}.$ (A20)

If we take

$\displaystyle\begin{bmatrix}\mu_{xx}&\mu_{xy}\\\
\mu_{yx}&\mu_{yy}\end{bmatrix}=\begin{bmatrix}A_{xzxz}&A_{xzyz}\\\
A_{yzxz}&A_{yzyz}\end{bmatrix},$ (A21)

the rank-4 viscosity can be reduced to the rank-2 shear viscosity. The bulk
viscosity disappears here due to the incompressibility [Eq. (A14)]. The zero
value of $\tau_{zz}=A_{zzxz}\frac{\partial v_{x}}{\partial
z}+A_{zzyz}\frac{\partial v_{y}}{\partial z}$ can be guaranteed by taking
$A_{zzxz}=A_{zzyz}=0$. A similar discussion on constructing such a
“generalized Darcy’s law” can also be found in a recent work arxiv22 , which
further considers the effect of gravity.

###  IV.4 Form variance of the Stokes equation and NS equations.

Here we investigate how the NS equations [Eq. (2)] and the Stokes equation
[Eq. (4)] change their form under curvilinear transformations. There are three
parts in Eq. (2): the pressure term, the inertial (or advection) term, and the
viscous stress term. With contravariant basis
$\left\\{\mathbf{g}^{i},\,\mathbf{g}^{j},\,\mathbf{g}^{k}\right\\}$ and
contravariant coordinates $\left\\{x^{i},\,x^{j},\,x^{k}\right\\}$, it’s easy
to see the pressure term is form-invariant:

$\nabla p=\dfrac{\partial p}{\partial x^{l}}\mathbf{g}^{l}.$ (A22)

The inertial term $(\mathbf{v}\cdot\nabla)\mathbf{v}$, which can be neglected
in Eq. (4), can be written as

$v^{i}\partial_{i}\left(v^{j}\mathbf{g}_{j}\right)=v^{i}\left(\partial_{i}v^{j}\right)\mathbf{g}_{j}+v^{i}v^{j}\partial_{i}\mathbf{g}_{j}.$
(A23)

The connection term $\partial_{i}\mathbf{g}_{j}$ is actually not the same in
different curvilinear coordinate systems. If this form change can’t be offset
by another in the divergence of viscous stress $\nabla\cdot{{\bm{\tau}}}$, the
NS equations can’t be form-invariant.

Now we check the form of
$\nabla\cdot{{\bm{\tau}}}=\nabla\cdot\mu\left(\nabla\mathbf{v}+\nabla\mathbf{v}^{\top}\right)$.
We can deal with this question from different perspectives.

Approach 1. The term $\nabla\cdot{\bm{\tau}}$ as a whole in curvilinear
coordinate systems is

$\displaystyle\nabla\cdot{\bm{\tau}}$
$\displaystyle=\mathbf{g}^{n}\cdot\dfrac{\partial}{\partial
x^{n}}\left(\tau^{ij}\mathbf{g}_{i}\otimes\mathbf{g}_{j}\right)$ (A24)
$\displaystyle=\dfrac{1}{\sqrt{g}}\partial_{i}\left(\sqrt{g}\tau^{ij}\right)\mathbf{g}_{j}+\tau^{ij}\Gamma^{n}_{ij}\mathbf{g}_{n}.$

Here the Christoffel symbol
$\Gamma^{n}_{ij}=\left(\partial_{i}\mathbf{g}_{j}\right)\cdot\mathbf{g}^{n}$
varies with the choice of coordinate system. Obviously, it and the connection
in inertial term do not cancel out. This is the easiest way to find the
general NS equations and the Stokes equation don’t meet the requirement of TO.

Approach 2. We assume that there is a rank-2 viscosity tensor $\bm{\mu}$,
although we’re not sure if it really exists. Then we focus on the generalized
form of Eq. (3):
$\nabla\cdot{\bm{\tau}}=\nabla\cdot\left({\bm{\mu}}\cdot\left(\nabla\mathbf{v}+\nabla\mathbf{v}^{\top}\right)\right)$
[Eq. (8)]. Recall
$\nabla\mathbf{v}\triangleq\bm{\epsilon}=\epsilon^{ij}\mathbf{g}_{i}\otimes\mathbf{g}_{j}$,
which is a little different from the strain rate tensor. We should notice
that, in the curvilinear coordinate system, $\epsilon^{ij}$ is not simply
equal to $\frac{\partial v^{l}}{\partial x^{j}}$, but involves Christoffel
symbols $\Gamma^{i}_{kl}$:

$\bm{\epsilon}=\nabla\mathbf{v}=g^{kj}\left({\frac{\partial v^{i}}{\partial
x^{k}}}+\Gamma^{i}_{kl}v^{l}\right)\mathbf{g}_{i}\otimes\mathbf{g}_{j}.$ (A25)

Finally we can find

$\nabla\cdot\left({\bm{\mu}}\cdot\left(\nabla\mathbf{v}+\nabla\mathbf{v}^{\top}\right)\right)=\dfrac{1}{\sqrt{g}}\partial_{i}\left(\sqrt{g}\mu^{ij}\left(\epsilon_{jm}+\epsilon_{mj}\right)\right)\mathbf{g}^{m}+\left(\mu^{ij}\left(\epsilon_{jm}+\epsilon_{mj}\right)\right)\partial_{i}\mathbf{g}^{m}.$
(A26)

Not surprisingly, the troublesome connection term still appears, and the
divergence of viscous stress is not form-invariant.

Approach 3. We start directly from Eq. (A1), which is expressed by the rank-4
viscosity tensor. In fact, writing the symmetric form
$\nabla\mathbf{v}+\nabla\mathbf{v}^{\top}$ in Eq. (A26) indicates
$A_{ijkl}=A_{ijlk}$ landau-6 . Again, we can rewrite the viscous stress term
using tensor components:

$\displaystyle\nabla\cdot{\bm{\tau}}$
$\displaystyle=\mathbf{g}^{n}\cdot\dfrac{\partial}{\partial
x^{n}}\left(A^{ijkl}\epsilon_{kl}\mathbf{g}_{i}\otimes\mathbf{g}_{j}\right)$
(A27)
$\displaystyle=\dfrac{1}{\sqrt{g}}\partial_{i}\left(\sqrt{g}A^{ijkl}\epsilon_{kl}\right)\mathbf{g}_{j}+A^{ijkl}\epsilon_{kl}\Gamma^{n}_{ij}\mathbf{g}_{n},$

and obtain the same conclusion as the previous two approaches.

Further, these approaches can help us find the conditions when TO can work in
Stokes flows. Since $\Gamma^{n}_{ij}=\Gamma^{n}_{ji}$, we can find one
possible condition to make Eqs. (A24) and (A27) look form-invariant:

$\tau^{ij}=-\tau^{ji}.$ (A28)

For a non-trivial case where $\bm{\epsilon}\neq\bm{0}$, we can see Eq. (A28)
actually requires

$A^{ijkl}=-A^{jikl}.$ (A29)

Under this condition, we can rewrite Eq. (A24) without the connection term as

$\dfrac{1}{\sqrt{g}}\partial_{i}\left(\sqrt{g}\tau^{ij}\right)\mathbf{g}_{j}-\dfrac{\partial
p}{\partial x^{k}}g^{kl}\mathbf{g}_{l}=0.$ (A30)

It seems difficult to judge the form invariance directly based on Eq. (A30).
Here, we apply the divergence operator on both sides of Eq. (A30) and have

$\partial_{j}\left(\sqrt{g}g^{ij}\partial_{i}p\right)+\partial_{j}\partial_{i}\left(\sqrt{g}\tau^{ij}\right)=0.$
(A31)

If the flow before transformation is incompressible, isotropic ad homogeneous
with a constant shear viscosity $\mu$, meaning $\tau^{ij}=\mu\epsilon^{ij}$,
then Eq. (A31) can be rearranged as

$\dfrac{1}{\sqrt{g}}\left(\partial_{j}\left(\dfrac{1}{\mu}\sqrt{g}g^{ij}\partial_{i}p\right)+\partial_{j}\partial_{i}\left(\sqrt{g}\epsilon^{ij}\right)\right)=0.$
(A32)

However, if $\mu\neq 0$, to satisfy Eq. (A28), we must require

$\epsilon^{ij}=-\epsilon^{ji},$ (A33)

and it’s easy to find the viscous term
$\partial_{j}\partial_{i}\left(\sqrt{g}\tau^{ij}\right)$ (or
$\partial_{j}\partial_{i}\left(\sqrt{g}\epsilon^{ij}\right)$) in Eq. (A31) (or
Eq. (A32)) should vanish (the so-called “drag-free” PRL2019 ). In addition,
Eq. (A33) can lead to $\epsilon_{ij}=-\epsilon_{ji}$ so the whole Eq. (A26)
also becomes a zero vector in this case. This is just a scenario when the
pressure satisfies the Laplace’s equation

$\nabla\cdot\nabla p=0,$ (A34)

or

$\nabla\cdot\left(\dfrac{1}{\mu}\nabla p\right)=0,$ (A35)

since $\mu$ is a constant. Eq. (A35) is mathematically form-invariant, and it
should be transformed to

$\nabla^{\prime}\cdot\left(\text{J}\text{J}^{\top}\det\text{J}^{-1}\dfrac{1}{\mu}\nabla^{\prime}p^{\prime}\right)=0$
(A36)

in the physical space. Since $g^{ij}=\text{J}_{ix}\text{J}_{xj}^{\top}$, we
can obtain the transformation rule for the viscosity again. However, to make
the governing equations reduce to Eq. (A36), we still have to seek for the
help of shallow geometries or the Darcy’s law in porous media PRL2011 .

###  IV.5 Deviations for invisible devices with different parameters

In this note, we compare the performance of the cloak, the rotator and the
concentrator when their corresponding transformations can have different
geometric parameters. Since the three devices are all expected to be visible,
we can use the following quantity to measure their invisibility effect:

$\frac{\left\langle\left|p^{\prime}-p\right|\right\rangle_{r^{\prime}>R_{2}}}{\Delta
p}=\frac{\int_{r^{\prime}>R_{2}}\left|p^{\prime}-p\right|\rm{d}r^{\prime}\rm{d}\theta^{\prime}}{\Delta
p\int_{r^{\prime}>R_{2}}\rm{d}r^{\prime}\rm{d}\theta^{\prime}}\approx\frac{\sum_{i=1}^{N}\left|p^{\prime}-p\right|}{N\Delta
p}.$ (A37)

The two formulas on the left side of the approximation symbol represent the
average absolute difference (deviation) between the background pressures in
the physical and virtual spaces. The formula on the right side gives how to
calculate this deviation based on the data of numerical simulations. $N$ is
the number of the nodes of finite-element method in the background area, i.e.,
$\\{\mathbf{r}:r>R_{2}\\}$ or $\\{\mathbf{r}^{\prime}:r^{\prime}>R_{2}\\}$.

Table 1: Average pressure deviation. We give the simulation results with the parameters in Fig. 2 and another set of parameters. | Cloak | Rotator | Concentrator
---|---|---|---
$h/L$ |  $R_{1}=2$ mm |  $R_{1}=1$ mm | $\theta=60\degree$ | $\theta=30\degree$ |  $R_{m}=2.5$ mm |  $R_{m}=2.9$ mm
$10^{-3}$ | $0.059\%$ | $0.0071\%$ | $0.041\%$ | $0.0056\%$ | $0.0015\%$ | $0.0061\%$
$10^{-2}$ | $0.59\%$ | $0.072\%$ | $0.31\%$ | $0.060\%$ | $0.0094\%$ | $0.035\%$
$10^{-1}$ | $4.0\%$ | $0.95\%$ | $1.7\%$ | $0.51\%$ | $0.10\%$ | $0.65\%$
$10^{-1}$111The governing equation Eq. (13) in $\\{\bm{r}^{\prime}:R_{1}\leqslant r^{\prime}\leqslant R_{2}\\}$ is replaced by $\nabla^{\prime}p^{\prime}+\frac{12}{h^{2}}{\bm{\mu}^{\prime}}\cdot\mathbf{v}^{\prime}=0$ [Eq. (12)]. | $0.019\%$ | $0.014\%$ | $0.058\%$ | $0.020\%$ | $0.0062\%$ | $0.034\%$

TABLE 1 shows the calculated deviations for the three invisible devices. Each
device corresponds to two columns of data, where the first column corresponds
to the results in Fig. 2, and the second column comes from a new set of
transformation parameters. For the cloak, we use a smaller $R_{1}$ in the new
simulations. For the rotator, we reduce the rotation angle $\theta$. The two
changes decrease the intensity of spatial transformations. Conversely, we use
a larger $R_{m}$ for the new concentrator which makes the transformation more
drastic. From TABLE 1, we can see the deviations increase with $h/L$ and some
even become non-negligible, indicating the failure of TO. The new parameters
reduce the deviations of the cloak and rotator while that of the new
concentrator becomes obvious now at $h/L=10^{-1}$. The last row of data gives
the deviations at $h/L=10^{-1}$ when the governing equation is the really
form-invariant one [Eq. (12)] in the device domain. Since $h/L=0.1$ is not
very large, such a modification is enough to generate negligible deviations.
When $h/L$ further becomes larger, we must use Eq. (12) in the whole
simulation domain to see the devices work well.

## References

  * (1) J. B. Pendry, D. Schurig, and D. R. Smith, Controlling electromagnetic fields, Science 312, 1780–1782 (2006).
  * (2) U. Leonhardt, Optical conformal mapping, Science 312, 1777–1780 (2006).
  * (3) A, Greenleaf, Y. Kurylev, M. Lassas, and G. Uhlmann, Cloaking devices, electromagnetic wormholes, and transformation optics, SIAM Rev. 51, 3–33 (2009).
  * (4) H. Chen, C. T. Chan, and P. Sheng, Transformation optics and metamaterials, Nature Mater. 9, 387–396 (2010).
  * (5) F. Sun, B. Zheng, H. Chen, W. Jiang, S. Guo, Y. Liu, Y. Ma, and S. He, Transformation optics: From classic theory and applications to its new branches, Laser Photonics Rev. 11, 1700034 (2017).
  * (6) M. Kadic, T. Bückmann, R. Schittny, and M. Wegener, Metamaterials beyond electromagnetism, Rep. Prog. Phys. 76, 126501 (2013).
  * (7) L. Xu and H. Chen, Transformation metamaterials, Adv. Mater. 33, 2005489 (2021).
  * (8) G. W. Milton, M. Briane, and J. R. Willis, On cloaking for elasticity and physical equations with a transformation invariant form, New J. Phys. 8, 248 (2006).
  * (9) L. D. Landau and E. M. Lifshitz, Fluid Mechanics, 2nd. (Butterworth-Heinemann, Oxford, 1987).
  * (10) S. A. Cummer and D. Schurig, One path to acoustic cloaking, New J. Phys. 9, 45 (2007).
  * (11) H. Chen and C. T. Chan, Acoustic cloaking in three dimensions using acoustic metamaterials, Appl. Phys. Lett. 91, 183518 (2007).
  * (12) Y. Chen, M. Zheng, X. Liu, Y. Bi, Z. Sun, P. Xiang, J. Yang, and G. Hu, Broadband solid cloak for underwater acoustics, Phys. Rev. B 95, 180104(R) (2017).
  * (13) H. Chen, J. Yang, J. Zi, and C. T. Chan, Transformation media for linear liquid surface waves, EPL 85, 24004 (2009).
  * (14) L. Han, S. Chen, and H. Chen, Water wave polaritons, Phys. Rev. Lett. 128, 204501 (2022).
  * (15) M. Farhat, S. Enoch, S. Guenneau, and A. B. Movchan, Broadband cylindrical acoustic cloak for linear surface waves in a fluid, Phys. Rev. Lett. 101, 134501 (2008).
  * (16) Y. A. Urzhumov and D. R. Smith, Fluid flow control with transformation media, Phys. Rev. Lett. 107, 074501 (2011).
  * (17) G. Dai, J. Shang, and J. Huang, Theory of transformation thermal convection for creeping flow in porous media: Cloaking, concentrating, and camouflage, Phys. Rev. E 97, 022129 (2018).
  * (18) J. Park, J. R. Youn, and Y. S. Song, Hydrodynamic metamaterial cloak for drag-free flow, Phys. Rev. Lett. 123, 074502 (2019).
  * (19) J. Park, J. R. Youn, and Y. S. Song, Metamaterial hydrodynamic flow concentrator, Extreme Mech. Lett., 42, 101061 (2021).
  * (20) J. Park, J. R. Youn, and Y. S. Song, Fluid-flow rotator based on hydrodynamic metamaterial, Phys. Rev. Applied 12, 061002 (2019).
  * (21) J. Park and Y. S. Song, Laminar flow manipulators, Extreme Mech. Lett. 40, 100908 (2020).
  * (22) For example, in the key step from Eq. (17) to Eq. (18) to obtain the form-invariance of the Stokes equation in the Supplementary Material of Ref. PRL2019 , spatially varying quantities (Jacobian and its determinant) and differential operators were directly swapped in order, which is mathematically forbidden.
  * (23) F. Tay, Y. Zhang, H. Xu, H. Goh, Y. Luo, and B. Zhang, A metamaterial-free fluid-flow cloak, Natl. Sci. Rev., nwab205 (2021).
  * (24) S. Yang, J. Wang, G. Dai, F. Yang, and J. Huang, Controlling macroscopic heat transfer with thermal metamaterials: Theory, experiment and application, Phys. Rep. 908, 1–65 (2021).
  * (25) O. Boyadjian, É. Boulais, and T. Gervais, Microfluidic surface shields: Control of flow and diffusion over sensitive surfaces, Phys. Rev. Applied 17, 014012 (2022).
  * (26) E. Boyko, V. Bacheva, M. Eigenbrod, F. Paratore, A. D. Gat, S. Hardt, and M. Bercovici, Microscale hydrodynamic cloaking and shielding via electro-osmosis, Phys. Rev. Lett. 126, 184502 (2021).
  * (27) H. S. Hele-Shaw, The flow of water, Nature 58, 34–36 (1898).
  * (28) H. A. Stone, A. D. Stroock, and A. Ajdari, Engineering flows in small devices: Microfluidics toward a lab-on-a-chip, Annu. Rev. Fluid Mech. 25, 381–411 (2004).
  * (29) A. D. Stroock, S. K. Dertinger, G. M. Whitesides, and A. Ajdari, Patterning flows using grooved surfaces, Anal. Chem. 74, 5306-5312 (2002).
  * (30) R. L. Panton, Incompressible Flow, 4th ed. (John Wiley & Sons, Inc., Hoboken, 2013).
  * (31) H. Chen and C. T. Chan, Transformation media that rotate electromagnetic fields, Appl. Phys. Lett. 90, 241105 (2007).
  * (32) M. Rahm, D. Schurig, D. A. Roberts, S. A. Cummer, D. R. Smith, and J. B. Pendry, Design of electromagnetic cloaks and concentrators using form-invariant coordinate transformations of Maxwell’s equations, Photonics Nanostruct. Fundam. Appl. 6, 87 (2008).
  * (33) Y. Huang, Y. Feng, and T. Jiang, Electromagnetic cloaking by layered structure of homogeneous isotropic materials, Opt. Express 15, 11133-11141 (2007).
  * (34) H. Chen and C. T. Chan, Electromagnetic wave manipulation by layered systems using the transformation media concept, Phys. Rev. B 78, 054204 (2008).
  * (35) C. Navau, J. Prat-Camps, and A. Sanchez, Magnetic energy harvesting and concentration at a distance by transformation optics, Phys. Rev. Lett. 109, 263903 (2012).
  * (36) R. Tao and K. Huang, Reducing blood viscosity with magnetic ﬁelds, Phys. Rev. E 84, 011905 (2011).
  * (37) B. Wang, T.-M. Shih, L. Xu, G. Dai, and J. Huang, Intangible hydrodynamic cloaks for convective flows, Phys. Rev. Applied 15, 034014 (2021).
  * (38) T. T. Al-Housseiny and H. A. Stone, Controlling viscous fingering in tapered Hele-Shaw cells, Phys. Fluids 25, 092102 (2013).
  * (39) E. O. Dias and J. A. Miranda, Taper-induced control of viscous fingering in variable-gap Hele-Shaw flows, Phys. Rev. E 87, 053015 (2013).
  * (40) A. B. Thompson, A. Juel, and A. L. Hazel, Multiple finger propagation modes in Hele-Shaw channels of variable depth, J. Fluid Mech. 746, 123–164 (2014).
  * (41) Y. Hua, C. Qian, H. Chen, and H. Wang, Experimental topology-optimized cloak for water waves, Mater. Today Phys. 27, 100754 (2022).
  * (42) X. Zhao, X. Hu, and J. Zi, Fast water waves in stationary surface disk arrays, Phys. Rev. Lett. 127, 254501 (2021).
  * (43) R. D. Smith and J. C. McWilliams, Anisotropic horizontal viscosity for ocean models, Ocean Modell. 5, 129–156 (2003).
  * (44) T. Khain, C. Scheibner, M. Fruchart, and V. Vitelli, Stokes flows in three-dimensional fluids with odd and parity-violating viscosities, J. Fluid Mech. 934, A23 (2022).
  * (45) J. Jost, Riemannian Geometry and Geometric Analysis, 7th ed. (Springer, Cham, 2017).
  * (46) D. Reynolds, G. M. Monteiro, and S. Ganeshan, Generalized Darcy’s law for parity odd three-dimensional fluids, arXiv:2112.03076.