11institutetext:  Information and Communications, Nagoya University,
Furo-cho, Chikusa-ku, Nagoya, Aichi 4648601, Japan
11email<EMAIL_ADDRESS>22institutetext: Graduate School of
Informatics, Nagoya University, Aichi, Japan 33institutetext: Graduate School
of Science and Technology, Nara Institute of Science and Technology, Nara,
Japan 44institutetext: Research Center for Medical Bigdata, National Institute
of Informatics, Tokyo, Japan 55institutetext: Department of Radiology, Keio
University School of Medicine, Tokyo, Japan 66institutetext: Department of
Radiology, Juntendo University, Tokyo, Japan

# COVID-19 Infection Segmentation from Chest CT Images Based on Scale
Uncertainty

Masahiro Oda 1122 Tong Zheng 22 Yuichiro Hayashi 22 Yoshito Otake 3344
Masahiro Hashimoto 55 Toshiaki Akashi 66 Shigeki Aoki 66 Kensaku Mori 221144

###### Abstract

This paper proposes a segmentation method of infection regions in the lung
from CT volumes of COVID-19 patients. COVID-19 spread worldwide, causing many
infected patients and deaths. CT image-based diagnosis of COVID-19 can provide
quick and accurate diagnosis results. An automated segmentation method of
infection regions in the lung provides a quantitative criterion for diagnosis.
Previous methods employ whole 2D image or 3D volume-based processes. Infection
regions have a considerable variation in their sizes. Such processes easily
miss small infection regions. Patch-based process is effective for segmenting
small targets. However, selecting the appropriate patch size is difficult in
infection region segmentation. We utilize the scale uncertainty among various
receptive field sizes of a segmentation FCN to obtain infection regions. The
receptive field sizes can be defined as the patch size and the resolution of
volumes where patches are clipped from. This paper proposes an infection
segmentation network (ISNet) that performs patch-based segmentation and a
scale uncertainty-aware prediction aggregation method that refines the
segmentation result. We design ISNet to segment infection regions that have
various intensity values. ISNet has multiple encoding paths to process patch
volumes normalized by multiple intensity ranges. We collect prediction results
generated by ISNets having various receptive field sizes. Scale uncertainty
among the prediction results is extracted by the prediction aggregation
method. We use an aggregation FCN to generate a refined segmentation result
considering scale uncertainty among the predictions. In our experiments using
199 chest CT volumes of COVID-19 cases, the prediction aggregation method
improved the dice similarity score from 47.6% to 62.1%.

###### Keywords:

COVID-19 Infection segmentation Scale uncertainty.

##  1 Introduction

Novel coronavirus disease 2019 (COVID-19) spread worldwide, causing many
infected patients and deaths. The total number of cases and deaths related to
COVID-19 are more than 212 million and 4.4 million in the world [1]. Because
of the rapid increase of COVID-19 patients, medical institutions suffer from a
human resources shortage. To prevent further infection, a quick inspection
method for COVID-19 infection is pressing required. Such quick inspection
enables providing appropriate treatments to patients and curbs the spread of
COVID-19. Reverse transcriptase-polymerase chain reaction (RT-PCR) testing is
used as an inspection method of COVID-19 cases. However, it takes some hours
to give a diagnosis result. Furthermore, its sensitivity is not high, ranging
from 42% to 71% [2]. As another choice of COVID-19 cases, CT image-based
diagnosis is helpful. The sensitivity of CT image-based COVID-19 diagnosis is
reported as 97% [3]. Furthermore, a CT scan takes only some minutes. A CT
image-based computer-aided diagnosis (CAD) system for COVID-19 is expected to
provide a quick and accurate diagnosis to patients. For such CAD systems, a
quantitative analysis method of the lung condition is essential. Ground-glass
opacities (GGOs) and consolidations are commonly found in the lung of viral
pneumonia cases, including COVID-19. We call them infection regions. Automatic
segmentation of them is an essential component of CAD systems.

Related Work of COVID-19 Segmentation: Previously, deep learning-based
automatic segmentation methods of infection regions from CT volumes of
COVID-19 cases were proposed [4, 5, 6, 7, 8]. Fan et al. [4] proposed an
infection region segmentation method using the Inf-Net. The Inf-Net utilizes
reverse attention and edge attention to learn features to differentiate
infection and other regions. However, because they employ 2D image-based
process, 3D positional information is not utilized in their segmentation
method. Other papers also employ 2D image-based process [5, 6, 7]. Yan et al.
[8] proposed a fully convolutional network (FCN) to segment infection and
normal regions in the lung. The FCN has contrast enhancement branches to
extract features of infection regions that have various intensities. However,
because contrast information of segmentation targets is not explicitly
provided to the FCN, the contrast enhancement branches’ contribution to
improving segmentation accuracy is limited.

Scale Uncertainty on Patch-based Process: Infection regions contain many small
regions. Segmentation processes hardly segment them from whole 2D slice image
or 3D volume as performed in the previous methods[4, 8]. To segment such small
regions, a patch-based approach is practical. Patch-based approach is commonly
employed in segmentation methods from images of large data size such as 3D CT
volume [9, 10, 11] or pathological images [12, 13, 14, 15]. The approach is
advantageous to perform deep learning-based segmentation under the limitation
of GPU memory size. In patch-based approaches, patch size is an essential
factor of segmentation accuracy. The patch size defines the size of the
receptive field of segmentation models. Also, original images or volumes can
be scaled before patch clipping to change the receptive field size. In
summary, (a) the resolution of original volume (VRes) and (b) the size of
patch (PSize) are essential factors for the segmentation accuracy in patch-
based approaches. In a multi-organ segmentation method [9], the use of a
relatively large PSize resulted in the achievement of high segmentation
accuracies among large organs (liver, spleen, and stomach). However, their
segmentation accuracy of small organs (artery, vein, and pancreas) was low.
Other paper [10] reported that the use of small PSize is effective for small
organ (artery) segmentation. VRes and PSize should be selected to patch covers
the segmentation target from their results.

In infection region segmentation, selecting appropriate VRes and PSize are
difficult because the sizes of infection regions are different for each
region. If we apply a segmentation process using multiple VRess and PSizes, we
can obtain multiple prediction maps having variation among them. The variation
can be considered as uncertainty among scales. The scale uncertainty
represents useful information to obtain an accurate segmentation result. Scale
uncertainty-aware aggregation process of multiple prediction maps is essential
for segmenting infection regions with various sizes.

Proposed Method and Contributions: We present an infection region segmentation
method from a chest CT volume of a COVID-19 patient. We developed a patch-
based FCN for infection region segmentation called infection segmentation
network (ISNet) to perform segmentation. Also, we propose a scale uncertainty-
aware aggregation method of prediction results. These methods enable the
segmentation of infection regions of various sizes. ISNet has multiple encoder
and a single decoder style structure. The use of the multiple encoders enables
feature extraction from infection regions with a variation of CT values. Deep
supervision is employed to improve the decoder’s ability to decode the
prediction result from feature value. ISNets having various receptive field
sizes are trained and used to generate prediction maps from the CT volume. The
scale uncertainty-aware prediction aggregation is applied to the multiple
prediction maps to generate a final segmentation result considering
uncertainty among the prediction results related to the receptive fields’
size.

The contributions of this paper are (1) proposal of the ISNet with multiple
encoders for feature extraction from infection regions that have a variation
of CT values and (2) proposal of the scale uncertainty-aware aggregation
method of prediction maps that are generated by segmentation models having a
various size of receptive fields. These methods improve the segmentation
accuracy of targets with significant variations in their intensity values and
sizes.

##  2 Method

The proposed method segments infection regions from a chest CT volume of a
COVID-19 patient. Set of patch volumes clipped from a CT volume is provided to
ISNet. VRes and PSize define the size of the receptive field of ISNet. Change
of the receptive field size causes variation on segmentation results (scale
uncertainty). The scale uncertainty contains valuable information to refine
segmentation results. We propose a scale uncertainty-aware aggregation process
of segmentation results, which ISNets segment on various VRess and PSizes. The
process generates a refined segmentation result.

###  2.1 Infection Region Segmentation by ISNet

Overview of Model: The structure of ISNet is shown in Fig. 1. ISNet has
multiple encoders and a single decoder. Multiple volumes are generated from an
input CT volume by applying CT value normalization by multiple value ranges to
improve the segmentation accuracy of infection regions with various CT values.
Patch volumes clipped from the volumes are input to ISNet. ISNet has multiple
encoders corresponding to the multiple inputs to extract features in the CT
value ranges selectively. The encoder has dense pooling connections [16] that
prevent loss of spatial information by pooling layers. We employ deep
supervision [17, 18] in the decoder to improve its decoding performance from
feature values.

Figure 1: Model structure of ISNet. Two encoders process patch volumes
generated by different intensity normalizations. Encoders have dense pooling
connections to bottleneck layer. Deep supervision is employed to evaluate
subscale outputs.

Multiple Range Normalized Patch: An input CT volume is converted to a volume
having an isotropic resolution in three dimensions. Then, the volume is scaled
to $v\times v\times z$ voxels maintaining the aspect ratio. The number of
voxels along the body axis $z$ differs for each CT volume depending on its
scanning range. CT values of infection regions distribute widely. CT values of
consolidations range from -300 to 100 H.U.. GGO has a lower and broader range
of CT values than the consolidations, ranging from -800 to 0 H.U.. CT value
normalizations by multiple ranges are adequate for such a target. We apply CT
value normalizations to the scaled CT volume using two CT value ranges,
including; wide range (WRange): [-1000 H.U., 950 H.U.] and narrow range
(NRange): [-1000 H.U., -400 H.U.]. Normalization results by the WRange are
suitable to observe high-intensity infection regions, including consolidations
and GGOs having high intensities. Normalization results by the NRange are
suitable to observe GGOs having low intensities. Samples of normalization
results are shown in Fig. 2. Two normalized volumes, including WRange volume
and NRange volume, are generated from this process. We clip patch volumes from
them at random positions. Patch volumes clipped from the WRange and NRange
volumes are described as ${\bf I}^{W}_{v,p},{\bf
I}^{N}_{v,p}\in\mathbb{R}^{p\times p\times p}$, respectively.

Figure 2: Multiple range normalization results. Visibility of high- and low-
intensity infection regions are high in WRange and NRange normalization
result, respectively.

Multiple Encoding Paths: Inputs of the ISNet are patch volumes. We use two
independent encoders to process ${\bf I}^{W}_{v,p}$ and ${\bf I}^{N}_{v,p}$.
Feature values extracted by the encoders are concatenated at the bottleneck
layer.

Pooling layers are commonly used in encoder, although it reduces spatial
information in feature maps. The bottleneck layer connected after the encoder
cannot receive enough spatial information. It causes segmentation results
having incorrect boundaries. To reduce the loss of spatial information in the
encoder, we adopt dense pooling connections [16] in the two encoders of ISNet.
The dense pooling connections provide spatial information at each resolution
in the encoder to the bottleneck layer. In the dense pooling connections,
mixed poolings [16] are used instead of max poolings to reduce the loss of
spatial information. Furthermore, we use dilated convolution [19] to utilize
sparsely-distributed features in convolution operations. Dilated convolution
block was implemented by connecting dilated convolutions of multiple dilation
rates in parallel to obtain multiple-scales convolution results. Some dilated
convolution blocks are inserted into ISNet.

Training: ISNet estimates a patch prediction volume ${\bf
P}_{v,p}\in\mathbb{R}^{p\times p\times p}$ of infection regions from two input
patch volumes ${\bf I}^{W}_{v,p}$ and ${\bf I}^{N}_{v,p}$. ISNet that performs
estimation from input patch volumes ($p\times p\times p$ voxels) clipped from
a volume ($v\times v\times z$ voxels) can be represented as a function
$f_{v,p}$. Estimation of a patch prediction volume is formulated as

${\bf P}_{v,p}=f_{v,p}({\bf I}^{W}_{v,p},{\bf
I}^{N}_{v,p};\boldsymbol{\theta}_{v,p}),$ (1)

where $\boldsymbol{\theta}_{v,p}$ is a parameter vector for infection region
segmentation. The parameter vector is optimized in a supervised training
process using CT volumes for training and their corresponding ground truth
volumes ${\bf G}\in\\{0,1\\}^{p\times p\times p}$, whose elements 1 and 0
represent voxels in target or background regions. We employ deep supervision
[17, 18] for two subscales outputs. Their patch prediction volumes are ${\bf
P}_{v,p}^{(1)}$ and ${\bf P}_{v,p}^{(2)}$. Their sizes are magnified to the
same size as ${\bf P}_{v,p}$. The loss function to train ISNet is defined as

$L=Dice({\bf G},{\bf P}_{v,p})+\sum^{2}_{i=1}Dice({\bf G},{\bf
P}_{v,p}^{(i)}),$ (2)

where $Dice$ is the dice loss between the ground truth volume and the patch
prediction volumes.

Prediction: Patch volumes clipped from a CT volume for prediction are given to
the trained ISNet $f_{v,p}$. The resulting patch prediction volumes are
reconstructed as the same size as the CT volume. The reconstructed prediction
volume is denoted as ${\bf R}_{v,p}\in\mathbb{R}^{v\times v\times z}$.

###  2.2 Scale Uncertainty-Aware Prediction Aggregation

The parameters $v$ and $p$ define the size of the receptive field of ISNet.
The size of the receptive field of ISNet has a relationship to its
segmentation accuracy. ISNets having various sizes of their receptive fields
are trained and perform predictions, and we obtain multiple prediction volumes
containing scale uncertainty from them. We utilize the scale uncertainty-aware
aggregation method of the prediction volumes. An aggregation function is
automatically trained based on each prediction volume’s contribution to a
segmentation result.

Figure 3: Model structure of aggregation FCN. FCN processes axial slices
obtained from prediction volumes. FCN outputs aggregation results from them.

We train ISNets using training cases on multiple value settings of $v$ and
$p$. Using the ISNets, multiple reconstructed prediction volumes are generated
from a CT volume for prediction. We perform aggregation of them using an
aggregation FCN. The structure of the aggregation FCN is shown in Fig. 3. The
FCN employs multiple 2D axial slice-based process. The FCN combines given
prediction results by considering variation among them and how each prediction
result contributes to generating a segmentation result. The FCN is trained
using prediction volumes and their corresponding ground truth volumes of
training cases. Generalized dice loss [20] is used to train the FCN. In a
prediction process, outputs (on 2D axial slices) of the trained aggregation
FCN are reconstructed to a 3D volume that has the same resolution and size as
the original CT volume and thresholded to generate a segmentation result.

##  3 Experiments and Results

We evaluated segmentation accuracy of the proposed method. We used 199 chest
CT volumes of COVID-19 patients provided by the Multi-national NIH Consortium
for CT AI in COVID-19 [21] via the NCI TCIA public website [22]. The
corresponding ground truth data of infection regions were also provided. We
conducted three-fold cross-validations in our evaluations. Averaged precision,
recall, and dice score among all CT volumes are used as the evaluation
criteria. Methods were implemented using Keras 2.2.4 and TensorFlow 1.14.0.
NVIDIA Tesla V100 GPU having 32GB memory $\times$ 1 was used to train and test
the methods.

###  3.1 Ablation and Comparative Study of ISNet

We confirmed the segmentation performance of infection regions using ISNet.
Techniques explained in 2.1 are used, including the deep supervision (DS) and
the multiple encoding paths (ME). We confirmed the effectiveness of the
technique in an ablation study. Dice scores of segmentation results obtained
using ISNet and ISNet without DS and ME were calculated. ISNets were trained
on parameter settings of $v=192$, $p=32$, minibatch size: 16, learning rate:
$10^{-5}$, and training epochs: 40. Adam was used as the optimization
algorithm. The segmentation result was generated by applying thresholding to
the prediction volume. Also, we compared dice scores of previous methods,
including 3D U-Net [23] and 3D U-Net having squeeze-and-excitation (SE) blocks
[24, 25] with ISNet. These previous methods were applied to perform patch-
based processes. The results are shown in Table 1. ISNet had a higher dice
score than the previous methods. Also, the use of DS and ME contributed to
improving the dice score.

Table 1: Segmentation accuracies of ISNet using deep supervision (DS) and multiple encoding paths (ME). Accuracies of previous methods are also shown. Method | Precision (%) | Recall (%) | Dice (%)
---|---|---|---
ISNet (proposed method) | 58.7 | 54.6 | 56.6
ISNet without DS and ME | 51.1 | 58.7 | 54.6
3D U-Net having SE blocks [24, 25] | 52.3 | 59.0 | 55.5
3D U-Net [23] | 52.5 | 55.6 | 54.0

###  3.2 Segmentation by Aggregation FCN

We applied the scale uncertainty-aware prediction aggregation to prediction
volumes generated by ISNets. 11 prediction volumes generated by ISNets trained
on parameter settings $(v,p)\in\\{(256,32),$ (256,64), (224,32), (224,64),
(192,32), (192,64), (160,32), (160,64), (128,32), (128,64), $(96,32)\\}$ were
used. The aggregation FCN was trained on parameter settings of minibatch size:
4, learning rate: $10^{-3}$, and training epochs: 5. Adam was used as the
optimization algorithm. Generated prediction volumes and aggregation results
from them are shown in Fig. 4. Segmentation accuracies of ISNets and
aggregation result are shown in Table 2. Accuracies were improved in all
criteria by the aggregation.

Figure 4: Prediction volumes generated by ISNets with parameters $(v,p)$ and aggregation results from them. While regions in prediction volumes have considerable variation, they are adequately aggregated. Table 2: Segmentation accuracies of ISNets and aggregation result from them. Method | Precision (%) | Recall (%) | Dice (%)
---|---|---|---
Aggregation FCN | 63.7 | 60.5 | 62.1
Mean $\pm$ S.D. of 11 ISNets (before aggregation) | $45.1\pm 12.3$ | $52.6\pm 3.7$ | $47.6\pm 8.1$

##  4 Discussion and Conclusions

Segmentation of the lung infection region is difficult because it has
significant variations in its CT values and shapes. We developed ISNet with
the multiple range normalized patch processing paths and the scale
uncertainty-aware prediction aggregation process to tackle infection
segmentation having such difficulties. ISNet achieved higher segmentation
accuracy than the previous methods, as shown in Table 1. Also, we confirmed
the effectiveness of techniques, including deep supervision and multiple
encoding paths in the ablation study. The scale uncertainty-aware prediction
aggregation process improved the dice similarity score of the segmentation
result. We used multiple prediction volumes generated by using multiple ISNets
with various receptive fields’ sizes. Because segmentation abilities and
effective segmentation target sizes are different among the ISNets, an
appropriate aggregation process of the prediction volumes can generate an
accurate segmentation result. The aggregation FCN with trainable aggregation
parameters was successfully built using training data. The evaluation result
obtained in the cross-validation proved that the trained aggregation FCN has a
high generalization ability to perform segmentation from prediction volumes.

This paper proposed a segmentation method of infection regions in the lung
from a CT volume of a COVID-19 patient. To segment infection regions having
variations of CT value and size, we proposed ISNet and the scale uncertainty-
aware prediction aggregation process. In our experiments, the aggregation
process improved segmentation accuracy from individual ISNet results. Future
work includes increasing variety of the receptive field sizes to process in
the prediction aggregation process and development of a CAD system for
COVID-19 diagnosis.

#### Acknowledgements.

Parts of this research were supported by the AMED Grant Numbers
18lk1010028s0401, JP19lk1010036, JP20lk1010036, JP20lk1010025, the NICT Grant
Number 222A03, the JST CREST Grant Number JPMJCR20D5, and the MEXT/JSPS
KAKENHI Grant Numbers 26108006, 17H00867, 17K20099.

## References

  * [1] Coronavirus Update, https://www.worldometers.info/coronavirus/. Last accessed 22 August 2021
  * [2] Simpson, S., Kay, F. U., Abbara, S., Bhalla, S., Chung, J. H., Chung, M., Henry, T. S., Kanne, J. P., Kligerman, S., Ko, J. P., Litt, H.: Radiological Society of North America Expert Consensus Document on Reporting Chest CT Findings Related to COVID-19: Endorsed by the Society of Thoracic Radiology, the American College of Radiology, and RSNA. Radiology: Cardiothoracic Imaging 2(2) (2020)
  * [3] Ai, T., Yang, Z., Hou, H., Zhan, C., Chen, C., Lv, W., Tao, Q., Sun, Z., Xia, L.: Correlation of Chest CT and RT-PCR Testing for Coronavirus Disease 2019 (COVID-19) in China: A Report of 1014 Cases. Radiology 296(2) (2020)
  * [4] Fan, D.-P., Zhou, T., Ji, G.-P., Zhou, Y., Chen, G., Fu, H., Shen, J., Shao, L.: Inf-Net: Automatic COVID-19 Lung Infection Segmentation From CT Images. IEEE Transactions on Medical Imaging 39(8), 2626–2637 (2020)
  * [5] Wang, G., Liu, X., Li, C., Xu, Z., Ruan, J., Zhu, H., Meng, T., Li, K., Ning, H., Zhang, S.: A Noise-Robust Framework for Automatic Segmentation of COVID-19 Pneumonia Lesions From CT Images. IEEE Transactions on Medical Imaging 39(8), 2653–2663 (2020)
  * [6] Zheng, B., Liu, Y., Zhu, Y., Yu, F., Jiang, T., Yang, D., Xu, T.: MSD-Net: Multi-Scale Discriminative Network for COVID-19 Lung Infection Segmentation on CT. IEEE Access 8, 185786–185795 (2020)
  * [7] Mahmud, T., Alam, M.J., Chowdhury, S., Ali, S.N., Rahman, M.M., Fattah, S.A., Saquib, M.: CovTANet: A Hybrid Tri-level Attention Based Network for Lesion Segmentation, Diagnosis, and Severity Prediction of COVID-19 Chest CT Scans. IEEE Transactions on Industrial Informatics (Early Access) (2020)
  * [8] Yan, Q., Wang, B., Gong, D., Luo, C., Zhao, W., Shen, J., Shi, Q., Jin, S., Zhang, L., You, Z.: COVID-19 Chest CT Image Segmentation - A Deep Convolutional Neural Network Solution. arXiv:2004.10987 (2020)
  * [9] Roth, H. R., Oda, H., Zhou, X., Shimizu, N., Yang, Y., Hayashi, Y., Oda, M., Fujiwara, M., Misawa, K., Mori, K.: An Application of Cascaded 3D Fully Convolutional Networks For Medical Image Segmentation. Computerized Medical Imaging and Graphics 66, 90–99 (2018)
  * [10] Oda, M., Roth, H.R., Kitasaka, T. et al.: Abdominal Artery Segmentation Method From CT Volumes Using Fully Convolutional Neural Network. International Journal of Computer Assisted Radiology and Surgery 14, 2069–2081 (2019)
  * [11] Kim, H., Jung, J., Kim, J. et al.: Abdominal Multi-organ Auto-segmentation Using 3D-patch-based Deep Convolutional Neural Network. Scientific Reports 10, 6204 (2020)
  * [12] Zhou, Y., Chang, H., Barner, K., Spellman, P., Parvin, B.: Classification of Histology Sections via Multispectral Convolutional Sparse Coding. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 3081–3088 (2014)
  * [13] Xu, Y., Jia, Z., Ai, Y., Zhang, F., Lai, M., Chang, E. I.: Deep Convolutional Activation Features For Large Scale Brain Tumor Histopathology Image Classification And Segmentation. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 947–951 (2015)
  * [14] Wang, D., Khosla, A., Gargeya, R., Irshad, H., Beck, A. H.: Deep Learning for Identifying Metastatic Breast Cancer. arXiv:1606.05718 (2016)
  * [15] Tokunaga, H., Teramoto, Y., Yoshizawa, A., Bise, R.: Adaptive Weighting Multi-field-of-view CNN For Semantic Segmentation in Pathology. IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 12589–12598 (2019)
  * [16] Playout, C., Duval, R., Cheriet, F.,: A Multitask Learning Architecture For Simultaneous Segmentation of Bright And Red Lesions in Fundus Images. MICCAI, LNCS, vol. 11071, pp.101–108 (2018)
  * [17] Zeng, G., Yang, X., Li, J., Yu, L., Heng, P. A., Zheng, G.: 3D U-net with Multi-level Deep Supervision: Fully Automatic Segmentation of Proximal Femur in 3D MR Images. Machine Learning in Medical Imaging (MLMI) LNCS, vol. 10541, pp.274–282 (2017)
  * [18] Dou, Q., Yu, L., Chen, H., Jin, Y., Yang, X., Qin, J., Heng, P.-A.: 3D Deeply Supervised Network For Automated Segmentation of Volumetric Medical Images. Medical Image Analysis 41, 40–54 (2017)
  * [19] Yu, F. Koltun, V.,: Multi-scale Context Aggregation by Dilated Convolutions. International Conference on Learning Representations (ICLR) (2016).
  * [20] Sudre, C. H., Li, W., Vercauteren, T., Ourselin, S., Cardoso, M. J.: Generalised Dice Overlap as a Deep Learning Loss Function for Highly Unbalanced Segmentations. International Workshop on Deep Learning in Medical Image Analysis (DLMIA) 10553, 240–248 (2017)
  * [21] An, P., Xu, S., Harmon, S. A., Turkbey, E. B., Sanford, T. H., Amalou, A., Kassin, M., Varble, N., Blain, M., Anderson, V., Patella, F., Carrafiello, G., Turkbey, B. T., Wood, B. J.: CT Images in Covid-19 [Data set]. The Cancer Imaging Archive (2020)
  * [22] Clark, K., Vendt, B., Smith, K., Freymann, J., Kirby, J., Koppel, P., Moore, S., Phillips, S., Maffitt, D., Pringle, M., Tarbox, L., Prior, F.: The Cancer Imaging Archive (TCIA): Maintaining and Operating a Public Information Repository. Journal of Digital Imaging 26(6), 1045–1057 (2013)
  * [23] Çiçek, Ö., Abdulkadir, A., Lienkamp, S. S., Brox, T., Ronneberger, O: 3D U-Net: Learning Dense Volumetric Segmentation from Sparse Annotation. MICCAI, LNCS, vol. 9901, pp.424–432 (2016)
  * [24] Roy, A. G., Navab, N., Wachinger, C.: Concurrent Spatial and Channel Squeeze & Excitation in Fully Convolutional Networks. MICCAI, LNCS, vol. 11070, pp.421–429 (2018)
  * [25] Rundo, L., Han, C., Nagano, Y., Zhang, J., Hataya, R., Militello, C., Tangherloni, A., Nobile, M. S., Ferretti, C., Besozzi, D., Gilardi, M. C., Vitabile, S., Mauri, G., Nakayama, H., Cazzaniga, P.: USE-Net: Incorporating Squeeze-and-Excitation blocks into U-Net for prostate zonal segmentation of multi-institutional MRI datasets. Neurocomputing 365, pp.31–43 (2019)


# Coherent structures in plane channel flow of dilute polymer solutions with
vanishing inertia

Alexander Morozov<EMAIL_ADDRESS>SUPA, School of Physics and
Astronomy, The University of Edinburgh, James Clerk Maxwell Building, Peter
Guthrie Tait Road, Edinburgh, EH9 3FD, United Kingdom

###### Abstract

When subjected to sufficiently strong velocity gradients, solutions of long,
flexible polymers exhibit flow instabilities and chaotic motion, often
referred to as elastic turbulence. Its mechanism differs from the familiar,
inertia-driven turbulence in Newtonian fluids, and is poorly understood. Here,
we demonstrate that the dynamics of purely elastic pressure-driven channel
flows of dilute polymer solutions are organised by exact coherent structures
that take the form of two-dimensional travelling waves. Our results
demonstrate that no linear instability is required to sustain such travelling
wave solutions, and that their origin is purely elastic in nature. We show
that the associated stress profiles are characterised by thin, filament-like
arrangements of polymer stretch, which is sustained by a solitary pair of
vortices. We discuss the implications of the travelling wave solutions for the
transition to elastic turbulence in straight channels, and propose ways for
their detection in experiments.

Close to the transition to chaos, a wide range of systems, spanning
magnetoelastic devices [1], electric circuits [2], chemical reactions [3],
flows of Newtonian fluids [4, 5], neural and myocardial [6] tissues, and
lasers [7], exhibits low-dimensional dynamics organised by unstable coherent
structures. The phase space vicinity of such structures is generally
attractive, with a (relatively) small number of unstable directions. In this
scenario, a chaotic trajectory amounts to a ‘pin-ball’-like motion in phase
space that spends a long time in close vicinity of coherent structures.
Knowledge of such structures thus allows for relatively accurate statistical
description of an otherwise unpredictable system [8].

In Newtonian parallel shear flows, the discovery of travelling waves and
periodic orbits organised around streamwise streaks and vortices has lead to
the prediction of the transitional values of the Reynolds number [9], revealed
the transient nature of turbulence near its onset [10, 11], and helped to
identify directed percolation as the universality class for the transition to
Newtonian turbulence [12, 13, 14]. Addition of small amounts of high-molecular
weight flexible polymers was shown to suppress Newtonian coherent structures
in parallel shear flows [15, 16], while higher levels of viscoelasticity
trigger elasto-inertial turbulence [17, 18, 19, 20], characterised by a
different set of coherent structures [21, 22, 23, 24].

At very low Reynolds numbers, in the presence of sufficiently strong velocity
gradients, dilute solutions of high-molecular weight flexible polymer
molecules exhibit a unique chaotic flow state often referred to as elastic
turbulence [25, 26]. Its origins differ from the inertia-dominated turbulence
in Newtonian fluids, and lie in the flow-induced stretch and orientation of
polymer molecules [27]. Despite being implicated as a major production-
limiting factor in polymer processing [28], elastic turbulence is still poorly
understood.

Here, we shed light on the transition to elastic turbulence in parallel shear
flows of polymer solutions, like flow through a plane channel or a straight
pipe. For the parameters relevant to typical experiments with dilute polymer
solutions, parallel shear flows are linearly stable [26, 29]. It has
previously been proposed theoretically [30, 31] that such flows exhibit a
direct, sub-critical transition to elastic turbulence, similar to their
Newtonian counterparts. These predictions are supported by experiments in
straight microfluidic channels that report the existence of strong velocity
fluctuations that appear sub-critically above the onset flow rate [32, 33, 34,
35, 36]. The emerging transition scenario suggests strong parallels between
purely elastic and Newtonian flows [30], implicating coherent structures in
organising the phase-space dynamics of both systems. Yet, very little is known
about purely elastic flow structures experimentally [36, 37], and we are not
aware of any numerical simulations of such flows to date. Despite the
aforementioned advances, purely elastic turbulence is still being commonly
rationalised in terms of the dynamics of individual polymer molecules in
imposed random flows [27]. In this Letter, we demonstrate numerically the
existence of two-dimensional travelling wave solutions (TWS) in parallel shear
flows of dilute polymer solutions at vanishingly small Reynolds numbers and
discuss how they are related to turbulent dynamics in such flows.

We perform direct numerical simulations of two-dimensional pressure-driven
channel flow of a model simplified Phan-Thien-Tanner (PTT) polymeric fluid
[38], chosen to capture the shear-thinning nature of dilute polymer solutions
[26]. We performed linear stability analysis and confirmed that the flow is
linearly stable for all parameters studied here (see [39] for detail). We
consider a straight channel formed by two infinite parallel plates, with $x$
and $y$ being Cartesian coordinates along the streamwise and gradient
directions, respectively. Equations are rendered dimensionless by using $d$,
$\mathcal{U}$, $d/\mathcal{U}$, and $\eta_{p}\mathcal{U}/d$, and
$(\eta_{s}+\eta_{p})\mathcal{U}/d$ as the units of length, velocity, time,
stress, and pressure, respectively. Here, $d$ is the channel half width,
$\eta_{s}$ and $\eta_{p}$ are the solvent and polymeric contributions to the
viscosity, and $\mathcal{U}$ is the maximum value of the laminar fluid
velocity of a Newtonian fluid with the viscosity $\eta_{s}+\eta_{p}$ at the
same value of the applied pressure gradient. In two spatial dimensions, the
dimensionless equations of motion are given by

$\displaystyle\frac{\partial{\bm{c}}}{\partial
t}+{\bm{v}}\cdot\nabla{\bm{c}}-\left(\nabla{\bm{v}}\right)^{T}\cdot{\bm{c}}-{\bm{c}}\cdot\left(\nabla{\bm{v}}\right)=\kappa\nabla^{2}{\bm{c}}$
$\displaystyle\qquad\qquad\qquad\qquad-\frac{{\bm{c}}-\mathbb{I}}{W\\!i}\Bigg{[}1-2\,\epsilon+\epsilon\,\mathrm{Tr}{\bm{c}}\Bigg{]},$
(1) $\displaystyle\frac{\partial{\bm{v}}}{\partial
t}+{\bm{v}}\cdot\nabla{\bm{v}}=-\nabla p+\frac{\beta}{Re}\nabla^{2}{\bm{v}}$
$\displaystyle\qquad\qquad\qquad\qquad+\frac{(1-\beta)}{Re\,W\\!i}\nabla\cdot{\bm{c}}+\begin{pmatrix}2/Re\\\
0\end{pmatrix},$ (2)
$\displaystyle\qquad\qquad\qquad\qquad\nabla\cdot{\bm{v}}=0,$ (3)

where, $p$ is the pressure, $\bm{c}$ is the polymer conformation tensor, and
$\bm{v}$ is the fluid velocity.

The parameter space of this model fluid is spanned by four dimensionless
numbers [26]: i) the Weissenberg number, $W\\!i=\lambda\,\mathcal{U}/d$, that
controls how strongly polymer molecules are stretched in a flow, and, thus,
the magnitude of the elastic stresses; here, $\lambda$ is the polymer
relaxation time. ii) the Reynolds number,
$Re=\rho\,\mathcal{U}d/(\eta_{s}+\eta_{p})$, that gives the relative
importance of inertia compared to the viscous stresses; here, $\rho$ is the
density of the fluid. Elastic turbulence is characterised by small values of
the Reynolds number and throughout this work we set $Re=10^{-2}$. iii) The
ratio of the solvent to the total viscosity of the solution,
$\beta=\eta_{s}/(\eta_{s}+\eta_{p})$, that acts as an indirect measure of the
polymer concentration in dilute polymer solutions. iv) The strength of shear-
thinning, $\epsilon$, that controls how fast elastic stresses grow with
$W\\!i$. To represent the weakly shear-thinning nature of dilute polymer
solutions, we set $\epsilon=10^{-3}$. We note that the choice of $\mathcal{U}$
as the velocity scale implies that the maximum value of the laminar streamwise
velocity profile can exceed unity due to shear thinning. The fluid velocity
obeys the no-slip boundary condition, $\bm{v}(x,y=\pm 1,t)=0$. The boundary
conditions for the conformation tensor [40] are obtained by requiring that
${\bm{c}}(x,\pm 1,t)$ is equal to the values obtained at the boundaries from
Eq.(Coherent structures in plane channel flow of dilute polymer solutions with
vanishing inertia) with $\kappa=0$.

Simulations are performed with an in-house MPI-parallel code developed within
the Dedalus spectral framework [41] on a rectangular domain
$\left[0,L_{x}\right]\times\left[-1,1\right]$, where $L_{x}=10$; periodic-
boundary conditions are applied in the $x$-direction. All physical fields are
represented by a spectral decomposition [42] with $256$ Fourier and $1024$
Chebyshev modes. We checked that this resolution is sufficient to achieve
numerical convergence. Time-iteration uses a four-stage, third-order implicit-
explicit Runge-Kutta method [43] with the timestep $dt=5\cdot 10^{-3}$.

The conformation diffusion term, $\kappa\nabla^{2}{\bm{c}}$, with
$\kappa=5\cdot 10^{-5}$, is added to Eq.(Coherent structures in plane channel
flow of dilute polymer solutions with vanishing inertia) of the main text to
stabilise the numerics. Its presence is motivated by kinetic theories of
dilute polymer solutions, where it stems from a mean-field approximation of
Brownian diffusion of individual polymer molecules [44], and is required to
ensure the formal existence of solutions to polymeric equations of motion
[45]. In the past, this term was often employed with a value of $\kappa$ that
was orders of magnitude larger than the values predicted by kinetic theories
[46], thus calling into question the validity of their results. We stress that
the value of $\kappa$ used in this work does not suffer from this drawback. It
was selected to be sufficiently large to stabilise the dynamics against short-
wave-length numerical instabilities [47], yet small enough to be consistent
with the kinetic theory values. Indeed, re-casting the kinetic theory
prediction in terms of our dimensionless parameters yields
$\kappa=D\lambda/d^{2}W\\!i$, where $D$ is the diffusion coefficient of a
polymer molecule in equilibrium. Using $D\sim 1\mu$m2/s, $\lambda\sim 10$s,
and $d\sim 100\mu$m, typical for the microfluidic experiments of Pan _et al._
[34], and setting $W\\!i\sim 30$, corresponding to the onset of sub-critical
solutions reported below, yields $\kappa\sim 3\cdot 10^{-5}$. While this
estimate can yield smaller values of the dimensionless diffusivity for other
fluids and wider channels, we have confirmed that our value of $\kappa$ can be
further decreased by using higher resolutions and smaller values of $dt$
without altering the results. Similar values of $\kappa$ were previously shown
[48] to yield results indistinguishable from the ones obtained with
$\kappa=0$.

All simulations are started from an initial condition comprising the laminar
profile and a small localised perturbation [49] to the $xx$-component of the
conformation tensor in the following form:

$\displaystyle\Delta\exp{\left[-\frac{25}{8}\left\\{\left(\frac{2x}{L_{x}}-1\right)^{2}+y^{2}\right\\}\right]},$
(4)

where the amplitude $\Delta$ is fixed to be $1\%$ of the maximum value of
$c_{xx}$ in the laminar state. We have verified that other initial conditions,
including small random perturbations added to the laminar profile, lead to the
same ultimate steady state; however, we found that perturbations in the form
presented above are the most efficient in triggering TWS. Simulations were
stopped when the absolute value of the time derivative of the kinetic energy
stayed below $10^{-10}$ for more than $10$ time units.

|
---
|
---
Figure 1:  Exact travelling-wave solutions. The trace of the conformation
tensor (colour) and the flow streamlines (solid lines) for $\beta=0.8$ and
$W\\!i=26$ (a), $W\\!i=40$ (b), and $W\\!i=80$ (c). The mean flow is from left
to right. The streamlines represent velocity deviation from the mean
streamwise profile ${\overline{v}}_{x}(y)$. d, The speed of the travelling
wave solutions determined by tracking the position of the maximum of
$\mathrm{Tr}\,{\bm{c}}$ as a function of time. e, The relative kinetic energy
of the travelling wave solutions tracing the upper branch of the corresponding
bifurcation from infinity. f, The saddle-node values $W\\!i_{sn}$ (black
symbols) and their rheological counterparts $W\\!i_{sn}^{(rheo)}$ (red
symbols) determined as the lowest value of the Weissenberg number at which the
solution can be sustained.

Above a critical value of the Weissenberg number, we observe the appearance of
purely elastic flow structures that are different from the laminar flow. In
Figs.1a-c we present examples of such structures for $\beta=0.8$ and
$W\\!i=26$, $40$, and $80$. The distinctive features of these states are the
thin, filament-like arrangements of the polymer stretch, which is proportional
to the trace of the conformation tensor $\bf c$. While at lower values of
$W\\!i$, the maximum polymer extension occurs at the channel walls, as in the
laminar profile, at higher $W\\!i$, it is localised around the centre of the
channel. The arrangements of the stress and velocity profiles in these states
move downstream with a constant velocity, i.e. they are steady-states in a co-
moving frame and, thus, represent travelling-wave solutions to the equations
of motion. In Fig.1d, we present their downstream speed $c$, measured by
tracking in time the streamwise position of
$\mathrm{max(Tr}\,\bm{c}\mathrm{)}$. Similar to their Newtonian counterparts
[50], purely elastic TWS move slower than the laminar velocity profile.

Since the laminar flow is linearly stable for the parameters considered here,
TWS appear through a sub-critical bifurcation from infinity [51]. Our direct
numerical simulations resolve the upper branches of the bifurcation diagram
(Fig.1e) that we report in terms of the total kinetic energy of the flow
compared to its laminar value at the same $W\\!i$ (see [39] for details).
Lower branches, which determine the strength of a finite-amplitude
perturbation necessary to destabilise the laminar flow [30], are usually
linearly unstable and cannot be determined by direct numerical simulations.
The bifurcation diagrams terminate at the saddle-node value of the Weissenberg
number, $W\\!i_{sn}$, determined as the lowest value of $W\\!i$ for which TWS
were observed. In Fig.1f we map the region of existence of TWS for various
values of $\beta$.

|  |
---|---|---
Figure 2:  Characterisation of the TWS velocity field for $\beta=0.8$ and
$W\\!i=80$. a, The velocity deviation from the mean profile
$(v_{x}(x,y)-{\overline{v}}_{x}(y),v_{y}(x,y))$ (vectors), and the out-of-
plane component of the vorticity
$\partial_{x}v_{y}(x,y)-\partial_{y}(v_{x}(x,y)-{\overline{v}}_{x}(y))$
(colour). b, Comparison between the laminar and the mean streamwise velocity
profiles.

To facilitate comparison with experiments, we express $W\\!i_{sn}$ in terms of
the rheological Weissenberg number $W\\!i_{sn}^{(rheo)}$ often used to report
the results of experiments dealing with shear-thinning fluids (see, for
instance, [52, 34, 53]). Since the relaxation time and the viscosity of such
fluids change with the local shear rate, the nominal value of $W\\!i$, which
is based on the relaxation time $\lambda$ measured at a single (low) shear
rate, is an overestimate of the strength of the elastic stresses. Instead, the
rheological Weissenberg number provides an estimate of the value of the
Weissenberg number in an Oldroyd-B fluid that would generate elastic stresses
of the same magnitude; here, we employ the definition as used by Pan _et al._
[34] (see also [39]). In Fig.1f, we report the onset of TWS in terms of the
rheological Weissenberg number $W\\!i_{sn}^{(rheo)}$. In the experiments of
Pan _et al._ [34], a sudden onset of large velocity fluctuations was observed
for $W\\!i^{(rheo)}>5.4$; at those flow rates, the total viscosity of the
solution was measured to be $\eta_{s}+\eta_{p}\sim 0.35$Pa$\cdot$s, while the
solvent viscosity $\eta_{s}\sim 0.2$Pa$\cdot$s, yielding $\beta\sim 0.57$ (see
Supplemental Material for [34]). As can be seen from Figs.1f and 1e, the sub-
critical nature of the transition, the value of $W\\!i_{sn}$, and the
amplitude of the jump of the velocity fluctuations around $W\\!i_{sn}$
reported by Pan _et al._ [34], are consistent with the onset of TWS found
here.

The exact solutions reported in our work resemble the ‘arrowhead’ structures
[54, 55, 48, 56, 57] found in elasto-inertial ($Re\gg 1$ and $W\\!i\gg 1$)
pressure-driven channel flows. There, they originate from a linear instability
recently discovered for $1-\beta\ll 1$, and extend sub-critically to lower
values of $Re$ and $W\\!i$ [48, 57]. This linear instability can be
continuously traced to the purely elastic regime [58, 57], where, similar to
its elasto-inertial analogue, it exists only in a narrow range $1-\beta\ll 1$.
While such a linear instability is probably not directly relevant for
experiments with dilute polymer solutions, as it is difficult to
simultaneously achieve $Re\ll 1$, $W\\!i\gg 1$, and $1-\beta\ll 1$, we
speculate that the non-linear state that originates from it directly connects
to TWS discovered here. Importantly, our results demonstrate that no linear
instability is required to sustain travelling waves in dilute polymer
solutions at experimentally relevant values of $\beta$, and that their origin
is purely elastic in nature.

Several prominent features of TWS are good potential candidates for
experimental detection. The velocity field associated with TWS is dominated by
a solitary pair of counter-rotating vortices arranged symmetrically around the
centreline of the channel, see Fig.2a. Such structures resemble ‘diwhirls’
previously reported in Taylor-Couette flows of dilute polymer solutions [59,
60]. In the absence of other prominent velocity structures, they might be
directly observable in particle image velocimetry, although their magnitude is
quite low when compared to the mean profile. The latter, defined as
${\overline{v}}_{x}(y)=\int_{0}^{L_{x}}v_{x}(x,y)dx/L_{x}$, appears to be the
most promising for experimental observation feature of the velocity field, and
is the main observable used by Pan et al. [34]. However, its deviations from
the laminar profile are localised around the centreline and are weak (Fig.2b).
This is in contrast with Newtonian parallel shear flows, where turbulent
production is largely associated with the walls [14, 50]. The associated
stresses, on the other hand, are about $10\%$ of their laminar values, have a
distinctive, filament-like spatial distribution, and should be detectable in
birefringence measurements. Finally, a more pronounced, though less specific,
signature of TWS is the pressure drop along the channel that can become
comparable to its laminar value (see Fig.S2 in [39]).

Two-dimensional TWS presented here are steady-states in a co-moving frame and
do not lead to time-dependent fluctuations associated with elastic turbulence.
To be relevant to the latter, they are expected to lose their stability for
sufficiently high $W\\!i$. While we have verified numerically that TWS are
stable with respect to small two-dimensional perturbations, their stability
with respect to three-dimensional ones is yet unknown. There are two possible
scenarios. First, the upper branch solutions in Fig.1e can be linearly stable
up to $W\\!i_{3D}$, which would mark the appearance of three-dimensional
coherent structures. In this scenario, TWS should be directly observable for
$W\\!i_{sn}<W\\!i<W\\!i_{3D}$, and their properties can be tested against our
predictions reported in Fig.1. Alternatively, the whole upper branch can be
unstable with respect to three-dimensional perturbations. This scenario would
mirror the transition to Newtonian turbulence in plane channel flows. There,
two-dimensional Tollmien-Schlichting TWS originate from a linear instability
at $Re=5772$ and extend sub-critically to lower values of $Re$ [61], while the
corresponding upper branch is unstable towards three-dimensional perturbations
for any value of $Re$ [62]. The close agreement between the onset of velocity
fluctuations reported by Pan _et al._ [34] and the onset of two-dimensional
TWS discovered here suggests that $W\\!i_{3D}$ is either close or equal to
$W\\!i_{sn}$. In this case, the flow properties would be determined by three-
dimensional structures that are yet to be discovered. Recent studies of
elasto-inertial turbulence suggest that such structures are only weakly three-
dimensional [17, 54, 55, 22, 24, 56], and we could expect the two-dimensional
flow features discussed above to be observable in three-dimensional channel
flows.

Finally, we note that a series of recent experiments from the Steinberg’s
group [37, 63, 64, 65] has suggested a potentially different transition
scenario that is rationalised in terms of _elastic Alfvén waves_[66]. While
our results are consistent with our earlier proposal [30, 31], and with the
experiments of the Arratia’s group [34, 35, 36], further experimental and
numerical work is required to understand the transition to elastic turbulence
in parallel shear flows. We note that some of the discrepancies between the
two sets of observations can be resolved if the aformentioned elastic waves
are interpreted as TWS discussed here.

## Acknowledgements

This work was partially funded by EPSRC (grant number EP/I004262/1). We would
like to thank EPSRC for the computational time made available on the ARCHER2
UK National Supercomputing Service (https://www.archer2.ac.uk) via the UK
Turbulence Consortium (EP/R029326/1). We thank Paulo Arratia, Martin Lellep,
Moritz Linkmann, Davide Marenduzzo, Jacob Page, and Rob Poole for helpful
discussions. We gratefully acknowledge the support we received from the
Dedalus team (https://dedalus-project.org), and Keaton Burns in particular.

## References

  * Ditto  _et al._ [1990] W. L. Ditto, S. N. Rauseo, and M. L. Spano, Phys. Rev. Lett. 65, 3211 (1990).
  * Hunt [1991] E. R. Hunt, Phys. Rev. Lett. 67, 1953 (1991).
  * Petrov  _et al._ [1993] V. Petrov, V. Gáspár, J. Masere, and K. Showalter, Nature 361, 240–243 (1993).
  * Tuckerman  _et al._ [2020] L. S. Tuckerman, M. Chantry, and D. Barkley, Annu. Rev. Fluid Mech. 52, 343 (2020).
  * Linkmann and Morozov [2015] M. F. Linkmann and A. Morozov, Phys. Rev. Lett. 115, 134502 (2015).
  * Garfinkel  _et al._ [1992] A. Garfinkel, M. L. Spano, W. L. Ditto, and J. N. Weiss, Science 257, 1230 (1992).
  * Roy  _et al._ [1992] R. Roy, T. W. Murphy, T. D. Maier, Z. Gills, and E. R. Hunt, Phys. Rev. Lett. 68, 1259 (1992).
  * Auerbach  _et al._ [1987] D. Auerbach, P. Cvitanović, J.-P. Eckmann, G. Gunaratne, and I. Procaccia, Phys. Rev. Lett. 58, 2387 (1987).
  * Waleffe [1997] F. Waleffe, Phys. Fluids 9, 883 (1997).
  * Hof  _et al._ [2006] B. Hof, J. Westerweel, T. M. Schneider, and B. Eckhardt, Nature 443, 59 (2006).
  * Avila  _et al._ [2011] K. Avila, D. Moxey, A. de Lozar, M. Avila, D. Barkley, and B. Hof, Science 333, 192 (2011).
  * Sano and Tamai [2016] M. Sano and K. Tamai, Nat. Phys. 12, 249 (2016).
  * Lemoult  _et al._ [2016] G. Lemoult, L. Shi, K. Avila, S. V. Jalikop, M. Avila, and B. Hof, Nat. Phys. 12, 254 (2016).
  * Barkley [2016] D. Barkley, J. Fluid Mech. 803, P1 (2016).
  * Stone  _et al._ [2002] P. A. Stone, F. Waleffe, and M. D. Graham, Phys. Rev. Lett. 89, 208301 (2002).
  * Dubief  _et al._ [2004] Y. Dubief, C. M. White, V. E. Terrapon, E. S. G. Shaqfeh, P. Moin, and S. K. Lele, J. Fluid Mech. 514, 271–280 (2004).
  * Samanta  _et al._ [2013] D. Samanta, Y. Dubief, M. Holzner, C. Schäfer, A. N. Morozov, C. Wagner, and B. Hof, Proc. Natl. Acad. Sci. U.S.A. 110, 10557–10562 (2013).
  * Choueiri  _et al._ [2018] G. H. Choueiri, J. M. Lopez, and B. Hof, Phys. Rev. Lett. 120, 124501 (2018).
  * Lopez  _et al._ [2019] J. M. Lopez, G. H. Choueiri, and B. Hof, J. Fluid Mech. 874, 699–719 (2019).
  * Yamani  _et al._ [2021] S. Yamani, B. Keshavarz, Y. Raj, T. A. Zaki, G. H. McKinley, and I. Bischofberger, Phys. Rev. Lett. 127, 074501 (2021).
  * Choueiri  _et al._ [2021] G. H. Choueiri, J. M. Lopez, A. Varshney, S. Sankar, and B. Hof, Proc. Natl. Acad. Sci. U.S.A. 118 (2021).
  * Shekar  _et al._ [2019] A. Shekar, R. M. McMullen, S.-N. Wang, B. J. McKeon, and M. D. Graham, Phys. Rev. Lett. 122, 124503 (2019).
  * Shekar  _et al._ [2020] A. Shekar, R. M. McMullen, B. J. McKeon, and M. D. Graham, J. Fluid Mech. 897, A3 (2020).
  * Shekar  _et al._ [2021] A. Shekar, R. M. McMullen, B. J. McKeon, and M. D. Graham, Phys. Rev. Fluids 6, 093301 (2021).
  * Groisman and Steinberg [2000] A. Groisman and V. Steinberg, Nature 405, 53 (2000).
  * Datta  _et al._ [2021] S. Datta, A. Ardekani, P. Arratia, A. Beris, I. Bischofberger, G. McKinley, J. Eggers, E. Lopez-Aguilar, S. Fielding, A. Frishman, M. Graham, J. Guasto, S. Haward, A. Shen, S. Hormozi, R. Poole, A. Morozov, V. Shankar, E. Shaqfeh, H. Stark, V. Steinberg, G. Subramanian, and H. Stone, arXiv preprint arXiv:2108.09841 (2021).
  * Steinberg [2021] V. Steinberg, Annu. Rev. Fluid Mech. 53, 27 (2021).
  * Denn [2001] M. M. Denn, Annu. Rev. Fluid Mech. 33, 265 (2001).
  * Castillo Sánchez  _et al._ [2022] H. A. Castillo Sánchez, M. R. Jovanović, S. Kumar, A. Morozov, V. Shankar, G. Subramanian, and H. J. Wilson, J. Non-Newtonian Fluid Mech. 302, 104742 (2022).
  * Morozov and van Saarloos [2007] A. N. Morozov and W. van Saarloos, Phys. Rep. 447, 112 (2007).
  * Morozov and van Saarloos [2019] A. Morozov and W. van Saarloos, J. Stat. Phys. 175, 554 (2019).
  * Bertola  _et al._ [2003] V. Bertola, B. Meulenbroek, C. Wagner, C. Storm, A. N. Morozov, W. van Saarloos, and D. Bonn, Phys. Rev. Lett. 90, 114502 (2003).
  * Bonn  _et al._ [2011] D. Bonn, F. Ingremeau, Y. Amarouchene, and H. Kellay, Phys. Rev. E 84, 10.1103/PhysRevE.84.045301 (2011).
  * Pan  _et al._ [2013] L. Pan, A. Morozov, C. Wagner, and P. E. Arratia, Phys. Rev. Lett. 110, 174502 (2013).
  * Qin and Arratia [2017] B. Qin and P. E. Arratia, Phys. Rev. Fluids 2, 083302 (2017).
  * Qin  _et al._ [2019] B. Qin, P. F. Salipante, S. D. Hudson, and P. E. Arratia, Phys. Rev. Lett. 123, 194501 (2019).
  * Jha and Steinberg [2020] N. K. Jha and V. Steinberg, arXiv:2009.12258 (2020).
  * Phan-Thien and Tanner [1977] N. Phan-Thien and R. I. Tanner, J. Non-Newtonian Fluid Mech. 2, 353 (1977).
  * [39] See Supplemental Material at [URL will be inserted by publisher], which includes Refs.[67, 68].
  * Liu and Khomami [2013] N. Liu and B. Khomami, J. Fluid Mech. 737, R4 (2013).
  * Burns  _et al._ [2020] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, Phys. Rev. Res. 2, 023068 (2020).
  * Canuto  _et al._ [2012] C. Canuto, M. Y. Hussaini, A. Quarteroni, A. Thomas Jr, _et al._ ,  _Spectral methods in fluid dynamics_ (Springer Science & Business Media, 2012).
  * Ascher  _et al._ [1997] U. M. Ascher, S. J. Ruuth, and R. J. Spiteri, Appl. Numer. Math. 25, 151 (1997).
  * Beris and Mavrantzas [1994] A. N. Beris and V. G. Mavrantzas, J. Rheol. 38, 1235 (1994).
  * Ardith W. El-Kareh and Leal [1989] A. W. Ardith W. El-Kareh and L. G. Leal, J. Non-Newtonian Fluid Mech. 33, 257 (1989).
  * Sureshkumar and Beris [1995] R. Sureshkumar and A. N. Beris, J. Non-Newtonian Fluid Mech. 60, 53 (1995).
  * Owens and Phillips [2002] R. G. Owens and T. N. Phillips,  _Computational rheology_ (World Scientific, 2002).
  * Page  _et al._ [2020] J. Page, Y. Dubief, and R. R. Kerswell, Phys. Rev. Lett. 125, 154501 (2020).
  * Schumacher and Eckhardt [2001] J. Schumacher and B. Eckhardt, Phys. Rev. E 63, 046307 (2001).
  * Graham and Floryan [2021] M. D. Graham and D. Floryan, Annu. Rev. Fluid Mech. 53, 227 (2021).
  * Rosenblat and Davis [1979] S. Rosenblat and S. H. Davis, SIAM J. Appl. Math. 37, 1 (1979).
  * McKinley  _et al._ [1996] G. H. McKinley, P. Pakdel, and A. Oeztekin, J. Non-Newton. Fluid Mech. 67, 19 (1996).
  * Schäfer  _et al._ [2018] C. Schäfer, A. Morozov, and C. Wagner, J. Non-Newtonian Fluid Mech. 259, 78 (2018).
  * Dubief  _et al._ [2013] Y. Dubief, V. E. Terrapon, and J. Soria, Phys. Fluids 25, 110817 (2013).
  * Sid  _et al._ [2018] S. Sid, V. E. Terrapon, and Y. Dubief, Phys. Rev. Fluids 3, 011301 (2018).
  * Dubief  _et al._ [2021] Y. Dubief, J. Page, R. R. Kerswell, V. E. Terrapon, and V. Steinberg, arXiv:2006.06770 (2021).
  * Buza  _et al._ [2021] G. Buza, J. Page, and R. R. Kerswell, arXiv:2107.06191 (2021).
  * Khalid  _et al._ [2021] M. Khalid, V. Shankar, and G. Subramanian, Phys. Rev. Lett. 127, 134502 (2021).
  * Groisman and Steinberg [1997] A. Groisman and V. Steinberg, Phys. Rev. Lett. 78, 1460 (1997).
  * Kumar and Graham [2000] K. A. Kumar and M. D. Graham, Phys. Rev. Lett. 85, 4056 (2000).
  * Herbert [1976] T. Herbert, in  _Proceedings of the Fifth International Conference on Numerical Methods in Fluid Dynamics, June 28-July 3, 1976, Twente University of Technology_ , edited by A. I. van de Vooren and P. J. Zandbergen (Springer-Verlag Berlin, 1976).
  * Orszag and Patera [1983] S. A. Orszag and A. T. Patera, J. Fluid Mech. 128, 347–385 (1983).
  * Jha and Steinberg [2021] N. K. Jha and V. Steinberg, Proc. Natl. Acad. Sci. U.S.A. 118 (2021).
  * Shnapp and Steinberg [2021] R. Shnapp and V. Steinberg, arXiv:2106.01817 (2021).
  * Li and Steinberg [2022] Y. Li and V. Steinberg, Elastic instability in a straight channel of viscoelastic flow without prearranged perturbations (2022).
  * Varshney and Steinberg [2019] A. Varshney and V. Steinberg, Nat. Commun. 10, 652 (2019).
  * Zhang  _et al._ [2013] M. Zhang, I. Lashgari, T. A. Zaki, and L. Brandt, J. Fluid Mech. 737, 249–279 (2013).
  * Poole  _et al._ [2019] R. J. Poole, M. Davoodi, and K. Zografos, Bull. - Br. Soc. Rheol. 60, 29 (2019).