# Enhancing Psychotherapy Counseling: A Data Augmentation Pipeline Leveraging
Large Language Models for Counseling Conversations

Jun-Woo Kim♣ Ji-Eun Han♣ Jun-Seok Koh  Hyeon-Tae Seo  Du-Seong Chang∗
KT
{jun_woo.kim, ji-eun.han, js.koh, ht.seo<EMAIL_ADDRESS>
♣ Equal contributions
* Corresponding author

###### Abstract

We introduce a pipeline that leverages Large Language Models (LLMs) to
transform single-turn psychotherapy counseling sessions into multi-turn
interactions. While AI-supported online counseling services for individuals
with mental disorders exist, they are often constrained by the limited
availability of multi-turn training datasets and frequently fail to fully
utilize therapists’ expertise. Our proposed pipeline effectively addresses
these limitations. The pipeline comprises two main steps: 1) Information
Extraction and 2) Multi-turn Counseling Generation. Each step is meticulously
designed to extract and generate comprehensive multi-turn counseling
conversations from the available datasets. Experimental results from both
zero-shot and few-shot generation scenarios demonstrate that our approach
significantly enhances the ability of LLMs to produce higher quality multi-
turn dialogues in the context of mental health counseling. Our pipeline and
dataset are publicly available here.

##  1 Introduction

In contemporary society, the prevalence of mental illness is rising, requiring
expert counseling. As large language models (LLMs) like OpenAI’s ChatGPT
demonstrate ease of access and potential for counseling, increasing numbers of
people are interested in using these tools as private counselors Lin et al.
(2023); De Choudhury et al. (2023).

However, existing AI-assisted chatbot services, typically designed for
everyday conversation, often fall short due to inadequate training and fail to
match the quality of responses provided by human experts. Another challenge is
the limited availability of multi-turn counseling datasets. Bots trained on
single-turn interactions are often incapable of offering practical solutions.
While single-turn datasets are more readily available, there is a scarcity of
high-quality multi-turn counseling data. Existing multi-turn psychotherapy
counseling datasets, such as HOPE Malhotra et al. (2022) and MEMO Srivastava
et al. (2022), are valuable resources. However, these datasets have not fully
utilized the unique counseling styles that every expert possesses. The
effectiveness of psychological counseling is enhanced when both the counselor
and the client have a better understanding of each other’s information. This
can be understood as the creation of high-quality counseling data when such
information is provided.

To address aforementioned challenges, we propose a pipeline that generates
multi-turn counseling conversations based on the client’s and
psychotherapist’s information. Recent research has shown that LLMs can be
applied to data augmentation tasks Zheng et al. (2022); Dai et al. (2023).
These models exhibit excellent capabilities in generating synthetic text data
that reflect various generation conditions. Therefore, we utilize this ability
to augment the source dataset. Our objective is to integrate experts’ distinct
counseling styles into the construction of a synthetic multi-turn counseling
dataset, resulting in a resource that is both more realistic and practically
applicable.

Our contributions are threefold:

  * •

We present a novel pipeline for augmenting psychotherapy multi-turn counseling
data augmentation using LLMs. This pipeline leverages the characteristics of
both clients and psychotherapists to generate practical counseling data.

  * •

We release an augmented dataset of multi-turn counseling chat that
incorporates details of the client’s mental illness and the psychotherapist’s
counseling characteristics.

  * •

We demonstrate the effectiveness of our data augmentation pipeline in
enhancing the performance of LLMs.

##  2 Related Work

The application of LLMs in psychological counseling has garnered significant
interest in recent years. Several studies are being conducted to enable LLMs
to replace the role of counselors in psychological counseling Chung et al.
(2023); Liu et al. (2023); Fu et al. (2023); Lai et al. (2023).

Building on the foundational advancements in LLMs, recent research has
increasingly focused on the development of specialized datasets to enhance the
performance of LLMs in psychological counseling. The effectiveness of these
models in counseling applications is significantly influenced by the quality
and relevance of the training data. Therefore, creating specialized datasets
that encapsulate various aspects of psychological counseling is crucial for
improving the models’ ability to generate contextually appropriate and
empathetic responses Inaba et al. (2024); Qiu et al. (2023); Chen et al.
(2023); Li et al. (2024); Bertagnolli (2020). High-quality datasets in
psychotherapy counseling can be likened to transcripts of effective
psychotherapy conversations. Research suggests that effective psychotherapy
stems from a good relationship between the psychotherapist and the client
Herman (1998); Tschuschke et al. (2022).

Our research extends these efforts by proposing a pipeline for transforming
single-turn psychological counseling dataBertagnolli (2020) into high-quality
multi-turn datasets and validating the effectiveness of the constructed data.
We adopt a multi-faceted approach that includes collecting information from
both counselors and clients, as well as gathering data on various mental
health disorders.

##  3 Preliminary

###  3.1 Task Definition

The augmentation is based on the source dataset from
CounselChat111https://towardsdatascience.com/counsel-chat-bootstrapping-high-
quality-therapy-data-971b419f33da. We selected this dataset as our source data
because it includes information such as responses from multiple therapists to
the same counseling session and the preference voting results for those
responses. Given the source dataset $D_{i}={(x_{i},y_{i},m_{i})}$, where
$x_{i}$ represents client’s utterance, $y_{i}$ represents therapist’s response
to the client’s utterance and $m_{i}$ represents client’s mental disorder. Our
aim is to augment the source data $D_{i}$ into a multi-turn dataset
$D_{i}^{\prime}={(x^{\prime}_{i},y^{\prime}_{i},m_{i},c_{i},t_{i})}$. Here,
$x^{\prime}_{i}=(x_{i}^{1},x_{i}^{2},\dots,x_{i}^{k})$ represents the
augmented client’s utterances derived from the base $x_{i}$ into $k$ multi-
turns. $y^{\prime}_{i}=(y_{i}^{1},y_{i}^{2},\dots,y_{i}^{k})$ represents the
augmented therapist’s responses derived from the base $y_{i}$ into $k$ multi-
turns. Therefore, the augmented $x_{i}^{1}$ and $y_{i}^{1}$ are $x_{i}$ and
$y_{i}$ from the source data, respectively. Additionally, $c_{i}$ contains the
client’s information, and $t_{i}$ is the therapist’s counseling
characteristics both of which are extracted from the source dataset.

###  3.2 Source Dataset Pre-processing

To maximize the utilization of therapist information, we preprocess the source
dataset. According to the latest South Korean government
report222https://www.ncmh.go.kr/mentalhealth/main.do, which identifies severe
stress, continuous depression, anxiety, and sleep disorders, as the most
prevalent mental illnesses in 2022, we have selected these conditions as our
focus for an augmented dataset.

We hypothesize that therapists with more recommendations have provided
appropriate and beneficial responses to clients. Therefore, we select the
expert who has received the most recommendations among the responses of
various experts to the client utterance previously identified for the mental
disorder. To verify the simplicity and feasibility of our approach, we decided
to focus on a single expert per client utterance. Since the information is
extracted from therapists who have provided good responses to specific mental
disorder-related questions, this extracted therapist information can be
valuable for data augmentation. The characteristics of the selected
therapists’ counseling sessions are utilized to augment the source dataset
into a high-quality multi-turn dataset.

Figure 1: Overview of the proposed data augmentation pipeline.

##  4 Method

We propose an augmentation pipeline for extending single-turn psychotherapy
counseling data into multi-turn counseling dialogue. The pipeline is shown in
Figure 1. Our pipeline consists of two steps: 1) Information extraction and 2)
Multi-turn dialogue generation. In the second phase, all four sub-prompts are
used to generate multi-turn counseling dialogue. Among various publicly
available LLMs, we utilize the instruction-tuned version of Llama3-70B AI@Meta
(2024) to generate multi-turn dialogue. Each step will be explained in detail
in the following sections.

###  4.1 Information Extraction

For the first step, we focus on extracting key information from the source
data. Given the client’s utterance $x_{i}$, the therapist’s response $y_{i}$,
and the client’s mental disorder $m_{i}$, we extract the inherent information
from both the client and therapist in the single-turn dialogue. Additionally,
we extract a description of the client’s mental disorder. This information is
crucial for deriving multi-turn dialogues from the source data, enabling for
the construction of more realistic dialogues.

###  4.2 Multi-turn Counseling Generation

In this step, we construct a prompt to extend single-turn counseling into
multi-turn counseling dialogue. It is composed of four sub prompts: 1)
Description prompt, 2) Condition prompt, 3) Information Prompt, and 4) Answer
Prompt.

####  4.2.1 Description Prompt

The description prompt outlines an overview of the content.

The following is a transcript of a chat between a psychotherapist and a client
about {client’s mental disorder}.

In the {client’s mental disorder} field, we place $m_{i}$ from the source
dataset.

####  4.2.2 Condition Prompt

The condition prompt sets the guidelines for the conditions that the model
should follow.

The client starts the conversation as [client] and the psychotherapist starts
the conversation as [psychotherapist]. Please use the dialog and speakers info
as a guide to continue your consultation like #format#. Never create anything
other than the #format# and just complete the “utterance” part.

This clarifies who initiates the conversation, what information the model
should utilize, and the generation format the model should follow.

####  4.2.3 Information Prompt

Using the information extracted from Section 4.1, we first provide the
client’s mental disorder information. This allows the model to generate a
counseling conversation that takes into account the client’s symptoms. Then,
we include a single-turn, $x_{i}$ and $y_{i}$, from the source data to direct
the model on what it should generate.

####  4.2.4 Answer Prompt

In the answer prompt, we instruct the model on the proper format for
generating the conversation.

#format#
[client]:“utterance”
[psychotherapist]: “utterance”

The model should generate appropriate utterances in the “utterance” field.
This structured format reduces confusion when the model generates its answer
and aids in the post-processing of generated text.

##  5 Augmented Dataset

The augmented dataset’s distribution of mental disorder categories, generated
by our pipeline, is presented in Table 1, which lists the number of cases for
each disorder.

Table 1: Mental disorder categories of augmented data Mental Disorder of Client | Number of Cases
---|---
Depression | 69
Anxiety | 45
Anger Management | 16
Trauma | 13

Depression is the most prevalent, with 69 cases, indicating it is the most
common issue addressed. Anxiety follows with 45 cases, showing a significant
presence but less than Depression. Anger Management and Trauma are less
common, with 16 and 13 cases, respectively. This distribution reflects the
original category distribution of the source dataset.

##  6 Experiment

We conduct zero-shot and few-shot experiments to demonstrate the effectiveness
of our pipeline for generating synthetic multi-turn counseling conversations.
For the zero-shot experiment, we did not incorporate any specific
psychological counseling data into the model and only configured the
psychotherapy role prompts. In contrast, for the few-shot experiment, we
utilized the dataset we constructed as input. The specific examples used in
the actual experiments are detailed in section 6.1.3 .

###  6.1 Experiment Details

####  6.1.1 Test Dataset

We select 70 dialogues from augmented data randomly. Table 2 presents the
mental disorder categories of clients in the test data.

Table 2: Mental disorder categories of test data Mental Disorder of Client | Number of Cases
---|---
Depression | 34
Anxiety | 22
Anger Management | 8
Trauma | 6

####  6.1.2 Baseline models

We use Llama2-7B, Llama3-70B as baseline models for generating multi-turn
counseling dialogues. We chose these models because they are open-sourced and
easily reproducible. We also compare these results with the Llama2-7B model
that has been fine-tuned on the Counsel Chat dataset, which is available on
HuggingFace 333https://huggingface.co/NadunAnjanaka/Llama-2-7b-chat-
Counsellor.

####  6.1.3 Experimental setting

We evaluate two experimental settings: 1) zero-shot multi-turn dialogue
generation (see Figure 2 for example), 2) few-shot multi-turn dialogue
generation (see Figure 3 for example). Results are compared within the same
mental disorder category. We consider these settings because well-constructed
synthetic multi-turn counseling data should improve the model’s few-shot
dialogue generation capabilities.

The following is a transcript of a chat between a psychotherapist and a client
about depression. The client starts the conversation as [client] and the
psychotherapist starts the conversation as [psychotherapist]. Please complete
new transcript about [Question].

[Question]
[client] I’m almost never happy. Half of the time, I don’t feel anything. I
find it easy to make myself feel nothing. I know I push people away because
it’s easier. I just want answers. I’m sick of feeling this way. It’s ruining
my relationships with people.
[psychotherapist] Figure 2: Example of zero-shot prompt The following is a
transcript of a chat between a psychotherapist and a client about depression.
The client starts the conversation as [client] and the psychotherapist starts
the conversation as [psychotherapist]. Please use the following [Example] as a
guide complete new transcript about [Question].

[Example]
[client] They don’t go away, and I feel like I’m going crazy. Does that ever
stop? Can it be a symptom of medication?
[psychotherapist] Since you realize that hearing voices in your head is not
usual for you, then definitely there is a problematic situation happening
within your awareness of who you are.if you recently started taking a new drug
or increased dosage of one you already were taking, and the voices started
shortly after, then yes, it is possible medication created your problem.Start
by telling whoever gave you the presecription, about the problem youŕe
having.”Crazy” has some flexibility as to whether someone is this way or
not.Certainly a very positive sign that youŕe not crazy, is that youŕe self-
aware of a problem within yourself. And, youŕe responsible toward yourself and
making effort to address this problem.Crazy people usually dont́ do
responsible behaviors.
[client] I’ve been taking the same medication for a while now, but the dosage
was increased a few weeks ago. Could that be the cause of the voices?
[psychotherapist] That’s a good point. The dosage increase could definitely be
a contributing factor. It’s possible that your body is reacting to the higher
dosage in a way that’s causing these symptoms. I would still recommend
reporting this to your prescribing doctor, as they can help you determine the
best course of action.

[Question]
[client] I have been dealing with depression and anxiety for a number of
years. I have been on medication, but lately my depression has felt worse. Can
counseling help?
[psychotherapist] Figure 3: Example of few-shot prompt

####  6.1.4 Automatic Evaluation

We conduct an automatic evaluation using GPT-4o, the latest model introduced
by OpenAI. The evaluation prompt is described in Figure 4. We follow the LLM
evaluation suggested by Chiang and Lee (2023) for automatic evaluation based
on prompting. The prompt involves rating two answers from different models on
a scale of 1 to 5, with higher values indicating better answers.

You’re an assistant who evaluates answers strictly from the psychotherapist’s
perspective about {mental disorder category}. Please rate [Answer 1] and
[Answer 2] for the consultation [Question], respectively. Rate the two answers
on a scale of 1-5, with higher values indicating better answers. Figure 4:
Example of evaluation prompt

###  6.2 Result

####  6.2.1 Evaluation of Zero-shot and Few-shot Performance

Table 3 presents the average scores evaluated by GPT-4o for the two models,
Llama2-7B-Chat and Llama3-70B-Instruct, in zero-shot and few-shot settings.
The average score in the zero-shot category is 3.814 for Llama2-7B-Chat,
slightly lower than the 4.042 for Llama3-70B-Instruct, indicating marginally
better performance by Llama3-70B-Instruct when no prior examples are provided.
In the few-shot setting, Llama2-7B-Chat has an average score of 4.557, while
Llama3-70B-Instruct achieves 4.785, demonstrating that both models perform
better when examples generated by our pipeline are provided.

Table 3: Average evaluation score evaluated by GPT-4o Model | Zeroshot Avg. | Fewshot Avg.
---|---|---
Llama2-7B-Chat | 3.814 | 4.557
Llama2-70B-Instruct | 4.042 | 4.785
Figure 5: Comparison of zero-shot and few-shot multi-turn counseling dialogue
generation performance for Llama2-7B-chat and Llama3-70B-Instruct. In the few-
shot setting, examples generated by our pipeline are used.

Figure 5 compares the win rates of two models, Llama2-7B-chat and
Llama3-70B-Instruct, in terms of their few-shot multi-turn generation
performance. Zero-shot wins are shown in green, few-shot wins in peach, and
ties in light grey. Each count reflects the win rate comparison between the
two models’ zero-shot and few-shot capabilities. Llama2-7B-chat recorded 14
zero-shot wins, 55 few-shot wins, and 1 Tie. Llama3-70B-Instruct, on the other
hand, had 11 zero-shot wins, 56 few-shot wins, and 3 ties. The results suggest
that Llama3-70B-Instruct performs slightly better in few-shot multi-turn
generation. In both Llama2 and Llama3 model, using few-shot examples
constructed by our pipeline contributes to generating better quality of multi-
turn counseling conversation. Using few-shot examples constructed by our
pipeline enhances the models’ ability to generate high-quality multi-turn
counseling conversations in both Llama2 and Llama3.

####  6.2.2 Evaluation of Zero-shot and Few-shot Performance based on Mental
Disorder Categories

Figure 6: Comparison of zero-shot and few-shot multi-turn counseling dialogue
generation performance across mental disorder categories for Llama2-7B-Chat
and Llama3-70B-Instruct.

Figure 6 shows a comparative analysis of win counts by category for
Llama2-7B-Chat (upper) and Llama3-70B-Instruct (lower). Win counts are
categorized into zero-shot (green), few-shot (peach), and tie (grey) for
mental disorders including Depression, Anxiety, Anger-management, and Trauma.
Llama2-7B-Chat demonstrates the highest win counts in the few-shot category,
with 28 wins in Depression and 15 in Anxiety. Similarly, Llama3-70B-Instruct
shows robust few-shot performance, recording 27 wins in Depression and 18 in
Anxiety. This comparison shows that both models are particularly effective in
few-shot scenarios across all disorders. The overall trend highlights the
superior performance of few-shot dialogue generation for both models in
addressing various psychological disorders.

####  6.2.3 Evaluation of Performance Based on Data Used

Table 4: Win rate and average evaluation score evaluated by GPT-4o Model | Win rate | Avg.
---|---|---
Llama2-7B-Chat-CounsellorBertagnolli (2020) | 0.014 | 1.828
Llama2-7B-Chat Few-shot | 0.985 | 4.528

The experimental setup involved comparing Llama2-7B-Chat-Counsellor and
Llama2-7B-Chat Few-shot, with the results illustrated in Table 4.

Table 4 demonstrates a clear comparison between the two models’ performance.
Llama2-7B-Chat-Counsellor achieved an average score of 1.828, while
Llama2-7B-Chat Few-shot significantly outperformed Llama2-7B-Chat-Counsellor
with an average score of 4.528. These results indicate that Llama2-7B-Chat
Few-shot is markedly more effective in generating high-quality, contextually
appropriate, and supportive responses in extended multi-turn interactions
compared to Llama2-7B-Chat-Counsellor.

This substantial difference in average scores underscores the effectiveness of
the pipeline proposed in our research for transforming single-turn data into
high-quality multi-turn datasets. The enhanced performance of Llama2-7B-Chat
Few-shot suggests that the comprehensive approach we employed, which includes
collecting detailed information from both counselors and clients and
incorporating data on various mental health disorders, is crucial for
improving the capabilities of LLMs in psychological counseling scenarios.

##  7 Conclusion

In this paper, we experimentally assess the utility of the pipeline we propose
by comparing multi-turn psychotherapy counseling dialogues generated using our
pipeline with those generated by models trained solely on original data. We
demonstrate that extracting implicit information from the original data to use
as input for multi-turn generation aids in producing high-quality multi-turn
dialogues. This is analogous to the way effective counseling, which alleviates
clients’ mental disorders, occurs when therapists and clients have a mutual
understanding of each other’s situations. This experiment provides practical
implications not only for researchers and practitioners in the psychotherapy
domain but also for those exploring domains with limited multi-turn dialogue
data, offering a baseline pipeline for data augmentation.

## References

  * AI@Meta [2024] AI@Meta.  Llama 3 model card.  2024\.
  * Bertagnolli [2020] Nicolas Bertagnolli.  Counsel chat: Bootstrapping high-quality therapy data.  by towardsdatascience. com. URL: https://towardsdatascience. com/c ounsel-chat-bootstrapping-high-quality-therapy-data-971b 419f33da (cit. on pp. 72, 96), 2020.
  * Chen et al. [2023] Yirong Chen, Xiaofen Xing, Jingkai Lin, Huimin Zheng, Zhenyu Wang, Qi Liu, and Xiangmin Xu.  Soulchat: Improving llms’ empathy, listening, and comfort abilities through fine-tuning with multi-turn empathy conversations.  In Findings of the Association for Computational Linguistics: EMNLP 2023, pages 1170–1183, 2023.
  * Chiang and Lee [2023] Cheng-Han Chiang and Hung-yi Lee.  Can large language models be an alternative to human evaluations?  In Anna Rogers, Jordan Boyd-Graber, and Naoaki Okazaki, editors, Proceedings of the 61st Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), pages 15607–15631, Toronto, Canada, July 2023. Association for Computational Linguistics.
  * Chung et al. [2023] Neo Christopher Chung, George Dyer, and Lennart Brocki.  Challenges of large language models for mental health counseling.  arXiv preprint arXiv:2311.13857, 2023.
  * Dai et al. [2023] Haixing Dai, Zhengliang Liu, Wenxiong Liao, Xiaoke Huang, Zihao Wu, Lin Zhao, Wei Liu, Ninghao Liu, Sheng Li, Dajiang Zhu, et al.  Auggpt: Leveraging chatgpt for text data augmentation.  arXiv preprint arXiv:2302.13007, 2023.
  * De Choudhury et al. [2023] Munmun De Choudhury, Sachin R Pendse, and Neha Kumar.  Benefits and harms of large language models in digital mental health.  arXiv preprint arXiv:2311.14693, 2023.
  * Fu et al. [2023] Guanghui Fu, Qing Zhao, Jianqiang Li, Dan Luo, Changwei Song, Wei Zhai, Shuo Liu, Fan Wang, Yan Wang, Lijuan Cheng, et al.  Enhancing psychological counseling with large language model: A multifaceted decision-support system for non-professionals.  arXiv preprint arXiv:2308.15192, 2023.
  * Herman [1998] Steven M Herman.  The relationship between therapist–client modality similarity and psychotherapy outcome.  The journal of psychotherapy practice and research, 7(1):56, 1998.
  * Inaba et al. [2024] Michimasa Inaba, Mariko Ukiyo, and Keiko Takamizo.  Can large language models be used to provide psychological counselling? an analysis of gpt-4-generated responses using role-play dialogues.  arXiv preprint arXiv:2402.12738, 2024.
  * Lai et al. [2023] Tin Lai, Yukun Shi, Zicong Du, Jiajie Wu, Ken Fu, Yichao Dou, and Ziqi Wang.  Psy-llm: Scaling up global mental health psychological services with ai-based large language models.  arXiv preprint arXiv:2307.11991, 2023.
  * Li et al. [2024] Anqi Li, Yu Lu, Nirui Song, Shuai Zhang, Lizhi Ma, and Zhenzhong Lan.  Automatic evaluation for mental health counseling using llms.  arXiv preprint arXiv:2402.11958, 2024.
  * Lin et al. [2023] Baihan Lin, Djallel Bouneffouf, Guillermo Cecchi, and Kush R Varshney.  Towards healthy ai: large language models need therapists too.  arXiv preprint arXiv:2304.00416, 2023.
  * Liu et al. [2023] June M Liu, Donghao Li, He Cao, Tianhe Ren, Zeyi Liao, and Jiamin Wu.  Chatcounselor: A large language models for mental health support.  arXiv preprint arXiv:2309.15461, 2023.
  * Malhotra et al. [2022] Ganeshan Malhotra, Abdul Waheed, Aseem Srivastava, Md Shad Akhtar, and Tanmoy Chakraborty.  Speaker and time-aware joint contextual learning for dialogue-act classification in counselling conversations.  In Proceedings of the fifteenth ACM international conference on web search and data mining, pages 735–745, 2022.
  * Qiu et al. [2023] Huachuan Qiu, Anqi Li, Lizhi Ma, and Zhenzhong Lan.  Psychat: A client-centric dialogue system for mental health support.  arXiv preprint arXiv:2312.04262, 2023.
  * Srivastava et al. [2022] Aseem Srivastava, Tharun Suresh, Sarah P Lord, Md Shad Akhtar, and Tanmoy Chakraborty.  Counseling summarization using mental health knowledge guided utterance filtering.  In Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pages 3920–3930, 2022.
  * Tschuschke et al. [2022] Volker Tschuschke, Margit Koemeda-Lutz, Agnes von Wyl, Aureliano Crameri, and Peter Schulthess.  The impact of clients’ and therapists’ characteristics on therapeutic alliance and outcome.  Journal of Contemporary Psychotherapy, 52(2):145–154, 2022.
  * Zheng et al. [2022] Chujie Zheng, Sahand Sabour, Jiaxin Wen, Zheng Zhang, and Minlie Huang.  Augesc: Dialogue augmentation with large language models for emotional support conversation.  arXiv preprint arXiv:2202.13047, 2022.


††thanks: These authors contribute equally to this work.††thanks: These
authors contribute equally to this work.††thanks: Corresponding author:
<EMAIL_ADDRESS>
# Manipulating non-Hermitian skin effect via electric fields

Yi Peng  Shenzhen Institute for Quantum Science and Engineering, Southern
University of Science and Technology, Shenzhen 518055, China  International
Quantum Academy, Shenzhen 518048, China  Guangdong Provincial Key Laboratory
of Quantum Science and Engineering, Southern University of Science and
Technology, Shenzhen 518055, China  Jianwen Jie  Shenzhen Institute for
Quantum Science and Engineering, Southern University of Science and
Technology, Shenzhen 518055, China  International Quantum Academy, Shenzhen
518048, China  Guangdong Provincial Key Laboratory of Quantum Science and
Engineering, Southern University of Science and Technology, Shenzhen 518055,
China  Dapeng Yu  Shenzhen Institute for Quantum Science and Engineering,
Southern University of Science and Technology, Shenzhen 518055, China
International Quantum Academy, Shenzhen 518048, China  Guangdong Provincial
Key Laboratory of Quantum Science and Engineering, Southern University of
Science and Technology, Shenzhen 518055, China  Yucheng Wang  Shenzhen
Institute for Quantum Science and Engineering, Southern University of Science
and Technology, Shenzhen 518055, China  International Quantum Academy,
Shenzhen 518048, China  Guangdong Provincial Key Laboratory of Quantum Science
and Engineering, Southern University of Science and Technology, Shenzhen
518055, China

###### Abstract

In non-Hermitian systems, the phenomenon that the bulk-band eigenstates are
accumulated at the boundaries of the systems under open boundary conditions is
called non-Hermitian skin effect (NHSE), which is one of the most iconic and
important features of a non-Hermitian system. In this work, we investigate the
fate of NHSE in the presence of electric fields by analytically calculating
the dynamical evolution of an initial bulk state and numerically computing the
spectral winding number, the distributions of eigenstates, as well as the
dynamical evolutions. We show abundant manipulation effects of dc and ac
fields on the NHSE, and that the physical mechanism behind these effects is
the interplay between the Stark localization, dynamic localization and the
NHSE. In addition, the finite size analysis of the non-Hermitian system with a
pure dc field shows the phenomenon of size-dependent NHSE. We further propose
a scheme to realize the discussed model based on an electronic circuit. The
results will help to deepen the understanding of NHSE and its manipulation.

Introduction.— Hermiticity of Hamiltonian has been regarded as a fundamental
requirement in standard quantum mechanics, and it ensures the conservation of
probability and limits energy-values to be real in isolated systems. However,
many systems, such as the nonequilibrium and open systems with gain and loss,
can be effectively described by non-Hermitian Hamiltonians. Especially in
recent years, non-Hermitian physics has attracted widespread attention in both
theory Ashida2020 ; Bender ; Hatano ; Lee2016 ; Shen2018 ; Gong2018 ; Yao2018
; Kunst2018 ; Yang2020 ; Lee2019 ; Kawabata2019 ; Longhi2019 ; Jiang2019 ;
Liu2020 ; Longwen2021 ; Bergholtz2021  and experiment Schomerus ; Luo2019 ;
Xiao2020 ; Weidemann2020 ; Helbig2020 ; XXZhang ; Zhang2021 ; Yi2021 . Various
unique features of non-Hermitian systems without any Hermitian counterparts
have been revealed, such as exceptional points and rings Debowski2001 ;
Wiersig2014 ; Hu2017 ; Chen2017 ; Hodaei2017 ; Zhang2019 ; Miri2019 ; Xu2017 ;
Cerjan2019 ; Gong2019 , enriched topological classifications Ueda2019 ; Magnea
; Wojcik2020 ; Li2021 ; Hu2021 ; Liu2019 , and non-Hermitian skin effect
(NHSE) Yao2018 ; Kunst2018 ; Lee2019 ; Helbig2020 ; Alvarez2018 ; Slager2020 ;
Zensen ; Yifei2020 ; Okuma2020 ; LinhuLi . NHSE, namely that a majority of
eigenstates are localized near the boundary under open boundary conditions
(OBC), is one of the most iconic properties of non-Hermitian systems. It
drastically reshapes the bulk-boundary correspondence principle and motivates
the establishment of generalized Brillouin zone Yao2018 ; Kunst2018 ; Yang2020
. The interplay between the NHSE and other fundamental phenomena (e.g.,
localization induced by external magnetic fields, defects, disorder and
quasiperiodic potentials Longhi2019 ; Jiang2019 ; Liu2020 ; Longwen2021 ;
LinhuLi ; Lu2021 ; Shao2021 ; Hughes2021 ; MBLUeda ) has also attracted
widespread attentions recently. On the other hand, electric field can induce
the Stark localization or dynamical localization, and is also a frequently
used fundamental method to manipulate other physical effects, since it is
easily realized and controlled. However, the effect of electric fields on NHSE
was not considered before.

Now we focus on how to manipulate NHSE by using electric fields. If the NHSE
can be fully suppressed in the modulation process, the non-Hermitian effect
may be eliminated, and the system may have the conservation of probability and
all the eigenvalues may become real, even though the non-Hermitian term
remains. Thus, manipulating NHSE is helpful for deepening our understandings
of non-Hermitian quantum mechanics and the differences between Hermitian and
non-Hermitian physics. Moreover, mastering how to control NHSE, we will be
able to obtain or remove it on demand. Therefore, manipulating NHSE also has
practical significance. In this work, we want to address whether electric
fields can manipulate NHSE, and furthermore, if they can, whether richer and
more interesting physics and applications will emerge in light of this.

Model and results.— We consider a one dimensional non-Hermitian system with
non-reciprocal hopping under the influence of electric fields, and the
Hamiltonian is written as

$\hat{H}=\sum_{n}(J_{L}|n\rangle\langle n+1|+J_{R}|n+1\rangle\langle
n|)+E(t)a\sum_{n}n|n\rangle\langle n|,$ (1)

where $|n\rangle$ is the Wannier state localized on the lattice site $n$,
$J_{L}(J_{R})$ represents the leftward (rightward) hopping amplitude, $a$ is
the lattice constant, being set as $1$ throughout this work, and
$E(t)=e\xi(t)$, with $e$ and $\xi(t)$ being the particle’s charge and external
electric field, respectively.

We can analytically confirm the existence or disappearance of NHSE by
investigating the motion of a particle, this is because the particle initially
localized in the bulk should move toward the boundary if the NHSE exists. We
firstly substitute an arbitrary time-dependent quantum state
$|\psi(t)\rangle=\sum_{m}C_{m}(t)|m\rangle$ into the Schrödinger equation
$i\partial_{t}|\psi(t)\rangle=\hat{H}(t)|\psi(t)\rangle$ and obtain the
equation of motion for the time-dependent amplitudes $C_{m}(t)$,

$\displaystyle
i\partial_{t}C_{m}(t)=J_{L}C_{m+1}(t)+J_{R}C_{m-1}(t)+mE(t)C_{m}(t).$ (2)

Here we set $\hbar=1$. By solving the Eq. (2), for arbitrary $E(t)$, we can
obtain the exact solutions,

$\displaystyle C_{m}(t)$ $\displaystyle=$
$\displaystyle\sum_{n}(-1)^{m-n}C_{n}(0)e^{-i\eta(t)n}\mathcal{J}_{m-n}\left(2\sqrt{J_{L}J_{R}\left[\mathcal{U}^{2}(t)+\mathcal{V}^{2}(t)\right]}\right)\left[\frac{J_{R}}{J_{L}}\frac{i\mathcal{V}(t)-\mathcal{U}(t)}{i\mathcal{V}(t)+\mathcal{U}(t)}\right]^{\frac{m-n}{2}}.$
(3)

Here $\mathcal{J}_{m-n}(x)$ is the $(m-n)$th order Bessel function of the
first kind,
$\mathcal{U}(t)=\int_{0}^{t}\cos\left[\eta(t)-\eta(t^{\prime})\right]dt^{\prime}$
and
$\mathcal{V}(t)=\int_{0}^{t}\sin\left[\eta(t)-\eta(t^{\prime})\right]dt^{\prime}$
with

$\eta(t)=\int_{0}^{t}E(t^{\prime})dt^{\prime}.$ (4)

This solution is valid for arbitrary initial bulk state $|\psi(0)\rangle$ and
arbitrary $E(t)$, and the calculation details are in the Supplementary
Materials SM . To simplify the expression, without loss of generality, we
consider a specific initial state that only occupy a single Wannier lattice
site $n_{0}$, namely $|\psi(0)\rangle=|n_{0}\rangle$, and then the probability
at any site $m$ after evolution time $t$, $\rho_{m}(t)=|C_{m}(t)|^{2}$, takes
SM ,

$\displaystyle\rho_{m}(t)=\mathcal{J}_{m-n_{0}}^{2}\left(2\sqrt{J_{L}J_{R}\left[u^{2}(t)+v^{2}(t)\right]}\right)\left(\frac{J_{R}}{J_{L}}\right)^{m-n_{0}},$
(5)

with

$\displaystyle u(t)=\int_{0}^{t}dt^{\prime}\cos\eta(t^{\prime}),\quad
v(t)=\int_{0}^{t}dt^{\prime}\sin\eta(t^{\prime}).$ (6)

In the absence of electric fields, i.e., $E(t)=0$, from Eq. (4) and Eq. (6),
we have $u(t)=t$ and $v(t)=0$, yielding
$\rho_{m}(t)=\mathcal{J}_{m-n_{0}}^{2}\left(2t\sqrt{J_{L}J_{R}}\right)(J_{R}/J_{L})^{m-n_{0}}$,
where $2t\sqrt{J_{L}J_{R}}$ linearly increases to infinity and thus
$\mathcal{J}_{m-n_{0}}(2t\sqrt{J_{L}J_{R}})$ tends to $0$ Bsel  [see Fig.
1(a)], such that $\rho_{m}/\rho_{m-1}\approx J_{R}/J_{L}$ when
$t\rightarrow\infty$ SM . Therefore, the system has right (left) boundary skin
mode when $J_{R}/J_{L}>1$ ($J_{R}/J_{L}<1$), which is consistent with previous
results Yao2018 ; Kunst2018 ; Lee2019 . The following sections will discuss
three cases: i) the pure dc field case, ii) the pure ac field case, and iii)
the dc-ac mixed field case.

Figure 1:  (a) Distributions of the zeroth, first and second order Bessel
function. We can see two characteristics used in the text: the amplitudes of
oscillation decrease with increasing $x$ and when $x\rightarrow\infty$,
$\mathcal{J}_{m}(x)\rightarrow 0$ for any $m$; $\mathcal{J}_{0}(0)=1$ and
$\mathcal{J}_{m}(0)=0$ when $m\neq 0$, so we have
$\mathcal{J}_{m}(0)=\delta_{m,0}$. (b) Localization-delocalization transition
characterized by winding number for finite lattices. The inset shows
transition dc field strength $E_{0}^{c}$ versus $1/L$, given $\gamma=0.77$
marked as dashed line in its parent figure. The distributions of eigenstates
with (c) $L=10$ and (d) $L=20$, and other parameters are $\gamma=0.185$ and
$E_{0}=0.3$ as marked by red star in (b). Dynamical evolution of a electron
started from the lattice center under (e) the weak dc field with $E_{0}=0.005$
and (f) the strong dc field with $E_{0}=0.5$, and other parameters are
$L=160$, $\gamma=0.769$.

The pure dc electric field case.— We firstly discuss the fate of NHSE in the
presence of a pure dc electric field, i.e., $E(t)=E_{0}$. From Eq. (4) and Eq.
(6), we have $u(t)=\sin(E_{0}t)/E_{0}$ and $v(t)=(1-\cos E_{0}t)/E_{0}$, and
then, Eq. (5) gives the probability:

$\displaystyle\rho_{m}(t)=\mathcal{J}_{m-n_{0}}^{2}\left(\frac{4\sqrt{J_{L}J_{R}}}{E_{0}}\sin\frac{E_{0}t}{2}\right)\left(\frac{J_{R}}{J_{L}}\right)^{m-n_{0}}.$
(7)

Note that $\sin(E_{0}t^{*}/2)=0$ at the time points $t^{*}=2\pi N/E_{0}$ with
$N=0,1,2,\cdots$. By using the properties of the Bessel function Bsel
$\mathcal{J}_{m-n_{0}}(0)=\delta_{m-n_{0},0}$ [see Fig. 1(a)] and
$(J_{R}/J_{L})^{m-n_{0}}=1$ when $m=n_{0}$, we have that
$\rho_{m=n_{0}}(t^{*})$ will oscillate back to $1$ whatever the initial site
$n_{0}$ is, and this phenomenon is called Stark localization Wannier1962 ,
which induces that the particle initially localized at the bulk does not move
toward the boundary. Therefore, the effect of the interplay between the Stark
localization and NHSE is that even a small dc field is sufficient to suppress
the NHSE.

The analytical results can be further confirmed by numerically calculating the
winding number (WN). To define the WN, we need to introduce the twist boundary
condition here, i.e., $\hat{H}(\Phi)=\hat{H}+J_{L}e^{i\Phi}|L\rangle\langle
1|+J_{R}e^{-i\Phi}|1\rangle\langle L|$, where $L$ is the system size and
$\Phi$ is the introduced phase factor, and then the WN reads

$w=\frac{1}{2\pi
i}\int_{0}^{2\pi}\partial_{\Phi}\ln\det\left[\hat{H}(\Phi)-\mathcal{E}_{c}\right]\mathrm{d}\Phi.$
(8)

$w=1$ ($w=0$) corresponds to the existence (non-existence) of NHSE under OBC
with eigenvalue around $\mathcal{E}_{c}$ Gong2018 ; Slager2020 ; Zensen ;
Okuma2020 , which is set to the algebra average of the spectrum in the
following calculation. For convenience, we set $J_{L}=J-\gamma/2$ and
$J_{R}=J+\gamma/2$ with $\gamma>0$ and $J=1$ as the unit of energy. Fig. 1 (b)
shows the transition of the existence-nonexistence of NHSE, obtained by
calculating the WN, which changes from $1$ to $0$ when the strength of dc
field $E_{0}$ increases cross the transition line with fixed $L$ from left to
right. It can be seen that the transition lines tend to $E_{0}=0$ with
increasing $L$, which is consistent with the analytical result.

The analytical expression of $\rho_{m}(t)$ can also tell us the finite size
effect, where interesting physics will emerge. Here, we consider
$J_{R}>J_{L}$, which makes the oscillation of the particle favor the right-
hand side of the initial position $n_{0}$, and the oscillation range is
approximate to $4J_{R}/|E_{0}|x_{\star}$ SM , where $x_{\star}$ only depends
on $J_{L}$ and $J_{R}$. If the distance between the right-side boundary and
$n_{0}$ is larger than $4J_{R}/|E_{0}|x_{\star}$, the particle will return
back to $n_{0}$ after a period of time, but if the distance is less than
$4J_{R}/|E_{0}|x_{\star}$, the particle will arrive at the boundary and then
stay there ever since. Thus, for fixed size $L$, there exist a critical
electric field strength $E^{c}_{0}$ that describes the transition of the
existence-nonexistence of NHSE and satisfies $4J_{R}/|E^{c}_{0}|x_{\star}=L$,
giving $E_{0}^{c}\propto 1/L$, as shown in the inset of Fig. 1 (b). When the
size exceed a critical value, the number of skin modes is about
$4J_{R}/|E_{0}|x_{\star}$, being independent of the size. It can be clearly
seen by comparing Fig. 1 (c) and (d), which show the distributions of
eigenstates with different sizes and same other parameters, and have the same
number of skin modes. Thus, with fixing $J_{L}$, $J_{R}$ and $E_{0}$, the
number of skin modes can be of the same order of magnitude as the total
eigenstate number for the system with small size and the NHSE exists. When the
size is large enough, the ratio of the number of skin modes to the total
eigenstate number will be insignificant and the NHSE will disappear,
suggesting that the NHSE is size-dependent, which is different from the
general NHSE. Moreover, we can increase the number of skin modes by decreasing
the electric field strength, and thus, we can control the appearance or
disappearance of NHSE for a finite size system, as shown in Fig. 1 (e) and
(f), where the system shows NHSE when $E_{0}=0.005$, but when $E_{0}=0.5$, the
NHSE disappear. The phenomenon of the size-dependent NHSE should widely exist
in the non-Hermitian systems with defective, disordered, quasiperiodic or
Stark potentials, or two coupled chains with dissimilar non-reciprocal
hoppings LinhuLi .

Figure 2:  (a) Red round dots are zeros of zeroth order of the Bessel function
$\mathcal{J}_{0}(E_{1}/\omega)$, which correspond to the condition of the
emergence of the dynamic localization under ac field driving. Dynamical
evolution of a particle initially localized at the lattice center with (b)
$E_{1}/\omega=2.405$, which is the first zero point of
$\mathcal{J}_{0}(E_{1}/\omega)$ and (c) $E_{1}/\omega=6.1$. Here we set $J=1$,
$\gamma=0.73$, and $\omega=0.46$.

The pure ac electric field case.— We then consider the monochromatic cosine
shape ac electric field $E(t)=E_{1}\cos(\omega t)$. It is a typical Floquet
driving system, with period $T=2\pi/\omega$. From Eq. (4) and Eq. (6), we have
$u(t)=\int_{0}^{t}dt^{\prime}\cos\left(\frac{E_{1}}{\omega}\sin\omega
t^{\prime}\right)$ and
$v(t)=\int_{0}^{t}dt^{\prime}\sin\left(\frac{E_{1}}{\omega}\sin\omega
t^{\prime}\right)$, and then, Eq. (5) gives the probability in the limit $t\gg
T$ SM ,

$\displaystyle\rho_{m}(t\gg
T)\approx\mathcal{J}_{m-n_{0}}^{2}\left(2t\sqrt{J_{L}J_{R}}\mathcal{J}_{0}\left(\frac{E_{1}}{\omega}\right)\right)\left(\frac{J_{R}}{J_{L}}\right)^{m-n_{0}}.$
(9)

When $\mathcal{J}_{0}({E_{1}}/{\omega})\neq 0$,
$2t\sqrt{J_{L}J_{R}}\mathcal{J}_{0}({E_{1}}/{\omega})$ increases linearly to
infinity and thus
$\mathcal{J}_{m-n_{0}}(2t\sqrt{J_{L}J_{R}}\mathcal{J}_{0}({E_{1}}/{\omega}))$
tends to zero, which is completely similar to the case without electric field,
suggesting that the NHSE is not affected by the ac field. For the special ac
field strength $E_{1}$ and frequency $\omega$ that satisfy
$\mathcal{J}_{0}({E_{1}}/{\omega})=0$, corresponding to the red round dots in
Fig. 2(a), due to $\mathcal{J}_{m-n_{0}}(0)=\delta_{m-n_{0},0}$, the particle
initially localized in the bulk will move around the initial position, and
thus, the NHSE will be suppressed, as shown in Fig. 2(b). This localization
phenomenon is called dynamic localization Dunlap1986 ; ZhaoXG ; Longhi2006 ;
Eckardt2009 , which is distinct from Anderson localization induced by random
disorder potential. For most cases, $\mathcal{J}_{0}({E_{1}}/{\omega})\neq 0$,
the particle will hop to the boundary eventually as demonstrated in Fig. 2(c).
To sum up, when only applying the ac electric field to the system, the NHSE is
not affected except for these special parameters of $E_{1}$ and $\omega$ that
causes the dynamic localization, suppressing the NHSE.

Figure 3: (a) The schematic phase diagram about the existence or non-existence
of NHSE under the mixed electric field driving. Red lines represent the
integer $E_{0}/\omega$. Dark green dots on the red lines correspond to
$E_{1}/\omega$ being the zeroes of $\mathcal{J}_{E_{0}/\omega}$, where dynamic
localization (DL) induced by the ac field occurs. The light green region
correspond to the Stark localization (SL). The dynamical evolutions of the
electron initialized on the lattice center under simultaneous driving of dc
and ac fields with (b) $E_{0}/\omega=0.5$ and $E_{1}/\omega=1.3$, (c)
$E_{0}/\omega=1$ and $E_{1}/\omega=3.832$ corresponding to the first zero of
$\mathcal{J}_{1}$, and (d) $E_{0}/\omega=1$ and $E_{1}/\omega=5.7$. Here we
fix $J=1$, $\gamma=0.73$ and $\omega=0.46$.

The dc+ac electric field case.— We finally consider the effect of dc-ac mixed
fields, i.e., $E(t)=E_{0}+E_{1}\cos(\omega t)$, in manipulating NHSE. By using
Eq. (4) and Eq. (6), we can calculate functions $u(t)$ and $v(t)$, whose
expressions look fair complicated SM , and find that when $E_{0}/\omega$ is
not an integer, all the terms in the expressions of $u(t)$ and $v(t)$ are
bounded oscillatory functions of time, meaning that the particle will
oscillate around the initial position. Therefore, the electric fields break
NHSE for non-integer $E_{0}/\omega$. When $E_{0}/\omega$ is an integer, the
probability in the limit of long evolution time SM  can be simplified to

$\rho_{m}(t\gg
T)\approx\mathcal{J}_{m-n_{0}}^{2}\left(2t\sqrt{J_{L}J_{R}}\mathcal{J}_{\frac{E_{0}}{\omega}}\left(\frac{E_{1}}{\omega}\right)\right)\left(\frac{J_{R}}{J_{L}}\right)^{m-n_{0}}.$
(10)

Similar to the discussions for Eq. (S31), the disappearance or existence of
NHSE depends on whether $E_{1}/\omega$ is one of the zeroes of the Bessel
function $\mathcal{J}_{E_{0}/\omega}$. The pure dc field can cause the Stark
localization, which suppresses NHSE. Then adding the ac field
$E_{1}\cos(\omega t)$ with $E_{0}/\omega$ being integers, the particle can
break though the localization barrier and move through the chain accompanied
by the photon absorption or emission. Thus, NHSE can exist only for integer
$E_{0}/\omega$ but except the situations that $E_{1}/\omega$ are the zeroes of
$\mathcal{J}_{E_{0}/\omega}$, which will induce the dynamic localization, as
discussed above. The effects of the dc-ac mixture fields on NHSE are
summarized in the Fig. 3(a). The light green region between the red lines and
the dark green dots on the red lines respectively correspond to the non-
integer $E_{0}/\omega$ and the zeroes of $\mathcal{J}_{E_{0}/\omega}$ with
$E_{0}/\omega$ being integers, which will induce the Stark localization and
dynamic localization and lead to the disappearance of NHSE, as shown in Fig.
3(b) and (c). The red lines correspond to the integer $E_{0}/\omega$, which
can break the bulk localization by photon assisted hopping, and the particle
initially localized in the bulk will move toward the boundary, as shown in
Fig. 3(d).

Electronic circuit’s realization.—The non-Hermitian model (1) can be simulated
by a classical electric circuit as depicted in Fig. 4, which consists of $L$
LC circuit units. Based on the Kirchhoff’s current law, we have

$I_{n}^{R}=I_{n}^{L}+I_{n}^{B},$ (11)

where $I_{n}^{B}$ is the current flows through the $n$th unit and $I_{n}^{L}$
($I_{n}^{R}$) is the current from the $(n-1)$th ($n$th) unit to the $n$th
($(n+1)$th) unit, and they satisfy

$\displaystyle
L_{n}\frac{\mathrm{d}{}I_{n}^{L}}{\mathrm{d}{}t}=V_{n}-V_{n-1},\quad
L_{n+1}\frac{\mathrm{d}{}I_{n}^{R}}{\mathrm{d}{}t}=V_{n+1}-V_{n},$ (12)
$\displaystyle
l_{n}\frac{\mathrm{d}}{\mathrm{d}{}t}\left[I_{n}^{B}-C_{n}\frac{\mathrm{d}{}V_{n}}{\mathrm{d}{}t}\right]=V_{n},\qquad\qquad$
(13)

where $V_{n}$ is the voltage on the node $n$. From Eqs. (11-13), one can
obtain,

$\frac{\mathrm{d}^{2}V_{n}}{\mathrm{d}{}t^{2}}=\frac{V_{n+1}-V_{n}}{C_{n}L_{n+1}}+\frac{V_{n-1}-V_{n}}{C_{n}L_{n}}-\frac{V_{n}}{C_{n}l_{n}}$
(14)

When choosing inductors with inductances $L_{n}=L_{0}g^{-n}$,
$l_{n}=L_{n}/(\Delta-anE_{0})$, and capacitors with capacitance
$C_{n}=C_{0}g^{n}$, Eq. (14) becomes

$\left(1+g+\Delta+\frac{1}{\omega_{0}^{2}}\frac{\mathrm{d}^{2}}{\mathrm{d}{}t^{2}}\right)V_{n}=V_{n-1}+gV_{n+1}+naE_{0}V_{n},$
(15)

where $\omega_{0}=1/\sqrt{C_{0}L_{0}}$. We make a transformation:
$V_{n}\rightarrow V_{n}e^{\pm{}i\omega{}t}$ with $n=1,2,\dots,L$ and
$\omega=\omega_{0}\sqrt{(1+g+\Delta)-E_{R}}$, then Eq. (15) becomes
$E_{R}V_{n}=V_{n-1}+gV_{n+1}+naE_{0}V_{n}$, which describes the dc case of the
model (1) with $J_{L}=1$ and $J_{R}=g$. By directly detecting the eigenvalues
and eigenstates through an elementary voltage measurement Helbig2020 , one can
detect the manipulation effect of dc field.

Figure 4: Schematic of the LC electronic circuit.

We note that our results can also apply to the manipulation of the NHSE
induced by the on-site dissipations. Recent experiment Yi2021  and theoretical
proposals Yangzs ; Zhoulh  suggest realizing and detecting the NHSE in the
dissipative ultracold atom systems, where the gradient fields can be easily
realized Arimondo2008 ; Wolfgang2013 ; Ibloch . Therefore, the control effect
of the electric fields on the NHSE can be detected in optical lattices.
Furthermore, this control effect can also be detected based on the photonic
quantum walk Xiao2020  and the sideband cooling setups in trapped ion systems
ionRMP ; ionWang .

Summary and discussion.— We have investigated the control effect of electric
fields on NHSE analytically and numerically. For the pure dc field case, in
the thermodynamic limit, a weak Stark localization induced by the dc field is
sufficient to win the competition with NHSE, so a non-zero dc field can
suppress the NHSE. When the system size is finite, the new interesting
phenomenon of the size-dependent NHSE will emerge, because the number of skin
modes is size-independent when the size exceed a critical value. For the pure
ac field case, only the special field strength $E_{1}$ and frequency $\omega$
that satisfy $\mathcal{J}_{0}({E_{1}}/{\omega})=0$ can suppress the NHSE due
to the dynamic localization. For the mixed field case, if $E_{0}/\omega$ is
not an integer, the NHSE will be suppressed by the Stark localization induced
by the dc field. For the integer $E_{0}/\omega$, NHSE can exist, because the
particle can move toward the boundary by the photon absorption or emission
except for the special case with $E_{1}/\omega$ being one of the zeroes of
$\mathcal{J}_{E_{0}/\omega}$, which causes the dynamic localization.

The control effects of electric fields on the NHSE are abundant, and moreover,
electric fields can be easily applied to a system. Thus, the manipulation
methods can be widely used in experiments and the fabrication of new devices.
For instance, based on the phenomenon of the size-dependent NHSE and the
sensitivity of the NHSE versus the field strength and frequency of ac fields
near the special situations that satisfy
$\mathcal{J}_{0}({E_{1}}/{\omega})=0$, by detecting the signals on the
boundary of a non-Hermitian system, one can carry out accurate measurements of
electric fields, which are important for many critical applications in science
and industry. As a second example, the NHSE can be used to design some
devices, such as directional amplifiers Metelmann ; Metelmann2 ; Peterson2017
; Barzanjeh ; McDonald2018 ; Malz2018 ; Harris2019 ; Wentan2021  and light
funnels Weidemann2020 . The appearance or disappearance of the directional
amplification and the funneling effect can be manipulated by using dc or ac
fields, such as by changing the strength or frequency of an added ac field.
Thus, one can design the switch of these devices by using electric fields.

###### Acknowledgements.

We thank L. Li and H. Jiang for valuable discussions. This work was supported
by the National Natural Science Foundation of China (Grants No. U1801661,
No.12104205, No.12104210), the Key-Area Research and Development Program of
Guangdong Province (Grant No. 2018B030326001) ,Guangdong Provincial Key
Laboratory (Grant No.2019B121203002).

## References

  * (1) Y. Ashida, Z. Gong, and M. Ueda, Non-Hermitian physics, Advances in Physics 69, 249 (2020).
  * (2) C. M. Bender and S. Boettcher, Real spectra in non-Hermitian Hamiltonians having PT-symmetry, Phys. Rev. Lett. 80, 5243 (1998); C. M. Bender, D. C. Brody, and H. F. Jones, Complex extension of quantum mechanics, Phys. Rev. Lett. 89, 270401 (2002).
  * (3) N. Hatano and D. R. Nelson, Vortex pinning and non-Hermitian quantum mechanics, Phys. Rev. B 56, 8651 (1997); N. Hatano and D. R. Nelson, Non-Hermitian delocalization and eigenfunctions, Phys. Rev. B 58, 8384 (1998).
  * (4) T. E. Lee, Anomalous edge state in a non-Hermitian lattice, Phys. Rev. Lett. 116, 133903 (2016).
  * (5) H. Shen, B. Zhen, and L. Fu, Topological band theory for non-Hermitian Hamiltonians, Phys. Rev. Lett. 120, 146402 (2018).
  * (6) Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Topological phases of non-Hermitian systems, Phys. Rev. X 8, 031079 (2018).
  * (7) S. Yao and Z. Wang, Edge states and Topological invariants of non-Hermitian syatems, Phys. Rev. Lett. 121, 086803 (2018); S. Yao, F. Song, and Z. Wang, Non-Hermitian Chern bands, Phys. Rev. Lett. 121, 136802 (2018).
  * (8) F. K. Kunst, E. Edvardsson, J. C. Budich, and E. J. Bergholtz, Biothogonal bulk-boundary correspondence in non-Hermitian systems, Phys. Rev. Lett. 121, 026808 (2018).
  * (9) Z. Yang, K. Zhang, C. Fang, and J. Hu, Non-Hermitian bulk-boundary correspondence and auxiliary generalized Brillouin zone theory, Phys. Rev. Lett. 125, 226402 (2020).
  * (10) C. H. Lee and R. Thomale, Anatomy of skin modes and topology in non-Hermitian systems, Phys. Rev. B 99, 201103 (2019).
  * (11) K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Symmetry and topology in non-Hermitian physics, Phys. Rev. X 9, 041015 (2019).
  * (12) S. Longhi, Topological Phase Transition in non-Hermitian Quasicrystals, Phys. Rev. Lett. 122, 237601 (2019); S. Longhi, Phase transitions in a non-Hermitian Aubry-Andre-Harper model, Phys. Rev. B 103, 054203 (2021).
  * (13) H. Jiang, L.-J. Lang, C. Yang, S.-L. Zhu, and S. Chen, Interplay of non-Hermitian skin effects and Anderson localization in non-reciprocal quasiperiodic lattices, Phys. Rev. B 100, 054301 (2019).
  * (14) Y. Liu, X. Jiang, J. Cao, and S. Chen, Non-Hermitian mobility edges in one-dimensional quasicrystals with parity-time symmetry, Phys. Rev. B 101, 174205 (2020); Y. Liu, Y. Wang, Z. Zheng, and S. Chen, Exact non-Hermitian mobility edges in one-dimensional quasicrystal lattice with exponentially decaying hopping and its dual lattice, Phys. Rev. B 103, 134208 (2021).
  * (15) L. Zhou, Floquet engineering of topological localization transitions and mobility edges in one-dimensional non-Hermitian quasicrystals, Phys. Rev. Research, 3, 033184 (2021).
  * (16) E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Exceptional topology of non-Hermitian systems, Rev. Mod. Phys. 93, 015005 (2021).
  * (17) H. Schomerus, Topologically protected midgap states in complex photonic lattices, Opt. Lett. 38, 1912 (2013).
  * (18) J. Li, A. K. Harter, J. Liu, L. de Melo, Y. N. Joglekar, and L. Luo. Observation of parity-time symmetry breaking transitions in a dissipative Floquet system of ultracold atoms, Nat. Commun. 10, 855 (2019).
  * (19) L. Xiao, T. Deng, K. Wang, G. Zhu, Z. Wang, W. Yi, and P. Xue, Non-Hermitian bulk-boundary correspondence in quantum dynamics, Nat. Phys. 16, 761 (2020).
  * (20) S. Weidemann, M. Kremer, T. Helbig, T. Hofmann, A. Stegmaier, M. Greiter, R. Thomale, and A. Szameit, Topological funneling of light, Science 368, 311 (2020).
  * (21) T. Helbig, T. Hofmann, S. Imhof, M. Abdelghany, T. Kiessling, L. Molenkamp, C. Lee, A. Szameit, M. Greiter, and R. Thomale, Generalized bulk-boundary correspondence in non-Hermitian topolectrical circuits, Nat. Phys. 16, 747 (2020).
  * (22) X.-X. Zhang and M. Franz, Non-Hermitian exceptional Landau quantization in electric circuits, Phys. Rev. Lett. 124, 046401 (2020).
  * (23) W. Zhang, X. Ouyang, X. Huang, X. Wang, H. Zhang, Y. Yu, X. Chang, Y. Liu, D.-L. Deng, and L.-M. Duan, Observation of non-Hermitian topology with nonunitary dynamics of solid-state spins, Phys. Rev. Lett. 127, 090501 (2021).
  * (24) Q. Liang, D. Xie, Z. Dong, H. Li, H. Li, B. Gadway, W. Yi, and B. Yan, Observation of non-Hermitian skin effect and topology in ultracold atoms, arXiv:2201.09478.
  * (25) C. Dembowski, H.-D. Gr$\ddot{a}$f, H. L. Harney, A. Heine, W. D. Heiss, H. Rehfeld, and A. Richter, Experimental observation of the topological structure of exceptional points, Phys. Rev. Lett. 86, 787 (2001); C. Dembowski, B. Dietz, H.-D. Gr$\ddot{a}$f, H. L. Harney, A. Heine, W. D. Heiss, and A. Richter, Encircling an exceptional point, Phys. Rev. E 69, 056216 (2004).
  * (26) J. Wiersig, Enhancing the sensitivity of frequency and energy splitting detection by using exceptional points: application to microcavity sensors for single-particle detection, Phys. Rev. Lett. 112, 203901 (2014).
  * (27) W. Hu, H. Wang, P. P. Shum, Y. D. Chong. Exceptional points in a non-Hermitian topological pump. Phys. Rev. B 95, 184306 (2017).
  * (28) W. Chen, Ş. K. özdemir, G. Zhao, J. Wiersig, and L. Yang, Exceptional points enhance sensing in an optical microcavity, Nature (London) 548, 192 (2017).
  * (29) H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia, R. El-Ganainy, D. N. Christodoulides, and M. Khajavikhan, Enhanced sensitivity at higher-order exceptional points, Nature (London) 548, 187 (2017).
  * (30) M. Zhang, W. Sweeney, C. W. Hsu, L. Yang, A. D. Stone, and L. Jiang, Quantum noise theory of exceptional point amplifying sensors, Phys. Rev. Lett. 123, 180501 (2019).
  * (31) M.-A. Miri and A. Alù, Exceptional points in optics and photonics, Science 363, eaar7709 (2019).
  * (32) Y. Xu, S.-T. Wang, and L.-M. Duan, Weyl exceptional rings in a three-dimensional dissipative cold atomic Gas, Phys. Rev. Lett. 118, 045701 (2017).
  * (33) A. Cerjan, S. Huang, M. Wang, K. P. Chen, Y. Chong, and M. C. Rechtsman, Experimental realization of a Weyl exceptional ring, Nat. Photonics 13, 623 (2019).
  * (34) C. H. Lee, L. Li, and J. Gong, Hybrid higher-order skin-topological modes in nonreciprocal systems, Phys. Rev. Lett. 123, 016805 (2019).
  * (35) K. Kawabata, S. Higashikawa, Z. Gong, Y. Ashida, and M. Ueda, Topological unification of time-reversal and particle-hole symmetries in non-Hermitian physics, Nat. Commun. 10, 297 (2019).
  * (36) U. Magnea, Random matrices beyond the Cartan classification, Journal of Physics A: Mathematical and Theoretical 41, 045203 (2008).
  * (37) C. C. Wojcik, X.-Q. Sun, T. Bzdu$\check{s}$ek, and S. Fan, Homotopy characterization of non-Hermitian Hamiltonians, Phys. Rev. B 101, 205417 (2020).
  * (38) Z. Li and R. S. K. Mong, Homotopical classification of non-Hermitian band structures, Phys. Rev. B 103, 155129 (2021).
  * (39) H. Hu and E. Zhao, Knots and non-Hermitian bloch bands, Phys. Rev. Lett. 126, 010401 (2021).
  * (40) C.-H. Liu and S. Chen, Topological classification of defects in non-Hermitian systems, Phys. Rev. B 100,144106(2019); C.-H. Liu, H. Jiang, and S. Chen, Topological classification of non-Hermitian systems with reflection symmetry, Phys. Rev. B 99, 125103 (2019).
  * (41) V. M. M. Alvarez, J. E. B. Vargas, and L. E. F. F. Torres, Non-hermitian robust edge states in one dimension: Anomalous localization and eigenspace condensation at exceptional points, Phys. Rev. B 97, 121401 (2018).
  * (42) D. S. Borgnia, A. J. Kruchkov, R.-J. Slager, Non-Hermitian boundary modes and topology, Phys. Rev. Lett. 124, 056802 (2020).
  * (43) K. Zhang, Z. Yang, and C. Fang, Correspondence between winding numbers and skin modes in non-Hermitian systems, Phys. Rev. Lett. 125, 126402 (2020); K. Zhang, Z. Yang, and C. Fang, Universal non-Hermitian skin effect in two and higher dimensions, arXiv:2102.05059.
  * (44) Y. Yi and Z. Yang, Non-Hermitian skin modes induced by on-site dissipations and chiral tunneling effect, Phys. Rev. Lett. 125, 186802 (2020).
  * (45) N. Okuma, K. Kawabata, K. Shiozaki, and M. Sato, Topological origin of non-Hermitian skin effects, Phys. Rev. Lett. 124, 086801 (2020).
  * (46) L. Li, C. H. Lee, S. Mu, and J. Gong, Critical non-Hermitian skin effect, Nat. Commun. 11, 5491 (2020); L. Li, C. H. Lee, and J. Gong, Impurity induced scale-free localization, Commun. Phys. 4, 42 (2021).
  * (47) M. Lu, X.-X. Zhang, and M. Franz, Magnetic suppression of non-Hermitian skin effects, Phys. Rev. Lett. 127, 256402 (2021).
  * (48) K. Shao, H. Geng, W. Chen, and D. Y. Xing, Interplay between non-Hermitian skin effect and magnetic field: skin modes suppression, Onsager quantization and $MT$ phase transition, arXiv:2111.04412.
  * (49) X.-Q. Sun, P. Zhu, and T. L. Hughes, Geometric response and disclination-induced skin effects in non-Hermitian systems, Phys. Rev. Lett. 127, 066401 (2021).
  * (50) R. Hamazaki, K. Kawabata, and M. Ueda, Non-Hermitian many-body localization, Phys. Rev. Lett. 123, 090603 (2019).
  * (51) See Supplemental Material for details on (I) deriving the analytical solutions of quantum dynamics; (II) discussing the analytical results; (III) the finite size effect in the pure dc field case; (IV) dynamical behavior in the absence of external fields. The Supplemental Materials includes the references Bsel ; Paris1984 ; Olver2010 .
  * (52) G. N. Watson, A Treatise on the Theory of Bessel Functions, Cambridge Mathematical Library (Cambridge University Press, 1995).
  * (53) G. H. Wannier, Dynamics of band electrons in electric and magnetic fields, Rev. Mod. Phys. 34, 645 (1962).
  * (54) D. H. Dunlap and V. M. Kenkre, Dynamic localization of a charged particle moving under the influence of an electric field, Phys. Rev. B 34, 3625 (1986); D. H. Dunlap and V. M. Kenkre, Dynamic localization of a particle in an electric field viewed in momentum space: Connection with Bloch oscillations, Phys. Lett. A 127, 438 (1988).
  * (55) X.-G. Zhao, R. Jahnke, and Q. Niu, Dynamic fractional stark ladders in dc-ac fields, Physics Letters A 202, 297 (1995); X.-G. Zhao, G. A. Georgakis, and Q. Niu, Rabi oscillations between Bloch bands, Phys. Rev. B 54, 5235 (1996).
  * (56) S. Longhi, M. Marangoni, M. Lobino, R. Ramponi, P. Laporta, E. Cianci, and V. Foglietti, Observation of dynamic localization in periodically curved waveguide arrays, Phys. Rev. Lett. 96, 243901 (2006).
  * (57) A. Eckardt, M. Holthaus, H. Lignier, A. Zenesini, D. Ciampini, O. Morsch, and E. Arimondo, Exploring dynamic localization with a Bose-Einstein condensate, Phys. Rev. A 79, 013611 (2009).
  * (58) S. Guo, C. Dong, F. Zhang, J. Hu, and Z. Yang, Theoretical prediction of non-Hermitian skin effect in ultracold atom systems, arXiv:2111.04220.
  * (59) L. Zhou, H. Li, W. Yi, and X. Cui, Engineering non-Hermitian skin effect with band topology in ultracold gases, arXiv:2111.04196.
  * (60) C. Sias, H. Lignier, Y. P. Singh, A. Zenesini, D. Ciampini, O. Morsch, and E. Arimondo, Observation of photon-assisted tuneling in optical lattices, Phys. Rev. Lett. 100, 040404 (2008).
  * (61) H. Miyake, G. A. Siviloglou, C. J. Kennedy, W. C. Burton, and W. Ketterle, Realizing the Harper Hamiltonian with laser-assisted tunneling in optical lattices, Phys. Rev. Lett. 111, 185302 (2013).
  * (62) M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Realization of the Hofstadter Hamiltonian with ultracold atoms in optical lattices, Phys. Rev. Lett. 111, 185301 (2013).
  * (63) D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Quantum dynamics of single trapped ions, Rev. Mod. Phys. 75, 281 (2003).
  * (64) Z. Wang, Y. Peng, Y. Lu, Y. Wang, and J. Jie (unpublished).
  * (65) A. Metelmann and A. A. Clerk, Nonreciprocal photon transmission and amplification via reservoir engineering, Phys. Rev. X 5, 021025 (2015).
  * (66) K. Fang, J. Luo, A. Metelmann, M. H. Matheny, F. Marquardt, A. A. Clerk, and O. Painter, Generalized non-reciprocity in an optomechanical circuit via synthetic magnetism and reservoir engineering, Nature Phys. 13, 465 (2017).
  * (67) G. A. Peterson, F. Lecocq, K. Cicak, R. W. Simmonds, J. Aumentado, and J. D. Teufel, Demonstration of efficient nonreciprocity in a microwave optomechanical circuit, Phys. Rev. X 7, 031001 (2017).
  * (68) S. Barzanjeh, M. Aquilina, and A. Xuereb, Manipulating the flow of thermal noise in quantum devices, Phys. Rev. Lett. 120, 060601 (2018).
  * (69) A. McDonald, T. Pereg-Barnea, and A. A. Clerk, Phase-dependent chiral transport and effective non-hermitian dynamics in a bosonic Kitaev-Majorana chain, Phys. Rev. X 8, 041031 (2018).
  * (70) D. Malz, L. D. Tóth, N. R. Bernier, A. K. Feofanov, T. J. Kippenberg, and A. Nunnenkamp, Quantum-limited directional amplifiers with optomechanics, Phys. Rev. Lett. 120, 023601 (2018).
  * (71) H. Xu, L. Jiang, A. A. Clerk, and J. G. E. Harris, Nonreciprocal control and cooling of phonon modes in an optomechanical system, Nature 568, 65 (2019).
  * (72) W.-T. Xue, M.-R. Li, Y.-M. Hu, F. Song, and Z. Wang, Simple formulas of directional amplification from non-Bloch band theory, Phys. Rev. B 103, L241408 (2021).
  * (73) R. B. Paris, An Inequality for the Bessel Function $J_{\nu}(\nu{}x)$, SIAM J. Math. Anal., 15, 203 (1984).
  * (74) F. W. J. Olver, D. W. Lozier, and R. F. Boisvert, NIST handbook of mathematical functions, (Cambridge University Press, 2010).

## Supplementary Material:
Manipulating non-Hermitian skin effect via electric field

In the Supplementary Materials, we first give the details of deriving the
analytical solutions of quantum dynamics. Then, we discuss the analytical
results and obtain the effects of the pure dc field, pure ac field and the
mixed field on the non-Hermitian skin effect. Finally, we discuss the finite
size effect in the pure dc field case and the dynamical behavior in the
absence of external fields.

##  I I. Analytical solutions of quantum dynamics

We derive the analytical solutions of quantum dynamics for our one-dimensional
non-Hermitian model driven by an arbitrary time-dependent electric field
$E(t)$. The Hamitonian are following,

$\displaystyle\hat{H}(t)=\sum_{m=-\infty}^{+\infty}\left(J_{L}|m\rangle\langle
m+1|+J_{R}|m+1\rangle\langle
m|\right)+E(t)\sum_{m=-\infty}^{+\infty}m|m\rangle\langle m|,$ (S1)

where $J_{R}$ and $J_{L}$ are the strengths of the leftforward and
rightforward hopping, respectively. We expand the time-dependent quantum state
as $|\psi(t)\rangle=\sum_{m=-\infty}^{+\infty}C_{m}(t)|m\rangle$ and
substitute it into Schrödinger equation
$i\partial_{t}|\psi(t)\rangle=\hat{H}(t)|\psi(t)\rangle$ (we set $\hbar=1$
throughout this Supplementary Material) to arrive at the equation of the time-
dependent amplitudes $C_{m}(t)=\langle m|\psi(t)\rangle$,

$\displaystyle
i\partial_{t}C_{m}(t)=J_{L}C_{m+1}(t)+J_{R}C_{m-1}(t)+E(t)mC_{m}(t).$ (S2)

In solving Eq. (S2), we transfer to the momentum space firstly by the discrete
Fourier transformation,

$\displaystyle C_{k}(t)=\sum_{m=-\infty}^{+\infty}e^{-ikm}C_{m}(t),$ (S3)

and then we can directly rewrite Eq. (S2) as

$\displaystyle\partial_{t}C_{k}(t)-E(t)\partial_{k}C_{k}(t)=-i\left(J_{L}e^{ik}+J_{R}e^{-ik}\right)C_{k}(t).$
(S4)

This partial differential equation can be transferred to an ordinary
differential equation by the variables transformation

$\displaystyle p=t,q=k+\eta(t),$ (S5)

with $\eta(t)=\int_{0}^{t}E(t^{\prime})dt^{\prime}$. This gives rise to
$C_{k}(t)=C(p,\tau)$ with $\tau=q-\eta(p)$ and the chain rule of the
derivation of $C(p,\tau)$ over $p$ results in the ordinary differential
equation,

$\displaystyle\frac{dC(p,\tau)}{dp}$ $\displaystyle=$
$\displaystyle\frac{\partial C(p,\tau)}{\partial p}\frac{\partial p}{\partial
p}+\frac{\partial C(p,\tau)}{\partial\tau}\frac{\partial\tau}{\partial p},$
(S6) $\displaystyle=$ $\displaystyle\frac{\partial C(p,\tau)}{\partial
p}-E(p)\frac{\partial C(p,\tau)}{\partial\tau},$ $\displaystyle=$
$\displaystyle-i\left(J_{L}e^{i\tau}+J_{R}e^{-i\tau}\right)C(p,\tau),$

where we applied Eq. (S4) to the last step. Eq. (S6) can be solved by
integrating $p$ over two sides and then

$\displaystyle C(p,\tau)$ $\displaystyle=$ $\displaystyle
C(0,\tau|_{p=0})e^{-i\int_{0}^{p}\left[J_{L}e^{i\left(q-\eta(p^{\prime})\right)}+J_{R}e^{-i\left(q-\eta(p^{\prime})\right)}\right]dp^{\prime}},$
(S7) $\displaystyle=$ $\displaystyle
C(0,\tau|_{p=0})e^{-i\int_{0}^{p}\left[(J_{L}+J_{R})\cos\left(q-\eta(p^{\prime})\right)+i(J_{L}-J_{R})\sin\left(q-\eta(p^{\prime})\right)\right]dp^{\prime}},$
$\displaystyle=$ $\displaystyle
C(0,q)e^{-i\int_{0}^{p}\left[J_{+}\cos(q-\eta(p^{\prime}))+iJ_{-}\sin(q-\eta(p^{\prime}))\right]dp^{\prime}},$

where $J_{\pm}=J_{L}\pm J_{R}$. Transferring Eq. (S7) back to the form with
variables $(k,t)$ yields the solution of Eq. (S4),

$\displaystyle
C_{k}(t)=C_{k+\eta(t)}(0)e^{-i\int_{0}^{t}\left[J_{+}\cos(k+\eta(t)-\eta(t^{\prime}))+iJ_{-}\sin(k+\eta(t)-\eta(t^{\prime}))\right]dt^{\prime}}.$
(S8)

To simplify the solutions, we define the following functions

$\displaystyle\mathcal{U}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}\cos\left[\eta(t)-\eta(t^{\prime})\right]dt^{\prime},$
(S9) $\displaystyle\mathcal{V}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}\sin\left[\eta(t)-\eta(t^{\prime})\right]dt^{\prime},$
(S10)

and the Eq. (S8) arrives at

$\displaystyle C_{k}(t)$ $\displaystyle=$ $\displaystyle
C_{k+\eta(t)}(0)e^{-i\left[J_{+}(\mathcal{U}\cos k-\mathcal{V}\sin
k)+iJ_{-}(\mathcal{U}\sin k+\mathcal{V}\cos k))\right]dt^{\prime}},$ (S11)
$\displaystyle=$ $\displaystyle
C_{k+\eta(t)}(0)e^{i\left(\overline{\mathcal{V}}\sin
k-\overline{\mathcal{U}}\cos k\right)}.$

where

$\displaystyle\overline{\mathcal{U}}(t)=J_{+}\mathcal{U}(t)+iJ_{-}\mathcal{V}(t),\overline{\mathcal{V}}(t)=J_{+}\mathcal{V}(t)-iJ_{-}\mathcal{U}(t).$
(S12)

We apply the following expansions, which is expanded by ordinary Bessel
functions $\mathcal{J}_{m}(x)$, to Eq. (S11) BselS ,

$\displaystyle e^{-i\overline{\mathcal{U}}\cos k}$ $\displaystyle=$
$\displaystyle\sum_{m=-\infty}^{+\infty}e^{-i\frac{m\pi}{2}}e^{imk}\mathcal{J}_{m}(\overline{\mathcal{U}}),$
(S13) $\displaystyle e^{i\overline{\mathcal{V}}\sin k}$ $\displaystyle=$
$\displaystyle\sum_{m=-\infty}^{+\infty}e^{imk}\mathcal{J}_{m}(\overline{\mathcal{V}}),$
(S14)

and transfer the solutions in Eq. (S11) back to the spatial space by the
discrete Fourier transformation,

$\displaystyle C_{m}(t)=\sum_{k}e^{ikm}C_{k}(t).$ (S15)

Thus the time-dependent amplitude in spatial space $C_{m}(t)$ are

$\displaystyle C_{m}(t)$ $\displaystyle=$
$\displaystyle\sum_{k}e^{ikm}C_{k+\eta(t)}(0)e^{i\left(\overline{\mathcal{V}}\sin
k-\overline{\mathcal{U}}\cos k\right)},$ (S16) $\displaystyle=$
$\displaystyle\sum_{k}\sum_{n,r=-\infty}^{+\infty}e^{ik(m+n+r)}C_{k+\eta(t)}(0)e^{-i\frac{n\pi}{2}}J_{n}(\overline{\mathcal{U}})J_{r}(\overline{\mathcal{V}}),$
$\displaystyle=$
$\displaystyle\sum_{k}\sum_{n,r=-\infty}^{+\infty}e^{i(k+\eta(t))r}C_{k+\eta(t)}(0)e^{-i\eta(t)r}e^{-i\frac{n\pi}{2}}J_{n}(\overline{\mathcal{U}})J_{r-n-m}(\overline{\mathcal{V}}),$
$\displaystyle=$
$\displaystyle\sum_{n,r=-\infty}^{+\infty}C_{r}(0)e^{-i\eta(t)r}e^{-i\frac{n\pi}{2}}J_{n}(\overline{\mathcal{U}})J_{r-n-m}(\overline{\mathcal{V}}),$
$\displaystyle=$
$\displaystyle\sum_{n,r=-\infty}^{+\infty}C_{r}(0)e^{-i\eta(t)r}e^{-i\frac{n\pi}{2}}J_{n}(\overline{\mathcal{U}})e^{-i(n+m-r)\pi}J_{n+m-r}(\overline{\mathcal{V}}),$
$\displaystyle=$
$\displaystyle\sum_{r=-\infty}^{+\infty}(-1)^{m-r}C_{r}(0)e^{-i\eta(t)r}\left[\sum_{n=-\infty}^{+\infty}J_{n}(\overline{\mathcal{U}})J_{n+m-r}(\overline{\mathcal{V}})e^{i\frac{n\pi}{2}}\right],$
$\displaystyle=$
$\displaystyle\sum_{r=-\infty}^{+\infty}(-1)^{m-r}C_{r}(0)e^{-i\eta(t)r}\left(\frac{\overline{\mathcal{V}}+i\overline{\mathcal{U}}}{\overline{\mathcal{V}}-i\overline{\mathcal{U}}}\right)^{\frac{m-r}{2}}J_{m-r}\left(\sqrt{\overline{\mathcal{U}}^{2}+\overline{\mathcal{V}}^{2}}\right).$

In the last step, we applied the Graf’s addition theorem of the Bessel
functions BselS . Through the relations in Eq. (S12), we can straightforwardly
obtain that

$\displaystyle\frac{\overline{\mathcal{V}}+i\overline{\mathcal{U}}}{\overline{\mathcal{V}}-i\overline{\mathcal{U}}}=\frac{J_{R}}{J_{L}}\frac{i\mathcal{V}-\mathcal{U}}{i\mathcal{V}+\mathcal{U}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\overline{\mathcal{U}}^{2}+\overline{\mathcal{V}}^{2}=2\sqrt{J_{L}J_{R}(\mathcal{U}^{2}+\mathcal{V}^{2})}=2\sqrt{J_{L}J_{R}(u^{2}+v^{2})}.$
(S17)

Those expressions simplify the solution of Eq. (S16) to

$\displaystyle C_{m}(t)$ $\displaystyle=$
$\displaystyle\sum_{n=-\infty}^{+\infty}(-1)^{m-n}C_{n}(0)e^{-i\eta(t)n}J_{m-n}\left(2\sqrt{J_{L}J_{R}\left[\mathcal{U}^{2}(t)+\mathcal{V}^{2}(t)\right]}\right)\left[\frac{J_{R}}{J_{L}}\frac{i\mathcal{V}(t)-\mathcal{U}(t)}{i\mathcal{V}(t)+\mathcal{U}(t)}\right]^{\frac{m-n}{2}}.$
(S18)

This solution is valid for arbitrary initial state and for any time-dependent
driven $f(t)$, and also can be specified to the system with finite long
chains.

##  II II. Discussions of the analytical results

Base on the exact analytical solution of quantum evolution in Eq. (S18), we
can study the properties of the skin-effect and the electric fields induced
Wannier Stark localization by choosing an appropriate initial state. Here, we
set the system in site $m=0$ as the specific initial state, i.e.,
$C_{m=0}(t=0)=1$, this is equivalent to set any other sites $m\neq 0$ as the
initial state for the infinity long chain. Therefore, we can simplify the
solution in Eq. (S18) and explicitly express the probability as,

$\displaystyle\rho_{m}(t)=|C_{m}(t)|^{2}=\mathcal{J}_{m}^{2}\left(2\sqrt{J_{L}J_{R}\left[u^{2}(t)+v^{2}(t)\right]}\right)\left(\frac{J_{R}}{J_{L}}\right)^{m},$
(S19)

where we used the relation $\mathcal{U}^{2}+\mathcal{V}^{2}=u^{2}+v^{2}$ and
here

$\displaystyle
u(t)=\int_{0}^{t}dt^{\prime}\cos\eta(t^{\prime}),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ v(t)=\int_{0}^{t}dt^{\prime}\sin\eta(t^{\prime}).$ (S20)

We note that Eq. (S19) will become Eq.(5) of the main text when we choose the
initial state localized at $n_{0}$, i.e., $C_{m=n_{0}}(t=0)=1$ instead of
$C_{m=0}(t=0)=1$. The exact solutions in Eq. (S19) are valid when
$J_{L}J_{R}\neq 0$ and project to Hermitian case when $J_{L}=J_{R}$. In the
Non-Hermitian cases $J_{L}\neq J_{R}$, the coefficients $(J_{R}/J_{L})^{m}$
are responsible to the skin effect. Eq. (S19) is normalized for the Hermitian
case, and for non-Hermitian cases, the addition normalization is needed.

Next, we consider the following electric fields

$\displaystyle E(t)=E_{0}+E_{1}\cos\omega t,$ (S21)

where $E_{0}$ and $E_{1}$ are the strengths for the dc and ac parts,
respectively. $\omega$ is the ac driven frequency. Under this driven, we have

$\displaystyle\eta(t)=E_{0}t+\frac{E_{1}}{\omega}\sin\omega t,$ (S22)

and then

$\displaystyle
u(t)=\int_{0}^{t}dt^{\prime}\cos\left(E_{0}t^{\prime}+\frac{E_{1}}{\omega}\sin\omega
t^{\prime}\right),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
v(t)=\int_{0}^{t}dt^{\prime}\sin\left(E_{0}t^{\prime}+\frac{E_{1}}{\omega}\sin\omega
t^{\prime}\right).$ (S23)

$u(t)$ and $v(t)$ are not bounded functions only in the absence of electric
fields, namely $E(t)=0$. In this case, we have $u(t)=t$ and $v(t)=0$, yielding
$\rho_{m}(t)=\mathcal{J}_{m}^{2}\left(2t\sqrt{J_{L}J_{R}}\right)(J_{R}/J_{L})^{m}$,
which gives purely skin effect due to the argument $x$ in $\mathcal{J}_{m}(x)$
is linearly increasing to infinity and then $\mathcal{J}_{m}(x)$ approaches to
vanish, such that only the probabilities of the edge sites survive by the
coefficients $(J_{R}/J_{L})^{m}$. This purely skin effect is shown in Fig.
S1(b) through the quantum dynamic and as a comparison, the Hermitian case is
shown in Fig. S1(a).

Figure S1:  The quantum dynamics in the absence of fields for Hermitian
case$\leavevmode\nobreak\ (a):J_{L}=J_{R}=1.0$, and for the Non-Hermitian
case$\leavevmode\nobreak\ (b):J_{L}=0.8,J_{R}=1.0$.

###  II.1 DC electric field

For this time-independent driven case $E(t)=E_{0}$, we have

$\displaystyle u(t)=\frac{\sin E_{0}t}{E_{0}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ v(t)=\frac{1-\cos E_{0}t}{E_{0}},$ (S24)

and the probability is

$\displaystyle\rho_{m}(t)=\mathcal{J}_{m}^{2}\left(\frac{4\sqrt{J_{L}J_{R}}}{E_{0}}\sin\frac{E_{0}t}{2}\right)\left(\frac{J_{R}}{J_{L}}\right)^{m}.$
(S25)

From Eq. (S25), the argument of $\mathcal{J}_{m}(x)$ is bounded and
oscillating along time when $E_{0}\neq 0$. It means that the probability
amplitude $\rho_{0}$ will oscillate back to 1 at the time points $t^{*}=2\pi
N/E_{0}$ with $N=0,1,2,\cdots$, due to the fact that
$\mathcal{J}_{m}(0)=\delta_{m0}$. This Wannier Stark localization will compete
with skin effect under the weakly driven shown in Fig. S2(a), and dominate the
dynamics for the strongly driven shown in Fig. S2(b).

Figure S2:  The quantum dynamics for the dc fields with the weakly
driven$\leavevmode\nobreak\ (a):E_{0}=0.03$, and the strongly
driven$\leavevmode\nobreak\ (b):E_{0}=1.0$. Other parameters are
$J_{L}=0.8,\leavevmode\nobreak\ J_{R}=1.0,\leavevmode\nobreak\ \omega=1.0$.

###  II.2 AC electric field

For this time-dependent driven case $E(t)=E_{1}\cos\omega t$, we have

$\displaystyle
u(t)=\int_{0}^{t}dt^{\prime}\cos\left(\frac{E_{1}}{\omega}\sin\omega
t^{\prime}\right),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
v(t)=\int_{0}^{t}dt^{\prime}\sin\left(\frac{E_{1}}{\omega}\sin\omega
t^{\prime}\right).$ (S26)

At each periodic cycles, i.e., $t^{*}=2\pi N/\omega$ with $N=0,1,2,\cdots$, we
have

$\displaystyle
u(t^{*})=\frac{t^{*}}{\pi}\int_{0}^{\pi}d\tau\cos\left(\frac{E_{1}}{\omega}\sin\tau\right)=t^{*}\mathcal{J}_{0}(\frac{E_{1}}{\omega}),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ v(t)=0,$ (S27)

and such that the probability at time $t^{*}$ reads

$\displaystyle\rho_{m}(t^{*})=\mathcal{J}_{m}^{2}\left(2t^{*}\sqrt{J_{L}J_{R}}\mathcal{J}_{0}(\frac{E_{1}}{\omega})\right)\left(\frac{J_{R}}{J_{L}}\right)^{m},$
(S28)

where shows that the system is localized only when $E_{1}/\omega$ is the zeros
of $J_{0}$, otherwise the skin effect will dominate the dynamics in the long
time. To explicitly present this, we define two bounded functions as

$\displaystyle\overline{u}(t)=u(t)-t\mathcal{J}_{0}(\frac{E_{1}}{\omega}),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \overline{v}(t)=v(t)-0.$ (S29)

Then we rewrite the probability as

$\displaystyle\rho_{m}(t)$ $\displaystyle=$
$\displaystyle\mathcal{J}_{m}^{2}\left(2t\sqrt{J_{L}J_{R}\left[\mathcal{J}_{0}^{2}(\frac{E_{1}}{\omega})+g_{1}(t)\mathcal{J}_{0}(\frac{E_{1}}{\omega})+g_{2}(t)\right]}\right)\left(\frac{J_{R}}{J_{L}}\right)^{m},$
(S30)

where $g_{1}(t)=2\overline{u}(t)/t$ and
$g_{2}(t)=(\overline{u}(t)^{2}+\overline{v}(t)^{2})/t^{2}$ are vanished at
long time, and then we have

$\displaystyle\rho_{m}(t\gg\frac{2\pi}{\omega})$ $\displaystyle\approx$
$\displaystyle\mathcal{J}_{m}^{2}\left(2t\sqrt{J_{L}J_{R}}\mathcal{J}_{0}(\frac{E_{1}}{\omega})\right)\left(\frac{J_{R}}{J_{L}}\right)^{m}.$
(S31)

In Fig. S3, we take the ratio $E_{1}/\omega$ as two non-zeros of
$\mathcal{J}_{0}(x)$ (a,c) where show the skin effect and as two zeros of
$\mathcal{J}_{0}(x)$ in (b,d) where show the dynamical localization.

Figure S3: The quantum dynamics for the ac fields with the
ratios$\leavevmode\nobreak\ (a):E_{1}/\omega=0.1,\leavevmode\nobreak\
(b):0.2402\leavevmode\nobreak\ (\text{1st zero of\leavevmode\nobreak\
}\mathcal{J}_{0}),\leavevmode\nobreak\ (c):3.0,\leavevmode\nobreak\
(d):5.52\leavevmode\nobreak\ (\text{2nd zero of\leavevmode\nobreak\
}\mathcal{J}_{0}(x))$. Other parameters are $J_{L}=0.8,\leavevmode\nobreak\
J_{R}=1.0,\leavevmode\nobreak\ \omega=1.0$. Figure S4:  The quantum dynamics
for the dc+ac fields with the ratios$\leavevmode\nobreak\
(a):E_{0}/\omega=1.0,\leavevmode\nobreak\
E_{1}/\omega=2.0,\leavevmode\nobreak\
(b):E_{0}/\omega=1.0,\leavevmode\nobreak\
E_{1}/\omega=3.832\leavevmode\nobreak\ (\text{1st zero of\leavevmode\nobreak\
}\mathcal{J}_{1}(x)),\leavevmode\nobreak\
(c):E_{0}/\omega=2.324,\leavevmode\nobreak\
E_{1}/\omega=2.0,\leavevmode\nobreak\
(d):E_{0}/\omega=3.0,\leavevmode\nobreak\ E_{1}/\omega=3.0$. Other parameters
are $J_{L}=0.8,\leavevmode\nobreak\ J_{R}=1.0,\leavevmode\nobreak\
\omega=1.0$.

###  II.3 DC+AC electric field

The results above show that even weakly dc field will push the system to
dynamical localization and while the ac field supports the skin effect except
the ratio $E_{1}/\omega$ touches the zeros of Bessel function
$\mathcal{J}_{0}(x)$. The competition between ac and dc fields in regard of
localization will significantly modify the behaviors above.

For this more general case, we firstly rewrite Eq. (S23) as ($\theta=\omega
t$)

$\displaystyle u(\theta)$ $\displaystyle=$
$\displaystyle\frac{1}{\omega}\int_{0}^{\omega
t}d\theta^{\prime}\cos\left(\frac{E_{0}}{\omega}\theta^{\prime}+\frac{E_{1}}{\omega}\sin\theta^{\prime}\right),$
(S32) $\displaystyle=$ $\displaystyle\frac{1}{\omega}\int_{0}^{\omega
t}d\theta^{\prime}\left[\cos\left(\frac{E_{0}}{\omega}\theta^{\prime}\right)\cos\left(\frac{E_{1}}{\omega}\sin\theta^{\prime}\right)-\sin\left(\frac{E_{0}}{\omega}\theta^{\prime}\right)\sin\left(\frac{E_{1}}{\omega}\sin\theta^{\prime}\right)\right],$
$\displaystyle v(\theta)$ $\displaystyle=$
$\displaystyle\frac{1}{\omega}\int_{0}^{\omega
t}d\theta^{\prime}\sin\left(\frac{E_{0}}{\omega}\theta^{\prime}+\frac{E_{1}}{\omega}\sin\theta^{\prime}\right),$
(S33) $\displaystyle=$ $\displaystyle\frac{1}{\omega}\int_{0}^{\omega
t}d\theta^{\prime}\left[\sin\left(\frac{E_{0}}{\omega}\theta^{\prime}\right)\cos\left(\frac{E_{1}}{\omega}\sin\theta^{\prime}\right)+\cos\left(\frac{E_{0}}{\omega}\theta^{\prime}\right)\sin\left(\frac{E_{1}}{\omega}\sin\theta^{\prime}\right)\right],$

Applying the following relations BselS ,

$\displaystyle\cos(z\sin\theta)=\mathcal{J}_{0}(z)+2\sum_{k=1}^{+\infty}\mathcal{J}_{2k}(z)\cos(2k\theta),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\sin(z\sin\theta)=2\sum_{k=0}^{+\infty}\mathcal{J}_{2k+1}(z)\sin[(2k+1)\theta],$
(S34)

and doing the integrating, we arrive at

$\displaystyle u(t)$ $\displaystyle=$
$\displaystyle\frac{\sin(E_{0}t)}{E_{0}}\mathcal{J}_{0}\left(\frac{E_{1}}{\omega}\right)+\sum_{k=1}^{+\infty}(-1)^{k}\mathcal{J}_{k}\left(\frac{E_{1}}{\omega}\right)\left[\frac{\sin(E_{0}-k\omega)t}{E_{0}-k\omega}+\frac{\sin(E_{0}+k\omega)t}{E_{0}+k\omega}\right],$
(S35) $\displaystyle v(t)$ $\displaystyle=$
$\displaystyle\frac{1-\cos(E_{0}t)}{E_{0}}\mathcal{J}_{0}\left(\frac{E_{1}}{\omega}\right)+\frac{2}{\omega}\sum_{k=1}^{+\infty}\mathcal{J}_{2k}\left(\frac{E_{1}}{\omega}\right)\frac{\frac{E_{0}}{\omega}-\frac{E_{0}}{\omega}\cos(E_{0}t)\cos(2k\omega
t)-2k\sin(E_{0}t)\sin(2k\omega t)}{(\frac{E_{0}}{\omega})^{2}-4k^{2}}$ (S36)
$\displaystyle+\frac{2}{\omega}\sum_{k=0}^{+\infty}\mathcal{J}_{2k+1}\left(\frac{E_{1}}{\omega}\right)\frac{-(2k+1)+(2k+1)\cos(E_{0}t)\cos\left[(2k+1)\omega
t\right]+\frac{E_{0}}{\omega}\sin(E_{0}t)\sin\left((2k+1)\omega
t\right]}{(\frac{E_{0}}{\omega})^{2}-(2k+1)^{2}}.$

These expressions look very complicated and are difficult in doing the
calculations continuously. Even so, we still can do some analysis depends on
whether the ratio $E_{0}/\omega$ is an integer or not. If $E_{0}/\omega$ is
not an integer, all the terms in Eqs. (S35-S36) are bounded oscillatory
functions of time such that the system will localize around the initial
states. It means that the dc part of electric field dominate over the ac part
and then indicate breaking of the skin effect. In the case of provided
$E_{0}/\omega$ is an integer, function $v(t)$ in Eq. (S36) will not deduce any
non-bounded functions and will have similar contributions as in the integer
case. However, the functions $u(t)$ in Eq. (S35) takes (set
$E_{0}=\nu\omega,\leavevmode\nobreak\ \nu\in Z$)

$\displaystyle u(t)$ $\displaystyle=$
$\displaystyle(-1)^{\nu}\mathcal{J}_{\nu}\left(\frac{E_{1}}{\omega}\right)t+\frac{\sin(E_{0}t)}{E_{0}}\mathcal{J}_{0}\left(\frac{E_{1}}{\omega}\right)+\sum_{k=1,k\neq\nu}^{+\infty}(-1)^{k}\mathcal{J}_{k}\left(\frac{E_{1}}{\omega}\right)\left[\frac{\sin(E_{0}-k\omega)t}{E_{0}-k\omega}+\frac{\sin(E_{0}+k\omega)t}{E_{0}+k\omega}\right],$
(S37)

and then the first term is a non-bounded and linearly increasing function of
time, such that at long time, we have

$\displaystyle
u(t\gg\frac{2\pi}{\omega})\approx(-1)^{\nu}\mathcal{J}_{\nu}\left(\frac{E_{1}}{\omega}\right)t,$
(S38)

and the probability

$\displaystyle\rho_{m}(t\gg\frac{2\pi}{\omega})$ $\displaystyle\approx$
$\displaystyle\mathcal{J}_{m}^{2}\left(2t\sqrt{J_{L}J_{R}}\mathcal{J}_{\frac{E_{0}}{\omega}}(\frac{E_{1}}{\omega})\right)\left(\frac{J_{R}}{J_{L}}\right)^{m},$
(S39)

which is merely replacing the order "0" in Eq. (S31) by $E_{0}/\omega$.

As a conclusion, in this "dc+ac" case, the system will be featured in
dynamical localization for all non-integer ratio $E_{0}/\omega$ (Fig. S4 (c))
and will be dominated by skin effect for all $E_{1}/\omega$ except zeros of
$\mathcal{J}_{E_{0}/\omega}(x)$ when the ratio $E_{0}/\omega$ takes integer
(Fig. S4 (a, b, d)). We show the parameters zone for the skin effect by the
solid lines in Fig.3(a) in main text and all the rest zone are for dynamical
localization.

Figure S5:  Dynamical evolution of a electron started from the lattice center
for the system described by Eq. (S1) with (a) $E_{0}=0.005$, $E_{1}=0$, (b)
$E_{0}=0.5$, $E_{1}=0$, (c) $E_{0}=0$, $E_{1}/\omega=2.405$, which is the
first zero point of $\mathcal{J}_{0}(E_{1}/\omega)$, (d) $E_{0}=0$,
$E_{1}/\omega=6.1$, (e) $E_{0}/\omega=0.5$, $E_{1}/\omega=1.3$, (f)
$E_{0}/\omega=1$, $E_{1}/\omega=5.7$. Except $J_{R}$ and $J_{L}$, (a),(b) and
Fig.1 (e),(f) of the main text have same parameters, (c),(d) and Fig.2 (c),(d)
of the main text have same parameters, (e),(f) and Fig.3 (b),(d) of the main
text have same parameters. Here $J_{R}=J_{L}=1$.

###  II.4 Hermitian case

To understand the physical picture mentioned in the main text more
intuitively, we here provide the time evolution of the corresponding Hermitian
case, i.e., $J_{R}=J_{L}$, as shown in Fig. S5. Figs. S5 (a),(b) show the dc
case, whose parameters are same with Figs.1 (e),(f) of the main text. The
previous discussions about the dc case can be directly applied to the
$J_{R}=J_{L}$ case, as shown in Fig. S5 (b). Although both Fig. S5(b) and
Fig.1(f) of the main text show the Stark localization, there also exist the
differences. From Fig. S5 (b), the particle oscillate around the initial
position, but when $J_{R}>J_{L}$, the oscillation of the particle is at the
right-hand side of the initial position [Fig.1(f)]. Comparing Fig. S5(c),
Fig.S5(e) and Fig.2(c), Fig.3 (b) of the main text, we can also clearly see
the phenomenon that the oscillation center move to the right of the initial
position when changing $J_{L}=J_{R}$ to $J_{L}<J_{R}$.

When $E_{0}$ is decreased, as discussed below, the oscillation amplitude will
increase. For a finite system, when the oscillation amplitude is larger than
the system size, the particle can arrive at the boundary, as shown in Fig. S5
(b). When $J_{R}>J_{L}$, the particle will arrive at the boundary and then
stay there ever since [Fig.1(e)], which corresponds to the existence of the
NHSE. The competition between the oscillation amplitude and the system size
induce the phenomenon of the size-dependent NHSE. Comparing Fig. S5(d),
Fig.S5(f) and Fig.2(d), Fig.3 (d) of the main text, it can be seen that the
extended states become the skin modes when changing $J_{L}=J_{R}$ to
$J_{L}<J_{R}$.

##  III III. Finite size effect for the pure dc case

In this section, we investigate the finite size effect in the pure dc field
case. According to the inequality proven by Paris Paris  for Bessel function,
the probability $\rho_{m}(t)$ is bounded from above for big enough $m$,

$\rho_{m}(t)\leq\mathcal{J}^{2}_{m-n_{0}}(m-n_{0})f_{\chi}^{2(m-n_{0})}(x_{t}).$
(S40)

Here $f_{\chi}(x_{t})=x_{t}e^{1-\chi{}x_{t}}$, $\chi=\sqrt{J_{L}/J_{R}}<1$,
and $x_{t}=\frac{4J_{R}}{(m-n_{0})|E_{0}|}\left|\sin\frac{E_{0}t}{2}\right|$.
The inequality holds if $4\sqrt{J_{L}J_{R}}\leq(m-n_{0})|E_{0}|$.
$f_{\chi}(x_{t})$ is monotonically increasing and bounded
$0\leq{}f_{\chi}(x_{t})\leq 1/\chi$, for $0\leq{}x_{t}\leq 1/\chi$. One can
calculate the derivation of $f_{\chi}$ to verify it.

$f^{\prime}_{\chi}(x_{t})=e^{1-\chi{}x_{t}}-\chi{}x_{t}e^{1-\chi{}x_{t}}=e^{1-\chi{}x_{t}}\left(1-\chi{}x_{t}\right).$
(S41)

And $f_{\chi}(x_{t}=1)=e^{1-\chi}>1$ since $\chi<1$. There exists
$0\leq{}x_{\star}<1$ such that $f_{\chi}(x_{\star})=1$ and
$0\leq{}f_{\chi}(x_{t})<1$ if $0\leq{}x_{t}<x_{\star}$. Note that $x_{\star}$
depends only on $J_{L}$ and $J_{R}$, and satisfies
$x_{\star}e^{1-\sqrt{J_{L}/J_{R}}x_{\star}}=1$. We also know
$|\mathcal{J}_{n}(\bullet)|\leq 1$ Olver . Thus

$\displaystyle\rho_{m}(t)$ $\displaystyle\leq$
$\displaystyle\mathcal{J}^{2}_{m-n_{0}}(m-n_{0})f_{\chi}^{2(m-n_{0})}\left(\frac{4J_{R}}{(m-n_{0})|E_{0}|}\left|\sin\frac{E_{0}t}{2}\right|\right)$
(S42) $\displaystyle\leq$
$\displaystyle\mathcal{J}^{2}_{m-n_{0}}(m-n_{0})\left[\underbrace{f_{\chi}\left(\frac{4J_{R}}{(m-n_{0})|E_{0}|}\right)}_{<1}\right]^{2(m-n_{0})}$
$\displaystyle\leq$
$\displaystyle\left[\underbrace{f_{\chi}\left(\frac{4J_{R}}{(m-n_{0})|E_{0}|}\right)}_{<1}\right]^{2(m-n_{0})}.$

As $m$ increases, $f_{\chi}\left({4J_{R}}/{(m-n_{0})|E_{0}|}\right)$ would
decrease further. Thus the decay rate of $\rho_{m}$ is faster than any
exponential decay rate as $m$ increases. Hence $\rho_{m}(t)$ would be
insignificant if $4J_{R}<(m-n_{0})|E_{0}|x_{\star}$. It is well known that
$\mathcal{J}_{-n}(\bullet)=(-1)^{n}\mathcal{J}_{n}(\bullet)$. For two sites
$2n_{0}-m$ (left side $<n_{0}$) and $m$ (right side $>n_{0}$) with equal
distance to $n_{0}$

$\frac{\rho_{2n_{0}-m}(t)}{\rho_{m}(t)}\approx\frac{\mathcal{J}_{n_{0}-m}^{2}\left(2t\sqrt{J_{L}J_{R}}\mathcal{J}_{0}(\frac{E_{1}}{\omega})\right)\left(\frac{J_{R}}{J_{L}}\right)^{n_{0}-m}}{\mathcal{J}_{m-n_{0}}^{2}\left(2t\sqrt{J_{L}J_{R}}\mathcal{J}_{0}(\frac{E_{1}}{\omega})\right)\left(\frac{J_{R}}{J_{L}}\right)^{m-n_{0}}}=\left(\frac{J_{R}}{J_{L}}\right)^{-2(m-n_{0})}.$
(S43)

Adding the fact that $J_{R}>J_{L}$, $\rho_{2n_{0}-m}(t)$ would be even less
significant, if $m>n_{0}$ and $4J_{R}<(m-n_{0})|E_{0}|x_{\star}$. Therefore,
in an infinite lattice, an electron would oscillate around its initial
position $n_{0}$ within a range of $4J_{R}/|E_{0}|x_{\star}$.

##  IV IV. Dynamical behavior of a particle on non-reciprocal tight-binding
model with non-Hermitian skin effect

In this section, we take the case without external field as an example to
discuss the dynamical behavior of a particle. When there is no electric field,
Eq. (S40) and Eq. (S43) would still be valid except $x_{t}=2tJ_{R}/(m-n_{0})$.
Eq. (S43) tells us that the electron would always favor the right-hand side of
$n_{0}$. Further, from Eq. (S41) we know that $f_{\chi}(x_{t})$ increases from
0 to $f^{\mathrm{max}}_{\chi}=\sqrt{J_{R}/J_{L}}$ as $x_{t}$ increases from 0
to $1/\chi$. Hence $f_{\chi}(x_{t})<1$ and thus the upper bound of $\rho_{m}$
(see Eq. (S40)) would be very small, when $x_{t}<x_{\star}$. Here
$f_{\chi}(x_{\star})=1$ and $x_{\star}<1$ as mentioned in the previous
section. Hence $\rho_{m}$ would be significant only if $x_{t}{\geq}x_{\star}$,
namely

$t\geq\frac{(m-n_{0})x_{\star}}{2J_{R}}.$ (S44)

Therefore, the electron reaches site $m$ around
$t_{a}=\frac{(m-n_{0})x_{\star}}{2J_{R}}$ and continues hopping rightward. Now
we consider the stable probability in a long time after the particle passing.
As discussed above, in the absence of external fields,
$\rho_{m}(t)=\mathcal{J}_{m-n_{0}}^{2}\left(2t\sqrt{J_{L}J_{R}}\right)(J_{R}/J_{L})^{m-n_{0}}$,
where $2t\sqrt{J_{L}J_{R}}$ is linearly increasing to infinity with increasing
time, and thus $\mathcal{J}_{m-n_{0}}(2t\sqrt{J_{L}J_{R}})$ tends to 0. We
consider $\rho_{m}/\rho_{m-1}$ when $t\rightarrow\infty$. For convenience, we
consider the long time averaged behavior after the particle has passed the
position $m$ for a long time. We notice that $\frac{\int^{\infty}_{t\gg
t_{a}}dt\mathcal{J}_{m-n_{0}}^{2}\left(2t\sqrt{J_{L}J_{R}}\right)}{\int^{\infty}_{t\gg
t_{a}}dt\mathcal{J}_{m-1-n_{0}}^{2}\left(2t\sqrt{J_{L}J_{R}}\right)}=1$ and so
$\rho_{m}/\rho_{m-1}=J_{R}/J_{L}$. Thus, for a system with non-Hermitian skin
effect, the probability distributions of the region that the particle has
passed through are not ergodic, which is obviously different from the system
without skin effect.

## References

  * (1) G. N. Watson, A Treatise on the Theory of Bessel Functions, Cambridge Mathematical Library (Cambridge University Press, 1995).
  * (2) R. B. Paris, An Inequality for the Bessel Function $J_{\nu}(\nu{}x)$, SIAM J. Math. Anal., 15, 203 (1984).
  * (3) F. W. J. Olver, D. W. Lozier, and R. F. Boisvert, NIST handbook of mathematical functions, (Cambridge University Press, 2010).