# Efficient CNN Building Blocks for Encrypted Data

Nayna Jain1,4, Karthik Nandakumar2, Nalini Ratha3, Sharath Pankanti5, Uttam
Kumar 1

###### Abstract

Machine learning on encrypted data can address the concerns related to privacy
and legality of sharing sensitive data with untrustworthy service providers,
while leveraging their resources to facilitate extraction of valuable insights
from otherwise non-shareable data. Fully Homomorphic Encryption (FHE) is a
promising technique to enable machine learning and inferencing while providing
strict guarantees against information leakage. Since deep convolutional neural
networks (CNNs) have become the machine learning tool of choice in several
applications, several attempts have been made to harness CNNs to extract
insights from encrypted data. However, existing works focus only on ensuring
data security and ignore security of model parameters. They also report high
level implementations without providing rigorous analysis of the accuracy,
security, and speed trade-offs involved in the FHE implementation of generic
primitive operators of a CNN such as convolution, non-linear activation, and
pooling. In this work, we consider a Machine Learning as a Service (MLaaS)
scenario where both input data and model parameters are secured using FHE.
Using the CKKS scheme available in the open-source HElib library, we show that
operational parameters of the chosen FHE scheme such as the degree of the
cyclotomic polynomial, depth limitations of the underlying leveled HE scheme,
and the computational precision parameters have a major impact on the design
of the machine learning model (especially, the choice of the activation
function and pooling method). Our empirical study shows that choice of
aforementioned design parameters result in significant trade-offs between
accuracy, security level, and computational time. Encrypted inference
experiments on the MNIST dataset indicate that other design choices such as
ciphertext packing strategy and parallelization using multithreading are also
critical in determining the throughput and latency of the inference process.

## Introduction

Deep neural networks have proven to be a promising technology due to their
ability to achieve competitive machine learning performance in diverse domains
including computer vision. Ever since AlexNet achieved a top-5 error rate of
16.4% using Convolutional Neural Networks (CNN) for the image classification
task on the ImageNet dataset of 1.2 million high-resolution images as part of
ILVSRC 2012 (Krizhevsky, Sutskever, and Hinton 2012) competition, various CNN
architectures (Simonyan and Zisserman 2014)(Szegedy et al. 2015)(Zeiler and
Fergus 2013) have been developed to improve accuracy and other performance
metrics. Though the continuous improvement in the performance of CNNs has
hastened their adoption in broader computer vision applications such as mobile
and embedded devices (Howard et al. 2017)(Phan et al. 2020), the large
computing/memory resource requirements of CNNs present practical challenges in
applying them in many applications. This is often overcome by outsourcing the
network training and/or inference computations to the Cloud, which is commonly
referred to as Machine Learning as a Service (MLaaS).

Figure 1: In a conventional Machine Learning as a Service (MLaaS) scenario,
both the data and model parameters are unencrypted. Existing works on secure
inference assume only the data is encrypted. Our proposed approach encrypts
both the data and model parameters.

In MLaaS, there are three parties: (i) the end-user, who requires predictions
on sensitive data, (ii) the cloud service provider, who has access to large
computational resources, and (iii) the model owner, who has developed the
inference model, trained using machine learning algorithms that may be
proprietary. In conventional MLaaS scenarios, the cloud service providers
typically have access to the inference models as well as the content of the
query data. In many applications, such access to data and the model is
undesirable because the query data/models could reveal sensitive information.
Moreover, sharing of such data may be prohibited under emerging privacy
regulations such as GDPR. This has led to the development of a variety of
privacy preserving machine learning algorithms. In this paper, we focus on a
specific class of privacy preserving machine learning algorithms that rely on
Homomorphic Encryption (HE).

Most of the approaches presented in the literature for encrypted inference
encrypt only the input data and do not encrypt the model parameters. This is
based on the assumption that the same entity owns the model and runs the cloud
service. However, access to the model parameters could lead to competitive
disadvantage for the model owner, when the model owner and cloud service
provider are separate entities (e.g., a start-up providing the MLaaS service
on a Amazon or Google cloud). Hence, it becomes imperative to protect both the
input data and model parameters in such scenarios (see Figure 1). (Bost et al.
2015) discussed about encrypted models in the context of Hyperplane Decision,
Navie Bayes and Decision Tree algorithm. In this work, we focus on the
scenario where the model owner and cloud service provider are separate
entities and the cloud provider is not trusted. Therefore, we encrypt the
model parameters111Note that the model architecture is still available to the
cloud service provider in the clear. in addition to the input data using the
end-user’s public key (as shown on Figure 1). The advantage of this approach
is that the cloud service provider cannot derive any benefit from the model
directly. The disadvantage of model parameter re-encryption using the specific
client’s public key is the increased computation and communication burden of
the model owner.

### Contributions

The objective of this work is to develop the component building blocks for
enabling generic CNN inference on encrypted data. We use the implementation of
the CKKS scheme available in the HElib library (Halevi and Shoup 2020) for all
our analysis. However, we note that the lack of a bootstrapping implementation
for CKKS makes it impossible to perform inferencing on any arbitrary CNN.
Hence, we design depth-constrained CNNs that operate within the multiplicative
depth limitations of the underlying leveled HE scheme. Our contributions are
four-fold:

  1. 1.

MLaaS scenario for convolution neural networks is proposed where both the
input data and model parameters are encrypted, thereby enabling a model owner
to provide the MLaaS service on an untrusted cloud.

  2. 2.

CNN inference on encrypted data using the CKKS scheme is presented, which
eliminates the need for careful model parameter quantization.

  3. 3.

Implementation and performance analysis of illustrative linear operators such
as convolution and matrix multiplication as well as non-linear operations such
as ReLU and maxpooling in the encrypted domain is provided.

  4. 4.

Various multithreading strategies were explored and the analysis revealed that
the optimal strategy involves a complex interplay between system resources and
the computations involved.

## Background and Related Work

### Homomorphic Encryption

Homomorphic encryption schemes enable a class of computations on the
ciphertexts (i.e., encrypted data) without decrypting them (Acar et al. 2018).
Let $\llbracket x\rrbracket$ denote the encryption of value $x$ using a public
key $pk$. A cryptosystem is said to be _fully homomorphic_ (FHE) (Gentry 2009)
if it enables computation of any arbitrary function
$f(x_{1},x_{2},\cdots,x_{d})$ in the encrypted domain. In other words, there
exists a function $g$ such that
$f(x_{1},x_{2},\cdots,x_{d})=\mathcal{D}\left(g(\llbracket
x_{1}\rrbracket,\llbracket x_{2}\rrbracket,\cdots\llbracket
x_{d}\rrbracket),sk\right)$, where $\mathcal{D}$ represents the corresponding
decryption mechanism of the cryptosystem using private/secret key $sk$.

Recently, a number of FHE schemes have been proposed based on the hardness of
the Ring Learning With Errors (RWLE) problem. The most well-known examples of
such schemes are BFV (Brakerski 2012)(Fan and Vercauteren 2012), BGV
(Brakerski, Gentry, and Vaikuntanathan 2014), TFHE (Chillotti et al. 2020),
and CKKS (Cheon et al. 2017) cryptosystems, all of which support both additive
and multiplicative homomorphism. Among the well-known FHE schemes, only CKKS
natively supports floating-point operations, albeit with approximate (limited
precision) computations. Most of the previous attempts at encrypted inference
coarsely quantize the weights of a neural network in order to satisfy the
requirements of the underlying HE schemes (Gilad-Bachrach et al. 2016; Lou and
Jiang 2019). In this paper, we explore the use of the CKKS scheme, which
eliminates the need for transformations of operands as integers.

FHE is typically obtained by combining a leveled HE scheme with a
bootstrapping operation. Since each HE computation increases the “noise” level
of a ciphertext, a leveled HE scheme has a (parameter-dependent) limit on the
computational (multiplicative) depth of the arithmetic circuit that can be
computed, while still allowing for meaningful decryption. This limitation is
overcome by using Gentry’s bootstrapping technique (Gentry 2009), which
“refreshes” the ciphertext and reduces its noise level (at the cost of relying
on circular security). While bootstrapping can enable computation of
arithmetic circuits with arbitrary multiplicative depth, it comes at a
prohibitively high computational cost. For some FHE constructs such as CKKS,
there is no existing implementation of true “Gentry-style” bootstrapping. This
leads to an intractable conundrum. While modern machine (deep) learning
algorithms rely on their depth to achieve good generalization performance
(Poggio et al. 2017), deeper networks are inherently more difficult to compute
in the encrypted domain without bootstrapping. Therefore, optimizing the
computational depth to achieve an acceptable trade-off between accuracy and
practical computational feasibility is one of the fundamental challenges in
encrypted deep learning.

### Secure Machine Learning Using FHE

Every supervised machine learning algorithm involves two steps: (i) learning
the model from training data; (ii) inference on new unknown samples. We
presume that an efficient model has already been trained from data available
in clear. Our focus is on how to carry out encrypted inferencing with minimal
leakage of information for all the parties in MLaaS while providing guaranteed
privacy.

Privacy preserving machine learning can have different challenges and threats
based on what is required to be private. In (Al-Rubaie and Chang 2019),
various types of attacks (_e.g._ , reconstruction attacks, model inversion
attacks, membership inference attacks) that could potentially reveal
private/proprietary information from data/model have been discussed. The work
also presents various cryptographic and non-cryptographic approaches to thwart
these attacks. Key techniques that have been explored to achieve privacy
preserving machine learning include differential privacy, secure multiparty
computation (SMPC), and homomorphic encryption. Homomorphic encryption has
been considered more suitable for the cloud-based applications, but its
practical feasibility is still in question. SMPC-based methods have been
generally deemed more practical, but it requires interaction between the
client and server resulting in large communication cost. Differential privacy
techniques lack the formal guarantees of security under assumptions of more
conservative attack models.

Most of the recent work in privacy-preserving machine learning (Gilad-Bachrach
et al. 2016)(Juvekar, Vaikuntanathan, and Chandrakasan 2018) (Chou et al.
2018)(Dathathri et al. 2019) has primarily focused on secure inference where
the model has already been trained in the clear. Once the model has been
trained, these methods use either SMPC or homomorphic encryption techniques to
enable inference. While CryptoNets (Gilad-Bachrach et al. 2016) rely only on
homomorphic encryption to enable inferencing over private input data, Gazelle
(Juvekar, Vaikuntanathan, and Chandrakasan 2018) uses a combination of
homomorphic encryption and garbled circuits to achieve two orders of magnitude
speed improvement. Faster CryptoNets (Chou et al. 2018) have demonstrated
improved performance over that of CryptoNets by leveraging sparsity properties
and reduced the amortized inference time from minutes to seconds. CHET
(Dathathri et al. 2019) provides an optimizing compiler for Fully Homomorphic
Neural Network Inference. It basically provides high level user framework to
automate the tuning of parameters for security as well as performance without
degrading the accuracy. The framework provides high level and low level
intermediate representations that can be linked with different homomorphic
libraries. Their results showed that automated homomorphic circuits from the
compiler outperform the circuits that were manually tuned. An encrypted
inference service built on top of TensorFlow using secure multiparty
computation is discussed in (Dahl et al. 2018). Encrypted model was considered
in (Bost et al. 2015) for Hyperplane Decision, Naive Bayes, Decision Tree
algorithms, but not convolution neural network.

### Convolutional Neural Network (CNN)

A CNN is a multi-layered neural network that is usually applied to image/video
data. Let $\mbox{\boldmath$X$}(\ell)$ and $\mbox{\boldmath$Z$}(\ell)$ denote
the input and output of layer $\ell$ ($\ell=1,2,\cdots,D$) in the network.
Each layer in the network performs specific arithmetic operations on the
outputs of the previous layer and passes the result to the subsequent layer,
i.e., $\mbox{\boldmath$X$}(\ell)=\mbox{\boldmath$Z$}(\ell-1)$ and
$\mbox{\boldmath$Z$}(0)$ is the raw input data. This structure allows higher-
level abstract features to be computed as non-linear functions of lower-level
features (starting with the raw data). The first few layers of a CNN generally
consists of the three types of operations.

  1. 1.

Convolution: This is a linear operation, where the input to the layer222For
the sake of convenience, we drop the layer index $\ell$ in the subsequent
discussion. ($X$) is convolved with a kernel ($W$) to generate a filtered
output $Y$. Let $\mbox{\boldmath$A$}_{i,j}$ denote the element in the $i$-th
row and $j$-th column of a matrix $A$. Suppose that the size of input $X$ is
$(M\times N)$ and the size of kernel $W$ is $(P\times Q)$. Typically, one can
pad the input with zeros to compute the filter responses at the edges. Let
$\widetilde{X}$ be the zero-padded version of the input $X$.

$\mbox{\boldmath$Y$}_{m,n}=\langle\mbox{\boldmath$\widehat{X}$}_{m,n},\mbox{\boldmath$\widehat{W}$}\rangle,$
(1)

where $\langle\mbox{\boldmath$a$},\mbox{\boldmath$b$}\rangle$ is the inner
product of two vectors $a$ and $b$, $\widehat{W}$ is the vectorized
(flattened) version of the kernel $W$, and
$\mbox{\boldmath$\widehat{X}$}_{m,n}$ is the vectorized version of the input
window (a window of size $P\times Q$ with
$\mbox{\boldmath$\widetilde{X}$}_{m,n}$ as the top-left element) selected from
the padded image. For the sake of simplicity, the stride length is assumed to
be $1$ in equation (1).

  2. 2.

Activation: In this step, a point-wise non-linear function is applied to the
filter responses. The commonly used activation functions are rectified linear
unit (ReLU), sigmoid, and hyperbolic tangent (tanh). In this paper, we
consider only the ReLU activation function, which is defined as follows:

$ReLU(a)=max(0,a)=\begin{cases}0,&\mbox{if }a\leq 0\\\ a,&\mbox{if
}a>0.\end{cases}$ (2)

  3. 3.

Pooling: This is typically used for dimensionality reduction, where multiple
responses within a neighborhood are pooled together. While max pooling is
often used for selecting the most dominant response in each neighborhood and
is a non-linear operation, mean pooling computes the average response in each
neighborhood and is a linear operation. The final output ($Z$) of the
activation and pooling layers can be modeled as a function $f$ of the filter
responses, i.e., $\mbox{\boldmath$Z$}=f(\mbox{\boldmath$Y$})$, where $f$ is
non-linear.

After the initial layers, the data is flattened into a vector and a few fully
connected (FC) layers are added to the network. The output of each node in the
FC layer is often computed by applying a non-linear activation function to the
weighted average of its inputs, which includes a bias term that always emits
value 1. This can be mathematically represented as
$\mbox{\boldmath$Z$}=f(\mbox{\boldmath$W$}\mbox{\boldmath$X$})$, where $W$ is
the weight matrix of the FC layer and $f$ is the activation function. The
final layer in a CNN is usually a softmax layer that provides an approximate
probability distribution over the class labels.

## Proposed Approach

In this work, we employ a typical CNN architecture as an illustrative example.
Though FHE schemes allow arbitrary computations on the encrypted data, two key
challenges need to be circumvented before applying a CNN on encrypted inputs.
The first challenge is the appropriate packing of input data to efficiently
make use of the Single Instruction Multiple Data (SIMD) operations. The next
challenge stems from the inherent inability of most FHE schemes to directly
compute non-linear functions (with the exception of polynomials). Often, the
non-linear function is approximated by an iterative algorithm or a polynomial
function. Since the approximation error is generally inversely proportional to
the computational depth (using more iterations or higher-degree polynomials
leads to less approximation error), the trade-off between computational depth
and accuracy must be carefully managed.

| (a)
---|---
(b) |
Figure 2: Ciphertext packing for batched inference. (a) Ciphertext
$\llbracket\mbox{\boldmath$X$}\rrbracket_{m,n}$ contains the encrypted values
of pixel $\llbracket x_{m,n}^{(\cdot)}\rrbracket$ from all the $K$ images in a
batch, where $K$ is the number of slots in the ciphertext. (b) Ciphertext
$\llbracket\mbox{\boldmath$W$}\rrbracket_{p,q}$ contains the encrypted value
of kernel weight $\llbracket w_{p,q}\rrbracket$ in all the $K$ slots.

#### Convolution Operator:

We employ the CKKS scheme for encrypted operations. The native plaintext in
the CKKS scheme is a polynomial in the cyclotomic ring, which enables packing
multiple plaintext values into different “slots” in the ciphertext. This
ciphertext packing enables parallelization of addition and multiplication
operations through SIMD operations. However, it must be emphasized that it is
not possible to randomly access values in the individual slots of the
ciphertext after packing. Since only a limited set of operations (e.g.,
rotation of values in the slots) are possible within a ciphertext, the
benefits of SIMD operations can be fully realized only when there is minimal
interaction between the slots. Therefore, we follow a more conservative and
straightforward approach to ciphertext packing that amortizes the
computational time by processing a batch of images in parallel.

Let $K$ be the number of slots available in the ciphertext. Given a batch of
$K$ images
$\mbox{\boldmath$X$}^{(1)},\mbox{\boldmath$X$}^{(2)},\cdots,\mbox{\boldmath$X$}^{(K)}$,
where each image is of size $M\times N$, we represent the encrypted image as a
ciphertext matrix $\llbracket\mbox{\boldmath$X$}\rrbracket$ of size $M\times
N$ as shown in Figure 2(a). The ciphertext
$\llbracket\mbox{\boldmath$X$}\rrbracket_{m,n}$ contains the encrypted values
of pixel $(m,n)$ from all the $K$ images, i.e.
$\llbracket\mbox{\boldmath$X$}\rrbracket_{m,n}=\left[\llbracket
x_{m,n}^{(1)}\rrbracket,\llbracket x_{m,n}^{(2)}\rrbracket,\cdots,\llbracket
x_{m,n}^{(K)}\rrbracket\right]$. Similarly, the encrypted convolution kernel
is also represented as a ciphertext matrix
$\llbracket\mbox{\boldmath$W$}\rrbracket$ of size $P\times Q$, where $P\times
Q$ is the kernel size333Encryption of the convolution kernel addresses the
more general scenario where the model provider may be different from the cloud
service provider.. However, as shown in Figure 2(b), the encrypted weight
value $\llbracket w_{p,q}\rrbracket$ is repeated across all the slots of the
ciphertext $\llbracket\mbox{\boldmath$W$}\rrbracket_{p,q}$.

While the above simple packing approach greatly increases the memory
requirement, it provides great simplicity in terms of implementing the
convolution operator. One can directly treat the encrypted image matrix
$\llbracket\mbox{\boldmath$X$}\rrbracket$ and the encrypted kernel
$\llbracket\mbox{\boldmath$W$}\rrbracket$ as equivalent of the plaintext
matrices $X$ and $W$, respectively, and compute the convolution results as the
inner product between the vectorized versions of the selected image window and
the kernel. Thus, equation (1) can be modified as follows:

$\llbracket\mbox{\boldmath$Y$}\rrbracket_{m,n}=\langle\llbracket\mbox{\boldmath$\widehat{X}$}\rrbracket_{m,n},\llbracket\mbox{\boldmath$\widehat{W}$}\rrbracket\rangle.$
(3)

When there is no padding involved, the above approach requires
$(M-P+1)(N-Q+1)PQ$ ciphertext multiplications (henceforth denoted as CT-CT
mult) and ciphertext additions (henceforth denoted as CT-CT add) and consumes
a multiplicative depth of $1$.

#### Activation and Pooling Operators:

In this paper, we implement a ReLU function using a polynomial approximation.
In (Gottemukkula 2020), a number of polynomial approximations have been
proposed for the ReLU function. Since one of our objectives is to minimize the
computational depth, we choose the following polynomial of degree 2 for our
approximation.

$g(u)=0.47+0.50u+0.09u^{2},u\in[-\sqrt{2},\sqrt{2}].$ (4)

A linear transformation may be required to limit $u$ within the appropriate
range before the above polynomial approximation is applied. The polynomial
approximation of ReLU consumes only a multiplicative depth of 1 and requires 1
CT-CT mult, 2 CT-PT mults (multiplication a ciphertext with a plaintext
operand), 2 CT-CT adds, and 1 CT-PT add (addition of a plaintext operand to a
ciphertext) operations.

In terms of pooling, we consider mean pooling operator in a $2\times 2$
neighborhood, which is a linear operation. While mean pooling does not involve
any ciphertext multiplication (CT-CT mults), it does require $4$ CT-CT adds
and $1$ CT-PT mult. Our activation and pooling operations can be summarized as
follows:

$\llbracket\mbox{\boldmath$Z$}\rrbracket_{i,j}=f\left(\llbracket\mbox{\boldmath$Y$}\rrbracket_{2i-1,2j-1},\llbracket\mbox{\boldmath$Y$}\rrbracket_{2i,2j-1},\llbracket\mbox{\boldmath$Y$}\rrbracket_{2i-1,2j},\llbracket\mbox{\boldmath$Y$}\rrbracket_{2i,2j}\right),$
(5)

where

$f(a,b,c,d)=(g(a)+g(b)+g(c)+g(d))/4.$ (6)

## Experimental Results

### Environment

For FHE, we chose IBM’s open-source library named HElib (Halevi and Shoup
2020) version 1.0.1. Since FHE based experiments require signficant compute
capability and efficiency, we chose the natural API interface of HElib library
through C++. The configuration details of the computing environment are listed
below:

  * •

System used: Virtual Machine ppc64le POWER9,

  * •

CPUs: 112 CPUs with 14 sockets, 1 core/socket, 8 threads per core

  * •

Memory: RAM - 511 GB,

  * •

Operating system: Fedora Core 32

### Impact of FHE parameters

Three key parameters define the operational characteristics of the FHE scheme
and determine the security level ($\lambda$). The plaintext space of CKKS is
the set of polynomials in the cyclotomic ring $\mathbb{Z}[X]/(\Phi_{m}(X))$,
where $\Phi_{m}(X)$ is the $m^{th}$ cyclotomic polynomial with degree given by
Euler’s totient $\phi(m)$. The number of slots ($K$) in the ciphertext is
given by $K=\phi(m)/2$. The second key parameter is $L$, which is the bitsize
of the modulus of a freshly encrypted ciphertext. Since we are using a leveled
HE scheme, $L$ determines the depth of the circuit that can be evaluated
without bootstrapping. Increasing $L$ allows more computations
(multiplications) to be performed before hitting the noise threshold. Finally,
the parameter $r$ determines the computational precision in the encrypted
domain. While larger values of $m$ (and hence $\phi(m)$) increase the security
level $\lambda$, larger values of $L$ and $r$ generally decrease the security.
In our experiments, we choose $m=2^{16}$ (corresponding to $\phi(m)=32768$ and
$K=16384$), $r=35$, and vary $L$ depending on the computational depth
required.

We first evaluate the impact of the depth parameter $L$ on the computational
time, security level, and ciphertext size. For measuring the computational
time, we consider the average time required for the adding/multiplying two
ciphertexts (filled with constant values in all the slots) over 1000 trials.
As shown in Figure 3, increasing $L$ increases the computational time and
ciphertext size in a linear fashion. Moreover, the rate of increase is much
higher for multiplication of two ciphertexts (CT-CT mult) than for the
addition of two ciphertexts (CT-CT add). Our analysis also indicates that for
the first multiplication to be successful, a minimum value of $L=200$ bits is
required and each additional multiplication in the circuit (increase in the
multiplicative depth by 1) requires an increase in the value of $L$ by
approximately 100 bits.

It must be highlighted that the impact of $L$ on the security estimate
($\lambda$) is non-linear. Note that evaluation of deeper (multiplicative
depth) circuits would require large values of $L$. This in turn would require
a large increase in the value of $m$ to maintain the same level of security,
consequently slowing down the computations significantly as well as increasing
the memory requirements (due to larger ciphertext sizes). Therefore, it is
critical to constrain the depth parameter $L$ to a reasonable value and design
the machine learning model to fit the depth constraint.

|  |
---|---|---
(a) | (b) | (c)
Figure 3: Impact of the depth parameter $L$ on (a) computational time, (b)
security level, and (c) ciphertext size. For this experiment, the value of $m$
and $r$ are set to $2^{16}$ and $35$, respectively. The computational time is
the average time required for adding/multiplying two ciphertexts over 1000
trials.

### Dataset and Inference Results

We use the MNIST dataset for our experiments. This dataset consists of $60000$
$28\times 28$ grayscale images of the $10$ digits (0-9), along with a test set
of $10000$ images. We designed a simple CNN model that consists of a single
convolutional layer with $28$ filters (each having a kernel size of $3\times
3$) without any padding, followed by _ActPool1_ layer (polynomial
approximation of ReLU and mean pooling), a flattening layer (whose output
dimension is $13\times 13\times 28=4732$), a single fully connected layer of
size $4732\times 10$, and a final softmax layer. The CNN is trained using the
Keras library for $10$ epochs using a sparse_categorical_crossentropy loss
function and ADAM optimizer with a learning rate of $0.01$. The test accuracy
of this simple CNN model on plaintext images was found to be $97.86$%, which
is only marginally lower than that of a CNN where the approximate ReLU
function is replaced with the standard ReLU function ($98.39$%).

Next, we focus on the problem of encrypted inference where both the inputs to
the CNN and weights of the network are encrypted. To implement inference based
on the above simple CNN model, the depth parameter $L$ was set to $600$ and
the values of $m=2^{16}$ and $r=35$ were retained. The conservative security
estimate based on these parameters is $128$ bits, which is obtained by using
an in-built library call in HElib. Note that the inference process does not
require the computation of the softmax function, which is monotonic in nature.
Thus, the user can decrypt the output of the fully connected layer to infer
the predicted class, which is the class with the maximum output at the FC
layer. Based on the above HE parameters, the number of slots ($K$) available
is more than the number of test images ($10000$). Thus, the complete inference
process could be executed in the encrypted domain in one batch without any
loss of test accuracy ($97.86$%).

The time required for the execution of each layer in the CNN during encrypted
inference is shown in Table 1. It must be emphasized that the reported timings
in Table 1 are for evaluating a single filter on a single thread. The results
for the complete network (all $28$ filters) are shown in the next sub-section
on multi-threading. As one would expect, the convolution layer is primary
bottleneck in terms of number of computations. On the other hand, the fully
connected (FC) layer presents a significant challenge in terms of memory
bottleneck. Even for a single filter, the FC layer requires holding
($(M-P+1)(N-Q+1)(R+1)/4$) ciphertexts in memory, where $R$ is the number of
classes (($R=10$ for MNIST). The number of ciphertexts gets inflated by a
factor of $28$ when all the filters are considered, which makes
parallelization more difficult.

Operation | Execution Time (in seconds) for Single Filter and Single Thread
---|---
Convolution | 487.4
Approximate ReLU | 102.1
Mean pooling | 16.9
Fully Connected | 123.4
Total (including overhead) | 812.6
Table 1: Execution time for different layers of the proposed simple CNN model.

### Parallelization using multi-threading

HElib supports thread level parallelism (inherited from the underlying NTL
library), which can speed up the execution time. Theoretically, the design of
the proposed simple CNN allows us to compute the response of each convolution
filter as well as the subsequent activation, pooling, and FC layers
independently. However, in practice this was not feasible due to memory
constraints. Moreover, our computing environment has 112 CPUs, but the model
has only 28 filters. Therefore, to make the best use of the available
resources we had to implement a nested multi-threading strategy.

In the nested multi-threading strategy, we spawn $F$ (we choose $F=14$ or
$F=28$) threads for the filters, and each filter thread in turn spawns $C$
convolution threads (we choose $C=1,3,5$ or $7$). Each convolutional thread
operates on a horizontal partition of the full image. For example, when $C=7$,
the $28\times 28$ image is partitioned into $7$ horizontal sub-images with
appropriate overlap (to avoid boundary issues). Since the approximate RELU is
a pixel-wise operation, it can be applied to each horizontal partition
independently. But since the mean pool operation could span across partitions,
we have to wait for all the convolution threads to complete before meanpool is
applied. Finally, due to memory constraints, we could not proceed with the FC
layer directly. Instead, we waited for all the filter threads to complete and
then spawn $H$ new threads, one for each of the $10$ classes (columns of the
FC layer matrix). As earlier, these $H$ class threads further spawn $J$
channel threads to parallelize the matrix multiplication task along rows of
the FC layer matrix.

Figure 4: Execution time for various number of filter and convolution threads.
In this experiment, we fix the number of class and channel threads to $10$ and
$1$, respectively.

Figure 4 shows the execution times for the above nested multi-threading
strategy. For this experiment, we fix $H=10$ and $J=1$ and vary $F$ and $C$.
From this figure, we can observe that having more filter threads generally
leads to faster execution. Since the maximum number of available threads is
112, we observe that utilizing only 70-80 threads provides the optimal results
because it achieves the best compromise between the number of threads and
thread efficiency. In particular, when $F=28$ and $C=3$, 84 threads are
utilized in total and it results in the lowest execution time of $864$
seconds. Similarly, when the number of channel threads $J$ is increased, the
best results were obtained when $J=7$ leading to a total execution time of
$561$ seconds.

Thus, the proposed multi-threading strategy leads to approximately $40$ times
improvement over the time required for computing all the $28$ filters on a
single thread ($28\times 812[Table1]=22736$ seconds). This is achieved by
utilizing $70-80$ threads on average. If the memory constraints can be
surmounted, there is scope for further refinement in the proposed multi-
threaded strategy.

It must be emphasized that the values reported in the above paragraphs are the
actual execution times and not the amortized time. Since $K$ images are
processed in parallel exploiting the SIMD mechanism, the amortized time per
image could be orders of magnitude lower. Since $K=16384$ in our experiments,
the amortized execution time for inference on a single image can be estimated
as $34$ milliseconds. Note that most of the reported literature on inference
based on encrypted MNIST data (e.g., (Gilad-Bachrach et al. 2016; Chou et al.
2018)) use the BFV encryption scheme available in the Microsoft SEAL library
implemented on a x86 platform. Hence, it is not possible to directly compare
the execution times reported in this work with those reported in the
literature. However, we do note that the wall-clock run time reported in
(Gilad-Bachrach et al. 2016) is 250 seconds, with only 10 $5\times 5$ filters
and a stride length of 2, which is comparable to our run time of 561 seconds
(with 28 filters and a stride length of 1). Moreover, unlike (Gilad-Bachrach
et al. 2016), the weights of the convolution and fully-connected layers are
encrypted in our proposed implementation.

It must be highlighted that the ciphertext packing strategy used in this work
has high latency for a single inference. This could be addressed using
different packing strategies as discussed in (Brutzkus, Gilad-Bachrach, and
Elisha 2019), which can reduce the latency at the cost of decreasing the
throughput. For example, for single image inference with a different packing
strategy (the entire image is encrypted within a single ciphertext), we were
able to achieve a total inference time of 8.8 seconds with encrypted model
parameters. When the model parameters are not encrypted, the latency can be
further reduced to 2.5 seconds, which is comparable to the latency of 2.2
seconds reported in (Brutzkus, Gilad-Bachrach, and Elisha 2019).

## Conclusions and Future Work

In this paper, we attempted to address the computation complexity in a
convolutional deep neural network for encrypted processing. The approach
presented in the paper describes a secure way to compute convolution, non-
linear activation and pooling layers using power of the CKKS FHE scheme. We
exploited two features to tame the complexity: SIMD ciphertext packing and
thread level parallelism. In addition, we empirically set the depth of the
circuit to further minimize the execution time. Using a recent method to
approximate ReLU, we were able to build a bootstrap-free ReLU function in FHE.
Since the max pooling layer requires bootstrapping, we replaced it with a mean
pooling layer. Overall, we show that it is possible to implement encrypted
inference in reasonable time by investing in the right CNN design and
parameter choices.

In future, we will proceed to integrate other layers with the end-goal of
support full inferencing and training on encrypted data. We also observe that
an overall latency of $561$ seconds could be too high for some applications.
The limitation of the ciphertext packing strategy used in this work is high
latency for a single inference and high memory requirement. This could be
addressed using different packing strategies as discussed in (Brutzkus, Gilad-
Bachrach, and Elisha 2019), which can reduce the latency while also decreasing
the throughput. Since all the ideas contained in (Brutzkus, Gilad-Bachrach,
and Elisha 2019) are practical HE implementation nuances that do not require
any change to the machine learning model, they can be readily applied to the
proposed approach. We leave this analysis for future work.

## Acknowledgement

The work was initiated when authors Nandakumar/Ratha/Pankanti were part of IBM
Research. We thank IBM, USA for providing us the computing facilities. Author
Kumar acknowledges Infosys Foundation for the financial support to IIIT-
Bangalore through the Infosys Foundation Career Development Chair Professor.

## References

  * Acar et al. (2018) Acar, A.; Aksu, H.; Uluagac, A. S.; and Conti, M. 2018.  A Survey on Homomorphic Encryption Schemes: Theory and Implementation.  _ACM Computing Surveys_ 51(4).
  * Al-Rubaie and Chang (2019) Al-Rubaie, M.; and Chang, J. 2019.  Privacy-Preserving Machine Learning: Threats and Solutions.  _IEEE Security & Privacy_ 17: 49–58.  doi:10.1109/MSEC.2018.2888775.
  * Bost et al. (2015) Bost, R.; Popa, R.; Tu, S.; and Goldwasser, S. 2015.  Machine Learning Classification over Encrypted Data.  doi:10.14722/ndss.2015.23241.
  * Brakerski (2012) Brakerski, Z. 2012.  Fully homomorphic encryption without modulus switching from classical GapSVP.  In _Annual Cryptology Conference_ , 868–886. Springer.
  * Brakerski, Gentry, and Vaikuntanathan (2014) Brakerski, Z.; Gentry, C.; and Vaikuntanathan, V. 2014.  (Leveled) fully homomorphic encryption without bootstrapping.  _ACM Transactions on Computation Theory (TOCT)_ 6(3): 1–36.
  * Brutzkus, Gilad-Bachrach, and Elisha (2019) Brutzkus, A.; Gilad-Bachrach, R.; and Elisha, O. 2019.  Low latency privacy preserving inference.  In _International Conference on Machine Learning_ , 812–821.
  * Cheon et al. (2017) Cheon, J. H.; Kim, A.; Kim, M.; and Song, Y. 2017.  Homomorphic encryption for arithmetic of approximate numbers.  In _International Conference on the Theory and Application of Cryptology and Information Security_ , 409–437. Springer.
  * Chillotti et al. (2020) Chillotti, I.; Gama, N.; Georgieva, M.; and Izabachène, M. 2020.  TFHE: fast fully homomorphic encryption over the torus.  _Journal of Cryptology_ 33(1): 34–91.
  * Chou et al. (2018) Chou, E.; Beal, J.; Levy, D.; Yeung, S.; Haque, A.; and Fei-Fei, L. 2018.  Faster CryptoNets: Leveraging Sparsity for Real-World Encrypted Inference.  _CoRR_ abs/1811.09953.  URL http://arxiv.org/abs/1811.09953.
  * Dahl et al. (2018) Dahl, M.; Mancuso, J.; Dupis, Y.; Decoste, B.; Giraud, M.; Livingstone, I.; Patriquin, J.; and Uhma, G. 2018.  _Private Machine Learning in TensorFlow using Secure Computation_.
  * Dathathri et al. (2019) Dathathri, R.; Saarikivi, O.; Chen, H.; Laine, K.; Lauter, K.; Maleki, S.; Musuvathi, M.; and Mytkowicz, T. 2019.  CHET: An Optimizing Compiler for Fully-Homomorphic Neural-Network Inferencing.  In _PLDI 2019_ , 142–156. ACM.  URL https://www.microsoft.com/en-us/research/publication/chet-an-optimizing-compiler-for-fully-homomorphic-neural-network-inferencing/.
  * Fan and Vercauteren (2012) Fan, J.; and Vercauteren, F. 2012.  Somewhat Practical Fully Homomorphic Encryption.  Cryptology ePrint Archive, Report 2012/144.  https://eprint.iacr.org/2012/144.
  * Gentry (2009) Gentry, C. 2009.  Fully homomorphic encryption using ideal lattices.  In _Proceedings of the forty-first annual ACM symposium on Theory of computing_ , 169–178.
  * Gilad-Bachrach et al. (2016) Gilad-Bachrach, R.; Dowlin, N.; Laine, K.; Lauter, K.; Naehrig, M.; and Wernsing, J. 2016.  Cryptonets: Applying neural networks to encrypted data with high throughput and accuracy.  In _International Conference on Machine Learning_ , 201–210.
  * Gottemukkula (2020) Gottemukkula, V. 2020.  Polynomial Activation Functions.  ICLR 2020 Open Review.  https://openreview.net/attachment?id=rkxsgkHKvH&name=original˙pdf.
  * Halevi and Shoup (2020) Halevi, S.; and Shoup, V. 2020.  _HElib - An Implementation of homomorphic encryption_.  https://github.com/shaih/HElib/ (accessed March 30, 2020).
  * Howard et al. (2017) Howard, A. G.; Zhu, M.; Chen, B.; Kalenichenko, D.; Wang, W.; Weyand, T.; Andreetto, M.; and Adam, H. 2017.  MobileNets: Efficient Convolutional Neural Networks for Mobile Vision Applications.  _CoRR_ abs/1704.04861.  URL http://arxiv.org/abs/1704.04861.
  * Juvekar, Vaikuntanathan, and Chandrakasan (2018) Juvekar, C.; Vaikuntanathan, V.; and Chandrakasan, A. 2018.  GAZELLE: A low latency framework for secure neural network inference.  In _27th $USENIX$ Security Symposium ($USENIX$ Security 18)_, 1651–1669.
  * Krizhevsky, Sutskever, and Hinton (2012) Krizhevsky, A.; Sutskever, I.; and Hinton, G. E. 2012.  ImageNet Classification with Deep Convolutional Neural Networks.  In Pereira, F.; Burges, C. J. C.; Bottou, L.; and Weinberger, K. Q., eds., _Advances in Neural Information Processing Systems_ , 1097–1105. Curran Associates, Inc.
  * Lou and Jiang (2019) Lou, Q.; and Jiang, L. 2019.  SHE: A Fast and Accurate Deep Neural Network for Encrypted Data.  In _Advances in Neural Information Processing Systems_ , 10035–10043.
  * Phan et al. (2020) Phan, H.; Huynh, D. T.; He, Y.; Savvides, M.; and Shen, Z. 2020.  MoBiNet: A Mobile Binary Network for Image Classification.  In _The IEEE Winter Conference on Applications of Computer Vision (WACV)_.
  * Poggio et al. (2017) Poggio, T.; Mhaskar, H.; Rosasco, L.; Miranda, B.; and Liao, Q. 2017.  Why and when can deep-but not shallow-networks avoid the curse of dimensionality: a review.  _International Journal of Automation and Computing_ 14(5): 503–519.
  * Simonyan and Zisserman (2014) Simonyan, K.; and Zisserman, A. 2014.  Very Deep Convolutional Networks for Large-Scale Image Recognition.  _arXiv 1409.1556_ .
  * Szegedy et al. (2015) Szegedy, C.; Liu, W.; Jia, Y.; Sermanet, P.; Reed, S.; Anguelov, D.; Erhan, D.; Vanhoucke, V.; and Rabinovich, A. 2015.  Going Deeper with Convolutions.  In _Computer Vision and Pattern Recognition (CVPR)_.  URL http://arxiv.org/abs/1409.4842.
  * Zeiler and Fergus (2013) Zeiler, M.; and Fergus, R. 2013.  Visualizing and Understanding Convolutional Neural Networks.  volume 8689.


# Effects of radial electric field on ion temperature gradient driven mode
stability

Ningfei Chen, Hanyuan Hu, Xiangyu Zhang, Shizhao Wei and Zhiyong Qiu111E-mail:
<EMAIL_ADDRESS>Institute for Fusion Theory and Simulation and Department of
Physics, Zhejiang University, Hangzhou 310027, China

###### Abstract

The local stability of ion-temperature gradient driven mode (ITG) in the
presence of a given radial electric field is investigated using gyrokinetic
theory and ballooning mode formalism with toroidal effect accounted for. It is
found that, zero frequency radial electric field induced poloidal rotation can
significantly stabilize ITG, while the associated density perturbation has
little effect on ITG stability due to the modification of finite-orbit-width
effect. However, the parallel mode structure is slightly affected due to the
evenly symmetric density modulation of ZFZF.

##  I Introduction

Drift waves (DWs) turbulence Horton (1999), driven by free energy associated
with plasma pressure gradients, are considered as candidates for inducing
anomalous plasma transport and degradation of confinement in magnetically-
confined fusion (MCF) devices. Ion-temperature gradient driven mode (ITG) is
one of the most intensively studied DWs due to its potential role in causing
anomalous ion thermal transport, which is much concerned in future fusion
reactors. ITG has two branches, i.e., a slab branch by the coupling of ion
parallel compression and diamagnetic drift, and a toroidal branch by the
coupling of diamagnetic drift with the unfavored curvature in the weak field
side Cheng and Chen (1980); Chen et al. (1991). In-depth understanding of the
mechanisms for ITG linear stability, nonlinear evolution and eventual
saturation, is needed for quantitative understanding of plasma confinement in
future tokamaks. Excitation of zonal flows (ZFs), is considered as an
important route for ITG self-regulation, and the regulation is achieved via
nonlinear excitation of ZFs by ITG via modulation instability as ITG amplitude
exceeds the threshold induced by frequency mismatch, which in turn, scatters
ITG into the linearly stable short radial wavelength regime Chen et al.
(2000); Zonca and Chen (2008).

ZFs are typically meso-scale radial corrugations with toroidally symmetric
($n=0$), and predominantly poloidally symmetric ($m\simeq 0$) scalar potential
fluctuation, and consist of zero-frequency ZF (ZFZF) Diamond et al. (2005) and
its finite frequency counter-part, geodesic acoustic mode (GAM) Winsor et al.
(1968); Zonca and Chen (2008). Here, $m$/$n$ are the poloidal/toroidal mode
numbers of the torus. The nonlinear interaction of ITG with ZFs are observed
in experiments Zhao et al. (2006); Conway et al. (2008); Lan et al. (2008);
Zhong et al. (2015); Melnikov et al. (2017), as well as in large scale
simulations Lin et al. (1998); Hahm et al. (2000); Liu et al. (2010); Liao et
al. (2016) where ITG and the associated transport are suppressed, and it is
also found the threshold on pressure gradient on ITG stability is up-shifted
as nonlinear effects are taken into account Dimits et al. (2000). Furthermore,
the radial electric field $E_{r}$ associated with large-scale mean flow, as
well as its gradient, is also observed to be related to turbulence suppression
and confinement improvement, and possibly the formation of transport barrier,
as well as transition from low- to high-confinement regime.

Several models are proposed to investigate the mechanism of ITG suppression by
ZFs, among which are $E\times B$ shear effect and radial envelope modulation.
In the $E\times B$ shearing model, a two-point nonlinear theory is proposed to
understand flow shear suppression of ITG turbulence in cylindrical and
toroidal geometry, and both high/low-frequency component of the radial
electric field is considered Hahm and Burrell (1995); Hahm et al. (1999). It
is found that a significant reduction of turbulence activity occurs when the
shearing rate $\omega_{S}$, which is proportional to
$\partial(v_{\theta}/r)/\partial_{r}$, with $v_{\theta}=E_{r}/B$ being the
radial electric field induced $E\times B$ drift velocity, exceeds the
decorrelation rate of the ambient turbulence. Additionally, compared to the
low-frequency $E_{r}$, the high-frequency component of $E_{r}$ typically plays
less significant role in turbulence suppression due to its oscillating nature.
On the other hand, the radial envelope modulation model may be referred to as
the scattering process, i.e., the potential well associated with $E_{r}$
scatters turbulence into linearly stable short radial wavelength domain. The
thresholds on the ITG amplitude for the ZFZF and/or GAM generation, however,
are comparable to each other Zonca and Chen (2008). Global descriptions are
needed in both models, which have previously been based on the local
description Romanelli and Zonca (1993). More specifically, the modulational
instability describing nonlinear interaction between ITG and ZFZF requires the
ITG linear dispersion relation with finite-$\theta_{k}$ modification to the
local one. Here, $\theta_{k}\equiv k_{r}/(n\partial_{r}q)$ with $k_{r}$ being
the radial envelope wavenumber, and $q$ being the safety factor. Investigation
of the interaction between turbulence and radial electric field in terms of
the parallel mode structure is not found in the literature to date. Typically,
this issue is treated in the gyrokinetic framework with toroidal effects
neglected Artun and Tang (1992); Gao et al. (2004).

In this work, a local model is proposed to figure out the local properties of
turbulence suppression by given $E_{r}$ and study the “linear” stability of
ITG in the presence of $E_{r}$. Here, “local” means the ITG eigenmode equation
is solved along the magnetic field lines, with toroidal effects and parallel
compression properly accounted for, while physics associated with radial
envelope is neglected systematically. Technically, this is achieved by
deriving an ITG governing eigenmode equation in the existence of the radial
electric field induced density modulation as well as poloidal rotation, which
is then solved in ballooning space for the ITG local dispersion relation and
mode structure Connor et al. (1979); Taylor (1976). For simplicity of
discussion and because of the temporal scale separation, ZFZF-type $E_{r}$ on
ITG local stability is investigated as an example Chen et al. (2000).
Analogous to the ZFZF, finite-frequency GAM effects on ITG stability may be
treated similarly if time scale separation between ITG and GAM is satisfied,
which is typically the case as we discuss in the final section. Finally, our
model may also shed light on turbulence suppression by the mean flow, whose
mechanism is not yet fully understood.

The rest of the paper is organized as follows. In section II, the ITG
eigenmode equation in the presence of a given radial electric field is derived
using gyrokinetic theory and ballooning mode representation. In section III,
the ITG stability is investigated assuming a radial electric field with zero
frequency, i.e., that of ZFZF, in both the short- and long-wavelength limit.
Summary and discussions are given in Section IV.

##  II General formalism

For simplicity of discussion while focusing on the main scope of the present
paper, we consider a tokamak with axisymmetric concentric circular magnetic
surface and straight field line, and a left-handed coordinate
$\left(r,\theta,\phi\right)$ is adopted, with $r$, $\theta$ and $\phi$ being
the minor radius, poloidal and toroidal angles of the torus, respectively. The
equilibrium magnetic field is given as
$\mathbf{B}=B_{0}[(1-\epsilon\cos\theta)\mathbf{e_{\phi}}+\epsilon\mathbf{e_{\theta}}/q]$,
where $\epsilon\equiv r/R$ is the inverse aspect ratio, $R$ is the on-axis
major radius and $q\equiv rB_{\phi}/(RB_{\theta})$. ITG modes generally have
ballooning structure with high mode numbers, and the characteristic scale of
equilibrium profile is generally much larger than the distance between
neighbouring mode rational surfaces. Consequently, the perturbed quantity can
be expressed as

$\displaystyle\delta\phi$ $\displaystyle=e^{in\phi-i\omega
t}\underset{j}{\sum}\hat{\Phi}\left(s-j\right)e^{-i\left(m_{0}+j\right)\theta}.$
(1)

Here, $s\equiv(r-r_{0})/\Delta r=nq-m_{0}$, $r_{0}$ denotes the reference
rational surface with $nq(r_{0})-m_{0}=0$, $\Delta r=1/(n\partial q/\partial
r)$ is the distance between neighboring mode rational surfaces, and
$\left|j\right|\ll m_{0}$ is an integer.

The gyrokinetic equation Taylor and Hastie (1968) is used to investigate the
stability of ITG turbulence in the presence of a given radial electric field.
Following Ref. Chen et al. (1991), we take the flat density gradient limit to
focus on effects of ion temperature gradient, i.e., assuming
$\eta_{i}=L_{ni}/L_{Ti}\rightarrow\infty$, with $L_{ni}=-n_{i}/(\partial
n_{i}/\partial r)$ and $L_{Ti}=-T_{i}/(\partial T_{i}/\partial r)$ being the
characteristic scale length of ion density and temperature nonuniformity,
respectively. The gyrokinetic equation for ion response to ITG can be written
as

$\displaystyle\left(\omega-
k_{\parallel}v_{\parallel}+\omega_{D}+k_{\theta}v_{\theta}\right)\delta
H_{I}=$
$\displaystyle\dfrac{e}{T_{i}}J_{0}\left(\omega+\omega_{*i}^{T}\right)F_{0i}\delta\phi_{I}$
$\displaystyle-i\frac{c}{B_{0}}\mathbf{b}\times\nabla\delta\phi_{I}\cdot\nabla\delta
H_{E}.$ (2)

Here, $k_{\parallel}\equiv\left(nq-m\right)/\left(qR\right)$ is the parallel
wavenumber, $\omega_{D}\equiv
2\omega_{d}C\left(x_{\perp}^{2}/2+x_{\parallel}^{2}\right)$ is the magnetic
drift frequency, with $\omega_{d}=k_{\theta}cT_{i}/(eBR)$,
$x_{\perp}=v_{\perp}/v_{ti}$ and $x_{\parallel}=v_{\parallel}/v_{ti}$ being
the ion perpendicular/parallel velocities normalized by thermal velocity
$v_{ti}=\sqrt{2T_{i}/m_{i}}$, respectively. $C=\cos\theta-
k_{r}\sin\theta/k_{\theta}$ is related to the curvature with $k_{r}$ and
$k_{\theta}=m_{0}/r_{0}$ being the radial/poloidal mode numbers. $\delta
H_{I}$ is the nonadiabatic ion response to ITG,
$J_{0}\left(k_{\perp}\rho_{i}\right)$ is the Bessel function of zero-index
accounting for Finite Larmor radius (FLR) effects,
$\rho_{s}=mv_{\perp,s}c/(eB)$ is the Larmor radius of species $s$, $F_{0i}$ is
the equilibrium ion distribution function, and
$\omega_{*i}^{T}=\omega_{*Ti}\left(x_{\perp}^{2}+x_{\parallel}^{2}-3/2\right)$
is the ion diamagnetic frequency in the flat density limit, with
$\omega_{*Ti}=k_{\theta}cT_{i}/(eBL_{Ti})$. Furthermore, the last term on the
left hand side accounts for the Doppler shift from the radial electric field
induced poloidal rotation, with
$v_{\theta}\equiv-c\mathbf{b}\times\partial_{r}\delta\phi_{E}/B_{0}$, while
the last term on the right hand side represents the perturbed diamagnetic term
associated with the density perturbation induced by the radial electric field
Chen et al. (2014). The last term, in fact, can be combined with the term
proportional to $\omega^{T}_{*i}$ , considering $F_{0}+\delta H_{E}$ to be the
renormalized equilibrium in the existence of the radial electric field. It is
worth noting that the two additional terms, i.e., radial electric field
induced poloidal rotation of ITG $k_{\theta}v_{\theta}\delta H_{I}$, and
variation along the magnetic field line induced by the density perturbation
associated with the radial electric field $\propto\delta\phi_{I}\delta H_{E}$,
can also be obtained from the perpendicular nonlinear term in nonlinear
gyrokinetic equation Frieman and Chen (1982), and thus, will be called
“nonlinear terms” in the following discussion for convenience, though the
radial electric field can also originate from linear effects, such as large
scale mean flow. Here, subscripts “E” and “I” represent quantities associated
with radial electric field and ITG, respectively. The dispersion relation can
be derived from charge quasi-neutrality condition

$\dfrac{eN_{0}\delta\phi}{T_{e}}+\left\langle\delta
H_{e}J_{0}\right\rangle=-\dfrac{eN_{0}\delta\phi}{T_{i}}+\left\langle\delta
H_{i}J_{0}\right\rangle,$ (3)

with $eN_{0}\delta\phi/T_{e}$ and $-eN_{0}\delta\phi/T_{i}$ being adiabatic
responses of electron and ion, respectively, and
$\left\langle\cdots\right\rangle$ representing velocity space integration. The
derivation follows closely the procedure of Ref. Guzdar et al. (1983). For
typical ITG fluctuation with
$k_{\parallel}v_{\parallel,e}\gg\omega\sim\omega_{*i}\gg\omega_{d}$,
$k_{\parallel}v_{\parallel,i}$, electrons responde adiabatically, i.e.,
$\delta H_{I,e}=0$. The nonadiabatic ion response can be derived as

$\displaystyle\delta
H_{I,i}\approx\dfrac{\Lambda}{\omega}\left[-\left(\dfrac{e}{T_{e}}+\dfrac{e}{T_{i}}\right)F_{0}\delta\phi_{E}+\delta
H_{E}\right]\delta\phi_{I}$
$\displaystyle+\dfrac{e}{T_{i}}J_{0}F_{0}\left(1+\dfrac{\omega_{*i}^{T}}{\omega}\right)\left(1+\dfrac{k_{\parallel}v_{\parallel}}{\omega}+\dfrac{k_{\parallel}^{2}v_{\parallel}^{2}}{\omega^{2}}-\dfrac{\omega_{D}}{\omega}\right)\delta\phi_{I}.$
(4)

Here, $\Lambda\equiv ick_{r}k_{\theta}/B_{0}$, and will be used in the rest of
the paper. The two terms in first bracket of equation (4) are the formal
nonlinear terms, and represent the effects associated, respectively, with the
potential and density fluctuation of $E_{r}$. The quasi-neutrality condition
of ITG is applied to simplify the first term. Substituting the ion and
electron response into quasi-neutrality condition (3), one then has the ITG
WKB dispersion relation

$\displaystyle\left\\{\dfrac{1}{\tau\left(1+\omega_{*Ti}/\omega\right)}+b_{\perp}-\dfrac{k_{\parallel}^{2}v_{ti}^{2}}{2\omega^{2}}+\dfrac{2\omega_{d}C}{\omega}+\dfrac{\Lambda}{\omega+\omega_{*Ti}}\right.$
$\displaystyle\left.\times\left[\left(1+\dfrac{1}{\tau}\right)\delta\phi_{E}-\dfrac{T_{i}}{e}\left\langle\dfrac{\delta
H_{E,i}}{N_{0}}\right\rangle\right]\right\\}\delta\phi=0,$ (5)

with $b_{\perp}\equiv k_{\perp}^{2}\rho_{i}^{2}/2$, and $k_{\perp}$ being the
perpendicular wavenumber. The first four terms of equation (5) constitute the
linear ITG dispersion relation, with the first three terms being respectively,
adiabatic electron response, the FLR effect (polarization) and parallel
compressibility, while the forth term related to magnetic drift peculiar in
toroidal configuration, resulting in coupling of neighbouring poloidal
harmonics. The last two terms are nonlinear modifications due to poloidal
rotation and density modulation associated with the radial electric field,
respectively. Noting $k_{\perp}^{2}=k_{\theta}^{2}-\partial^{2}/\partial
r^{2}$, the eigenmode equation in real space for $j$-th poloidal harmonics can
be derived as

$\displaystyle\left(b_{\theta}\hat{s}^{2}\dfrac{d^{2}}{dz^{2}}-\dfrac{1}{\tau\left(1+\omega_{*Ti}/\omega\right)}-b_{\theta}+\dfrac{k_{\parallel}^{2}v_{ti}^{2}}{2\omega^{2}}\right)\hat{\Phi}_{z}$
$\displaystyle=\dfrac{\omega_{d}}{\omega}\left[\hat{\Phi}_{z+1}+\hat{\Phi}_{z-1}+\hat{s}\dfrac{d}{dz}\left(\hat{\Phi}_{z+1}-\hat{\Phi}_{z-1}\right)\right]$
$\displaystyle+\dfrac{\Lambda}{\left(\omega+\omega_{*Ti}\right)}\left[\left(1+\dfrac{1}{\tau}\right)\delta\phi_{E}-\left\langle\dfrac{T_{i}\delta
H_{E,i}}{eN_{0}}\right\rangle\right]\hat{\Phi}_{z}.$ (6)

Here, $\hat{s}\equiv r\left(\partial q/\partial r\right)/q$ is the magnetic
shear, $\tau\equiv T_{e}/T_{i}$, $b_{\theta}\equiv
k_{\theta}^{2}\rho_{ti}^{2}/2$, $z\equiv s-j=nq-m$ is the normalized distance
to the mode rational surface. The first term on the right-hand side of
equation (6) comes from the curvature drift induced coupling between
neighbouring poloidal harmonics. Moreover, the term proportional to
$\left\langle\delta H_{E,i}\right\rangle$ may also have poloidal-dependence,
and causes additional toroidal coupling. For instance, GAM with
$\omega\gg\omega_{tr,i}$, is characterised by up-down anti-symmetric
($\propto\sin\theta$) density fluctuation Zonca and Chen (2008), while ZFZF
with $\omega\ll\omega_{tr,i}$, has $\cos\theta$-type density fluctuation Chen
et al. (2014). Here $\omega_{tr,i}\equiv v_{\parallel,i}/(qR)$ is the
circulating ion transit frequency. Equation (6) can be analyzed using the
ballooning mode formalism framework Connor et al. (1979), which is
accomplished by taking
$\Phi\left(\eta\right)=\int\hat{\Phi}\left(z\right)\exp\left(-i\eta
z\right)dz$, with $\eta$ being the extended poloidal angle along the magnetic
field lines. The ITG eigenmode equation in ballooning space reads

$\displaystyle\dfrac{d^{2}\Phi\left(\eta\right)}{d\eta^{2}}+q^{2}\Omega^{2}b\left(\dfrac{\tau\Omega}{1+\tau\Omega\epsilon_{Ti}^{1/2}}+b\left(1+\hat{s}^{2}\eta^{2}\right)\right.$
$\displaystyle+\dfrac{2}{\Omega}\left(\cos\eta+\hat{s}\eta\sin\eta\right)+\left(1+\tau\right)\Delta_{E}\delta\phi_{E}$
$\displaystyle\left.-\tau\Delta_{E}\left\langle\dfrac{T_{i}\delta
H_{E,i}}{eN_{0}}\right\rangle\right)\Phi\left(\eta\right)=0,$ (7)

where $\Omega\equiv\omega/(\tau\sqrt{\omega_{*Ti}\omega_{d}})$, $b\equiv\tau
b_{\theta}/\sqrt{\epsilon_{Ti}}$, $\epsilon_{Ti}\equiv L_{Ti}/R$ and
$\Delta_{E}\equiv\Lambda/[\omega_{*Ti}\sqrt{\epsilon_{Ti}}\left(1+\tau\Omega\sqrt{\epsilon_{Ti}}\right)]$.
Equation (7) is general and can be applied to study the nonlinear modification
of any given radial electric field to ITG stability, with the nonlinear
modifications by the radial electric field accounted for by the last two
terms. In this work, as a proof of principle demonstration, we will consider
ZFZF-type stationary radial electric field, while the effects of energetic
particle induced GAM (EGAM)/GAM can be investigated straightforwardly
following the same approach if the GAM/EGAM frequency is smaller than ITG
growth rate. It is natural to take the dominant $m=0,1$ components of
nonadiabatic ion response Qiu et al. (2008),

$\displaystyle\delta
H_{E,i}=\dfrac{eF_{0}\overline{\delta\phi_{Z}}}{T_{i}}\left(1+i\lambda\cos\theta\right),$
(8)

and $m=0$ component of $\delta\phi_{Z}$, i.e., $\overline{\delta\phi_{Z}}$,
where
$\left(\overline{\cdots}\right)\equiv\int_{0}^{2\pi}(\cdots)d\theta/2\pi$
represents surface averaged quantity. The higher order $m=1$ density
perturbation of ZFZF is included, to account for its unique role in inducing
periodic modification to the ITG eigenmode potential well along the magnetic
field line, that determines the condition for ITG stability. Here,
$\lambda=\omega_{Dr}/\omega_{tr,i}$ represents the finite drift orbit width
effect, with $\omega_{Dr}=2\omega_{dr}(x_{\perp}^{2}/2+x^{2}_{\parallel})$
being the magnetic drift frequency associated with normal curvature and
$\omega_{dr}\equiv k_{r}cT_{i}/(eBR)$.

##  III Effects of zero-frequency $E_{r}$ on ITG linear stability

With the specified expression of $\delta H_{E,i}$ presented in equation (8),
equation (7) can be written as

$\displaystyle\dfrac{d^{2}\Phi}{d\eta^{2}}+q^{2}\Omega^{2}b\left(\dfrac{\tau\Omega}{1+\tau\Omega\epsilon_{Ti}^{1/2}}+b\left(1+\hat{s}^{2}\eta^{2}\right)+\Delta_{E}\overline{\delta\phi_{Z}}\right.$
$\displaystyle\left.+\dfrac{2}{\Omega}\left(\cos\eta+\hat{s}\eta\sin\eta\right)+\Delta_{E}^{\prime}\overline{\delta\phi_{Z}}\cos\eta\right)\Phi=0,$
(9)

where $\Delta_{E}\overline{\delta\phi_{Z}}$ represents modification due to the
electrostatic potential, and
$\Delta_{E}^{\prime}\overline{\delta\phi_{Z}}\cos\eta$ results from the $m=1$
density perturbation of ZFZF, with $\Delta_{E}^{\prime}\equiv-i\tau
A\left(2\omega_{dr}/\omega_{t,i}\right)\Delta_{E}$ and $A\equiv\left\langle
F_{0}(x_{\perp}^{2}/2+x_{\parallel}^{2})/x_{\parallel}\right\rangle$. Note
that, the velocity space integral $A$ vanishes for typical Maxwellian
distribution, but finite value for non-even symmetric distribution, e.g.,
shifted-Maxwellian distribution with non-zero average parallel velocity due to
auxiliary current drive. Equation (9) will be investigated in both the short-
and long- wavelength limits, corresponding to strong and moderate ballooning
cases, respectively, as investigated Refs. Chen et al. (1991); Guzdar et al.
(1983). The two limiting parameter regimes, can be studied by taking $b\sim
O(1)$ and $b\ll 1$, respectively, as we shown after equation (10) that, the
mode width in ballooning space, is proportional to $1/\sqrt{b}$ (and thus,
$\sim O(\sqrt{b})$ radially).

###  III.1 Short-wavelength limit

In the short-wavelength limit, i.e., $b\sim O\left(1\right)$, the
eigenfunction is strongly localized in ballooning space Chen et al. (1991).
Thus, strong coupling approximation can be adopted by taking $\cos\eta\approx
1-\eta^{2}/2$ and $\sin\eta\approx\eta$ Taylor (1976). Note that, the
assumption underlying the above strong coupling approximation is that the mode
is localized around $\eta=0$, and the introduction of the $\cos\eta$-type
periodic modulation does not affect the validity of the assumption due to the
even-symmetric of it. The eigenmode equation then becomes

$\displaystyle\dfrac{d^{2}\Phi}{d\eta^{2}}+q^{2}\Omega^{2}b\left(\dfrac{\tau\Omega}{1+\tau\Omega\epsilon_{Ti}^{1/2}}+b+\dfrac{2}{\Omega}+\Delta_{E}\overline{\delta\phi_{Z}}+\Delta_{E}^{\prime}\overline{\delta\phi_{Z}}\right)$
$\displaystyle\left.+\left(b\hat{s}^{2}+\dfrac{2\hat{s}-1}{\Omega}-\dfrac{\Delta_{E}^{\prime}\overline{\delta\phi_{Z}}}{2}\right)\eta^{2}\right)\Phi=0,$
(10)

which can be rewritten as a standard Weber equation with the most unstable
ground eigenmode being given by $\delta\phi=\exp(-\sigma\eta^{2})$ with

$\displaystyle\sigma$ $\displaystyle=$
$\displaystyle\dfrac{q^{2}\Omega^{2}b}{2}\left(\dfrac{\tau\Omega}{1+\tau\Omega\epsilon_{Ti}^{1/2}}+\dfrac{2}{\Omega}+b+\Delta_{E}\overline{\delta\phi_{Z}}+\Delta_{E}^{\prime}\overline{\delta\phi_{Z}}\right).$

The half width of the lowest eigenmode in $\eta$ space is proportional to
$1/\sqrt{b}$. The corresponding dispersion relation is

$\displaystyle
q^{2}\Omega^{2}b\left[\dfrac{\tau\Omega}{1+\tau\Omega\epsilon_{Ti}^{1/2}}+\dfrac{2}{\Omega}+b+\Delta_{E}\overline{\delta\phi_{Z}}+\Delta_{E}^{\prime}\overline{\delta\phi_{Z}}\right]^{2}$
$\displaystyle+\left(b\hat{s}^{2}+\dfrac{2\hat{s}-1}{\Omega}-\dfrac{\Delta_{E}^{\prime}\overline{\delta\phi_{Z}}}{2}\right)=0.$
(11)

The dispersion relation is similar to corresponding linear result Chen et al.
(1991), except terms proportional to $\Delta_{E}$ and $\Delta_{E}^{\prime}$
originate from the contribution of radial electric field induced poloidal
rotation and density fluctuation, respectively. The dependence of ITG growth
rate and real frequency on the radial electric field are solved from the
theoretical dispersion relation equation (11), which are then compared with
the numerical solution of equation (9), and good agreement between analytical
and numerical results are obtained, as shown in Fig. 1a and 1b, respectively.
It is found in equation (11) that the $m=1$ density perturbation is of order
$k_{r}\rho_{i}q$ compared to $m=0$ potential fluctuation. Thus, it is obvious
that $m=1$ density fluctuation affects the dispersion relation slightly. The
ITG growth rate decreases significantly with increasing $\delta\phi_{Z}$. We
then analyze the contribution of the radial electric field induced poloidal
rotation and density modulation on ITG stability, by turning off the
corresponding terms in equation (9). It is shown in Fig. 2 that, when the
$E_{r}$ induced poloidal rotation is kept while the density perturbation is
turned off, the ITG growth rate is almost the same as that with both effects
properly accounted for; while as only the $E_{r}$ induced density perturbation
is kept, the ITG growth rate is slightly affected by the scalar potential. We
thus conclude that the reduction of the growth rate is mainly due to the
potential fluctuation (poloidal rotation). Besides, it is found that the ITG
growth rate is of order $C_{s}/L_{T}$, which is much larger than
characteristic GAM/EGAM frequency $\sim C_{s}/R$, hence our analysis can also
be applied to modulation of ITG by EGAM/GAM. The mode structure is also shown
in Fig. 3, and it is clearly seen that the mode structure peaks at $\eta=0$
and the even symmetry is not broken, resulting from the even-symmetric period
modulation introduced by the density fluctuation of ZFZF ($\propto\cos\theta$)
as denoted by $\Delta^{\prime}_{E}$. We note that, equation (9) can be further
simplified, by substituting the quasi-neutrality condition of ZFZF into
equation (7) to replace the last term proportional to $\langle\delta
H_{E,i}\rangle$. This process will introduce $O(k^{2}_{r}\rho^{2}_{i})$
uncertainty since it is $\langle J_{0}\delta H_{E,i}\rangle$ in the quasi-
neutrality condition instead of the $\langle\delta H_{E,i}\rangle$ in equation
(7). In this case, effects induced by the $m=1$ density modulation of ZFZF
cannot be larger than $O(k^{2}_{r}\rho^{2}_{i})$, and is thus weak, as shown
by our numerical results.

(a)

(b)

Figure 1: The dependence of normalized growth rate (a) and real frequency (b)
of ITG, which are normalized by $C_{s}/L_{Ti}$, on the normalized ZFZF
intensity $e\delta\phi_{Z}/T_{i}$. The squares represent the theoretical
result given by Eq. (11) while circles are numerical result of Eq. (9). Here,
$C_{s}^{2}\equiv 2T_{e}/m_{i}$ is the sound speed, $\epsilon_{Ti}=0.2$, $b=1$
and $A=1$. Figure 2: Numerical result of integrated and separated effects of
potential and density fluctuation of ZFZF. Figure 3: The mode structure of the
lowest eigenmode when $e\delta\phi_{Z}/T_{i}=0.1$. Here, $\epsilon_{Ti}=0.2$,
$b=1$ and $A=1$.

###  III.2 Long-wavelength limit

For typical tokamak plasmas, strong coupling approximation is usually a crude
constraint. In more general cases, $b\ll 1$ (long-wavelength limit) is
satisfied, and strong coupling approximation no longer holds. In the long-
wavelength limit, there are two branches, i.e., toroidal branch and slab
branch. We are more concerned about the toroidal branch Chen et al. (1991),
which is characterized by fast variation over connection length scale
($\eta\sim O\left(1\right)$) and a superimposed slowly varying envelope over
secular scale. The self-consistent ordering is given by balancing parallel
compressibility and adiabatic electron response, which results in
$\Omega=O\left(b^{-1/3}\right)$. Taking
$\Phi\left(\eta\right)=C_{0}\left(\eta_{1}\right)\cos\eta/2+S_{0}\left(\eta_{1}\right)\sin\eta/2$
with $\eta_{1}\equiv\hat{\epsilon}\eta$, and $\hat{\epsilon}=b^{1/3}$ denoting
slow variation in $\eta$, the eigenmode equations can be derived from
vanishing coefficients of $\sin\eta/2$ and $\cos\eta/2$:

$\displaystyle\dfrac{dS_{0}}{d\eta_{1}}+\left[\dfrac{b^{2/3}q^{2}\Omega^{3}\tau}{1+\tau\Omega\epsilon_{Ti}^{1/2}}-\dfrac{1}{4b^{1/3}}+q^{2}\Omega^{2}b^{2/3}\Delta_{E}\overline{\delta\phi_{Z}}\right.$
$\displaystyle\left.+\dfrac{1}{2}q^{2}\Omega^{2}b^{2/3}\Delta_{E}^{\prime}\overline{\delta\phi_{Z}}\right]C_{0}+q^{2}\Omega
b^{1/3}\hat{s}\eta_{1}S_{0}=0,$ (12)
$\displaystyle\dfrac{dC_{0}}{d\eta_{1}}-\left[\dfrac{b^{2/3}q^{2}\Omega^{3}\tau}{1+\tau\Omega\epsilon_{Ti}^{1/2}}-\dfrac{1}{4b^{1/3}}+q^{2}\Omega^{2}b^{2/3}\Delta_{E}\overline{\delta\phi_{Z}}\right.$
$\displaystyle\left.+\dfrac{1}{2}q^{2}\Omega^{2}b^{2/3}\Delta_{E}^{\prime}\overline{\delta\phi_{Z}}\right]S_{0}-q^{2}\Omega
b^{1/3}\hat{s}\eta_{1}C_{0}=0.$ (13)

Equations (12) and (13) can be cast into a Weber equation for $C_{0}$ and
$S_{0}$. The dispersion relation for the most unstable ground eigenstate is

$\displaystyle\Omega^{3}+\dfrac{\Lambda\overline{\delta\phi_{Z}}}{\tau\omega_{*T}\epsilon_{Ti}^{1/2}}\left(1-\dfrac{i\tau
Ak_{r}\rho_{i}q}{2}\right)\Omega^{2}-\dfrac{\epsilon_{Ti}^{1/2}}{4bq^{2}}\Omega$
$\displaystyle=$ $\displaystyle\dfrac{1}{4bq^{2}\tau}.$

Here, terms in the bracket of Eq. (LABEL:eq:14) comes from the $m=0$ component
of radial electric field induced poloidal rotation and $m=1$ density
modulation, while other terms originate from linear dispersion relation Chen
et al. (1991). It is noteworthy that the $m=1$ component of density
perturbation enters and affects the ITG dispersion relation, while the
$\sin\theta$-type density modulation of EGAM/GAM has no influence on the
dispersion relation, possibly due to the odd symmetry of the density
modulation of GAM/EGAM with $\omega_{G}\gtrsim\omega_{tr,i}$. The dependence
of ITG growth rate and real frequency on scalar potential of the radial
electric field are solved from the theoretical dispersion relation, which are
then compared with the numerical solution of Eq. (9), and good agreement are
obtained, as shown in Fig. 4. An artificially small $b=0.01$ is adopted to
separate different scales, although this is not the most relevant parameter
regime for ITG stability. As shown in Fig. 5, the $E_{r}$ induced poloidal
rotation is the main reason for the reduction of the ITG growth rate, as
clarified by our theoretical analysis; while its density perturbation has weak
effect on ITG stability. The mode structure is shown in Fig. 6, which is wider
than that in the short-wavelength limit, as demonstrated by our analysis.

(a)  (b)

Figure 4: The growth rate (a) and real frequency (b) of ITG versus the
normalized ZFZF intensity $e\delta\phi_{Z}/T_{i}$. The squares represent the
analytical result given by Eq. (LABEL:eq:14) while circles represent numerical
result of Eq. (9). Here, $\epsilon_{Ti}=0.2$, $b=0.01$ and $A=1$. Figure 5:
Numerical result of integrated and separated effects of ZFZF potential and
density fluctuation. Figure 6: The mode structure of the most unstable mode
when $e\delta\phi_{Z}/T_{i}=0.1$. Here, $\epsilon_{Ti}=0.2$, $b=0.01$ and
$A=1$.

##  IV Conclusion and Discussion

In this paper, a governing equation is formulated to investigate the ITG
“linear” stability in the presence of a given radial electric field, using
gyrokinetic equation and ballooning mode representation. The effects of the
radial electric field on ITG linear stability consist of $E_{r}$-induced
poloidal rotation and density fluctuation, and their integrated and separated
contribution to ITG stability are studied both theoretically and numerically.

Here, ZFZF is presented as an example for clarity of discussion. For the
adopted ZFZF-like radial electric field with $\omega\ll\omega_{tr,i}$, we
found that the poloidal rotation is the main reason for the significant
reduction of the ITG growth rate both in the short- and long-wavelength limit.
In contrast, the up-down symmetric density perturbation, which peaks at the
un-favourable curvature region, have weak suppression effect on ITG turbulence
in both short- and long-wavelength limit. The extension of our ITG stability
analysis to include EGAM/GAM with frequency being much higher than ion transit
frequency, and thus up-down anti-symmetric density perturbation
($\propto\sin\theta$), is straightforward, since in that case the GAM/EGAM
frequency is smaller than ITG growth rate.

The present work, motivated to understand local properties and mechanism of
ITG stability in the existence of a given $E_{r}$ in ballooning space, found
that, the radial electric field always plays a stabilizing role on ITG. This
model can also be applied to study the possible cross-scale interaction
between ITG and AEs, with the former being microscopic and the latter being
macro- or meso- scale, mediated by ZFZF. In such a two-pray one-predator
system, effects of ZFZF generated by one turbulence can be considered as
“passive” or “equilibrium” for another turbulence. Additionally, it is of
great interest and importance to investigate the nonlinear modulation of ITG
by EGAM/GAM, which can be excited by AEs (internally) and neutral beam
injection (externally), and can act as an active control of ITG turbulence.
But, in contrast to ZFZF, the time scale separation of EGAM/GAM with ITG is
not necessarily always satisfied, and depends on specific experimental
conditions. Hence, fully nonlinear process, i.e., DW-GAM nonlinear evolution,
should be considered to account for comparable frequency between ITG and GAM,
and will be investigated in a future publication.

## Acknowledgements

This work is supported by the National Key R&D Program of China under Grant
No. 2017YFE0301900, and the National Science Foundation of China under grant
No. 11875233. The original idea of investigating ITG stability in the
existence of EGAM induced periodic modification to the potential well was
provided by Prof. Liu Chen (Zhejiang University and University of California,
Irvine).

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Horton (1999) W. Horton, Rev. Mod. Phys. 71, 735 (1999).
  * Cheng and Chen (1980) C. Z. Cheng and L. Chen, The Physics of Fluids 23, 1770 (1980).
  * Chen et al. (1991) L. Chen, S. Briguglio, and F. Romanelli, Physics of Fluids B: Plasma Physics 3, 611 (1991).
  * Chen et al. (2000) L. Chen, Z. Lin, and R. White, Physics of Plasmas 7, 3129 (2000).
  * Zonca and Chen (2008) F. Zonca and L. Chen, EPL (Europhysics Letters) 83, 35001 (2008).
  * Diamond et al. (2005) P. H. Diamond, S.-I. Itoh, K. Itoh, and T. S. Hahm, Plasma Physics and Controlled Fusion 47, R35 (2005).
  * Winsor et al. (1968) N. Winsor, J. L. Johnson, and J. M. Dawson, The Physics of Fluids 11, 2448 (1968).
  * Zhao et al. (2006) K. J. Zhao, T. Lan, J. Q. Dong, L. W. Yan, W. Y. Hong, C. X. Yu, A. D. Liu, J. Qian, J. Cheng, D. L. Yu, et al., Phys. Rev. Lett. 96, 255004 (2006).
  * Conway et al. (2008) G. D. Conway, C. Troster, B. Scott, K. Hallatschek, and the ASDEX Upgrade Team, Plasma Physics and Controlled Fusion 50, 055009 (2008).
  * Lan et al. (2008) T. Lan, A. D. Liu, C. X. Yu, L. W. Yan, W. Y. Hong, K. J. Zhao, J. Q. Dong, J. Qian, J. Cheng, D. L. Yu, et al., Physics of Plasmas 15, 056105 (2008).
  * Zhong et al. (2015) W. Zhong, Z. Shi, Y. Xu, X. Zou, X. Duan, W. Chen, M. Jiang, Z. Yang, B. Zhang, P. Shi, et al., Nuclear Fusion 55, 113005 (2015).
  * Melnikov et al. (2017) A. Melnikov, L. Eliseev, S. Lysenko, M. Ufimtsev, and V. Zenin, Nuclear Fusion 57, 115001 (2017).
  * Lin et al. (1998) Z. Lin, T. S. Hahm, W. W. Lee, W. M. Tang, and R. B. White, Science 281, 1835 (1998).
  * Hahm et al. (2000) T. S. Hahm, K. H. Burrell, Z. Lin, R. Nazikian, and E. J. Synakowski, Plasma Physics and Controlled Fusion 42, A205 (2000).
  * Liu et al. (2010) F. Liu, Z. Lin, J. Q. Dong, and K. J. Zhao, Physics of Plasmas 17, 112318 (2010).
  * Liao et al. (2016) X. Liao, Z. Lin, I. Holod, B. Li, and G. Y. Sun, Physics of Plasmas 23, 122305 (2016).
  * Dimits et al. (2000) A. M. Dimits, G. Bateman, M. A. Beer, B. I. Cohen, W. Dorland, G. W. Hammett, C. Kim, J. E. Kinsey, M. Kotschenreuther, A. H. Kritz, et al., Physics of Plasmas 7, 969 (2000).
  * Hahm and Burrell (1995) T. S. Hahm and K. H. Burrell, Physics of Plasmas 2, 1648 (1995).
  * Hahm et al. (1999) T. S. Hahm, M. A. Beer, Z. Lin, G. W. Hammett, W. W. Lee, and W. M. Tang, Physics of Plasmas 6, 922 (1999).
  * Romanelli and Zonca (1993) F. Romanelli and F. Zonca, Physics of Fluids B: Plasma Physics 5, 4081 (1993).
  * Artun and Tang (1992) M. Artun and W. M. Tang, Physics of Fluids B: Plasma Physics 4, 1102 (1992).
  * Gao et al. (2004) Z. Gao, J. Q. Dong, and H. Sanuki, Physics of Plasmas 11, 3053 (2004).
  * Connor et al. (1979) J. Connor, R. Hastie, and J. Taylor, Proceedings of the Royal Society of London series a-mathematical physical and engineering sciences 365, 1 (1979).
  * Taylor (1976) J. Taylor, Plasma Physics and Controlled Nuclear Fusion Research 2, 323 (1976).
  * Taylor and Hastie (1968) J. B. Taylor and R. J. Hastie, Plasma Physics 10, 479 (1968).
  * Chen et al. (2014) L. Chen, Z. Qiu, and F. Zonca, EPL (Europhysics Letters) 107, 15003 (2014).
  * Frieman and Chen (1982) E. A. Frieman and L. Chen, The Physics of Fluids 25, 502 (1982).
  * Guzdar et al. (1983) P. N. Guzdar, L. Chen, W. M. Tang, and P. H. Rutherford, The Physics of Fluids 26, 673 (1983).
  * Qiu et al. (2008) Z. Qiu, L. Chen, and F. Zonca, Plasma Physics and Controlled Fusion 51, 012001 (2008).