# Cyber-Physical Testbed Integrating RTAC with RTDS for Game-Theoretic
Topology Control Under Load Altering Attacks

Alaa Selim, _Student Member, IEEE_ , Junbo Zhao, _Senior Member, IEEE_ A.
Selim and J. Zhao are with the Department of Electrical and Computer
Engineering, University of Connecticut (e-mail<EMAIL_ADDRESS>junbo@uconn.edu).

###### Abstract

This paper introduces a cyber-physical testbed that integrates the Real-Time
Digital Simulator (RTDS) with the Real-Time Automation Controller (RTAC) to
enhance cybersecurity in electrical distribution networks. Focused on
addressing vulnerabilities to cyber attacks, our testbed employs an advanced
control algorithm developed in Python and executed through real-time
controllers. Key to our approach is the seamless integration of the host PC
machine, RTDS, and RTAC via the Modbus protocol. We present a game theory-
based topology control strategy as an effective response to cyber attacks,
specifically targeting smart meters. This testbed validates the efficacy of
our method in fully eliminating voltage violations due to load altering
attacks, showcasing substantial advancements in smart grid cybersecurity via
the innovative use of RTDS and RTAC to simulate and counteract complex cyber
threats.

###### Index Terms:

Cybersecurity, Cyber-Physical Systems, Game Theory, Cyber Attack Mitigation,
Real-Time Digital Simulator (RTDS), Real-Time Automation Controller (RTAC)

##  I Introduction

With the advent of smart grid technology, the integration of cyber systems
with physical grid operations has dramatically transformed the landscape of
the power sector. This fusion has ushered in an era of enhanced efficiency and
sophisticated management capabilities, as evident in the widespread adoption
of advanced metering infrastructure (AMI) and intelligent control systems. As
the sector continues to evolve amidst increasing digitalization, the focus
shifts towards developing resilient systems to not only withstand cyber
threats but also recover swiftly from them. For instance, integrating utility
devices like switches and capacitors introduces both complexity and an
opportunity for increased grid resilience. By leveraging the synergy between
cyber advancements, such as real-time data analytics, and physical grid
enhancements, such as adaptive protection schemes, we can innovate and improve
grid operations.

The deployment of advanced control algorithms through real-time controllers
stands as a cornerstone for validating the effectiveness and feasibility of
resiliency solutions. This paper specifically concentrates on harnessing the
versatility and power of Python, a high-level programming language known for
its efficiency and readability, to address critical power system problems by
enabling the use of advanced machine learning and model-based algorithms [1].
We aim to bridge the gap between theoretical algorithm development and
practical application by integrating Python with real-time simulation
environments, as shown in Fig. 1.

Figure 1: Conceptual model

Several methodologies are available to achieve seamless integration between
Python and real-time controllers, each with its own characteristics and
suitability. While API interfaces and socket communication are common methods,
they present certain limitations in the context of real-time operations. API
interfaces [2], for instance, often lack the immediacy required for real-time
data exchange, as they typically operate over HTTP protocols not optimized for
low latency. Socket communication [3], on the other hand, while potentially
faster and more direct, introduces complexities in data packet management and
necessitates a more in-depth understanding of network protocols. In contrast
to these methods, Modbus protocol [4] emerges as a more suitable candidate for
our integration needs. Despite being slower than socket communication in raw
speed, Modbus compensates with its remarkable stability, standardized
communication format, and ease of implementation. Its simplicity and
robustness, essential for industrial applications, outweigh the speed
advantage of socket communication, especially considering the latter’s
increased complexity and potential for errors in data handling.

TABLE I: Comparison of Cyber-Physical System Papers in the literature Criteria |  [5] |  [6] |  [7] |  [8] |  [9] |  [10] |  Our Study
---|---|---|---|---|---|---|---
Use of RTDS |  |  |  $\checkmark$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$
Python Integration with PC |  |  |  |  |  |  |  $\checkmark$
Hardware-In-the-Loop |  $\checkmark$ |  |  $\checkmark$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$
Advanced Algorithms |  |  |  $\checkmark$ |  |  $\checkmark$ |  $\checkmark$ |  $\checkmark$
Use of RTAC |  $\checkmark$ |  $\checkmark$ |  |  $\checkmark$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$
Topology control |  $\checkmark$ |  |  |  |  $\checkmark$ |  $\checkmark$ |  $\checkmark$
Use of industrial protocol |  $\checkmark$ |  $\checkmark$ |  |  $\checkmark$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$
Cyber attack focus |  |  |  |  |  |  |  $\checkmark$

With the successful integration of Python with real-time controllers such as
Real-Time Automation Controller (RTAC) and Real-Time Digital Simulator (RTDS)
via Modbus, we are now strategically positioned to tackle pressing challenges
in power systems, particularly cyber attacks targeting smart meters. Such
attacks pose a significant threat, as they can lead to system instability and
violations of operational standards, undermining the reliability and security
of the power grid. The sophistication of these cyber attacks requires an
equally advanced response – a task well-suited to Python’s capabilities. By
deploying these Python-developed algorithms on real-time controllers, we can
effectively monitor smart meter operations, detect anomalous behaviors
indicative of cyber attacks, and initiate immediate countermeasures to
maintain system stability and compliance with regulatory standards.

The landscape of smart grid cybersecurity is continually evolving, as
highlighted by the comprehensive analysis in [11], which focuses on the
complexities of man-in-the-middle attacks altering network topology data. In
the context of this expanding field, our research emerges as a crucial next
step, concentrating specifically on the cybersecurity of grid meters. Our work
diverges from the traditional approach of network data manipulation, instead
directing attention toward the vulnerabilities inherent in metering systems
and the application of topology control as an innovative mitigation technique.
This shift in focus is driven by the need to address the unique challenges
posed by meter-centric cyber threats, which have not been extensively explored
in existing literature.

In the realm of cyber attack mitigation algorithms, numerous sophisticated
methods have been explored, including Deep Reinforcement Learning (DRL) and
Model Predictive Control (MPC) as in [12]. Recognizing the diverse landscape
of these approaches, this paper introduces a novel perspective by implementing
a game theory-based approach. Our method focuses on identifying the optimal
topology changes to mitigate cyber threats effectively. By analyzing and
predicting the outcomes of various scenarios, our algorithm strategically
reconfigures grid topology, offering a proactive stance in safeguarding the
power system’s integrity against cyber vulnerabilities.

In this study, we deploy an advanced cyber-physical testbed to systematically
address the complexities of cyber attacks in power systems, leveraging a game
theory-based mitigation approach. The efficacy and sophistication of our
testbed are critically assessed through a detailed comparative analysis with
contemporary cyber-physical testbeds, as delineated in Table I. Our testbed
excels in ease of use and efficient integration of PC-implemented advanced
Python algorithms, making it highly effective for cyber attack mitigation in
power systems. The main contributions of this paper can be summarized as
follows:

  1. 1.

Advanced Cyber-Physical Testbed: The study effectively utilizes a state-of-
the-art cyber-physical testbed, showcasing its capabilities in simulating and
evaluating game theory-based algorithms against cyber threats.

  2. 2.

Novel Cyber Attack Scenario: Introduction of a novel and elusive cyber attack
on smart meters, designed to challenge detection methods and cause significant
system voltage violations, underscoring the complexity of modern cybersecurity
threats.

  3. 3.

Game Theory-based Topology Control Strategy: Proposal of a strategic,
minimalistic topology control method, leveraging game theory for efficient
cyber attack mitigation with reduced switching actions, demonstrating an
effective approach to maintaining system stability.

##  II Designing a cyber Attack Vector for Voltage Violations

Informed by insights from [13], we develop the cyber attack model based on a
targeted false data injection attack scenario. These attacks are recognized
for their stealth, targeting specific state variables to manipulate the power
grid’s smart meters. The $A\textit{v}$ vector in our model exemplifies this
approach by strategically manipulating phase-wise loads to induce maximum
voltage violations. This underscores the importance of enhancing security
measures in power systems to mitigate such stealthy and targeted attacks. This
vector targets the three phases of the system, represented as:

$A\textit{v}=\begin{pmatrix}A\textit{v}_{A}\\\ A\textit{v}_{B}\\\
A\textit{v}_{C}\end{pmatrix}$ (1)

where $A\textit{v}_{A}$, $A\textit{v}_{B}$, and $A\textit{v}_{C}$ correspond
to phase-wise load alterations. The attack is executed through a controlled
manipulation of a selected subset of smart meters, as shown in Fig. 2. The
methodology involves:

  * •

Overloading Phase A: Increasing load in phase A, denoted as $A\textit{v}_{A}$,
to create a high-load scenario.

  * •

Underloading Phases B and C: Concurrently reducing loads in phases B
($A\textit{v}_{B}$) and C ($A\textit{v}_{C}$) to exacerbate system imbalances
and vice versa.

  * •

Adaptive Monitoring: Continuous observation of network responses, followed by
these adjustments in $A\textit{v}_{A}$, $A\textit{v}_{B}$, and
$A\textit{v}_{C}$ for probable cyber attack scenario.

Figure 2: Attack vector

The aim is to quickly induce voltage violations by exploiting specific
vulnerabilities in the network. The goal of the attack vector is to
precipitate the maximum voltage violations with minimal interventions. Under
constraints, ensure minimal changes in load settings and adherence to network
operational standards. This attack model can be represented by an optimization
problem as follows:

$\displaystyle\max_{\\{A\textit{v}_{A,t},A\textit{v}_{B,t},A\textit{v}_{C,t}\\}_{t=1}^{T}}V_{vio,t}$
(2) $\displaystyle\text{s.t.
}\sum_{t=1}^{T}(A\textit{v}_{A,t}+A\textit{v}_{B,t}+A\textit{v}_{C,t})\leq
A\textit{v}_{max}$ (3) $\displaystyle
A\textit{v}_{A,t}=\begin{cases}A\textit{v}_{A,t-1}\pm\gamma_{A},&\text{if
}V_{vio,t}<n_{tar}\\\ A\textit{v}_{A,t-1},&\text{if }V_{vio,t}\geq
n_{tar}\end{cases}$ (4) $\displaystyle
A\textit{v}_{B,t}=\begin{cases}A\textit{v}_{B,t-1}\pm\gamma_{B},&\text{if
}V_{vio,t}<n_{tar}\\\ A\textit{v}_{B,t-1},&\text{if }V_{vio,t}\geq
n_{tar}\end{cases}$ (5) $\displaystyle
A\textit{v}_{C,t}=\begin{cases}A\textit{v}_{C,t-1}\pm\gamma_{C},&\text{if
}V_{vio,t}<n_{tar}\\\ A\textit{v}_{C,t-1},&\text{if }V_{vio,t}\geq
n_{tar}\end{cases}$ (6)
$\displaystyle\gamma(V_{vio},n_{tar})=\begin{cases}\alpha,&\text{if
}V_{vio}=0\\\ \beta(V_{vio}),&\text{if }0<V_{vio}<n_{tar}\end{cases}$ (7)
$\displaystyle\beta(V_{vio})=k-\delta\times V_{vio},\text{ where }k,\delta>0$
(8) $\displaystyle A\textit{v}_{A,t}\in[P_{D_{A}}\cdot n_{min},P_{D_{A}}\cdot
n_{max}],\quad\forall t$ (9) $\displaystyle
A\textit{v}_{B,t}\in[P_{D_{B}}\cdot n_{min},P_{D_{B}}\cdot
n_{max}],\quad\forall t$ (10) $\displaystyle
A\textit{v}_{C,t}\in[P_{D_{C}}\cdot n_{min},P_{D_{C}}\cdot
n_{max}],\quad\forall t$ (11) $\displaystyle n_{min}>0$ (12) $\displaystyle
n_{max}=f(V_{vio}-V_{vio_{init}})$ (13)

where:

  * •

$A\textit{v}_{A,t},A\textit{v}_{B,t},A\textit{v}_{C,t}$: Attack vectors on
phases A, B, and C at time $t$.

  * •

$V_{vio,t}$: Number of voltage violations defined per [14] standard at time
$t$.

  * •

$A\textit{v}_{max}$: Maximum allowed sum of attack vectors across phases.

  * •

$\gamma_{A},\gamma_{B},\gamma_{C}$: Incremental change values for phases A, B,
and C.

  * •

$n_{tar}$: Target number of voltage violations.

  * •

$\gamma(V_{vio},n_{tar})$: Function defining change parameter based on current
violations and target.

  * •

$\alpha,\beta(V_{vio})$: Parameters in the $\gamma$ function, where $\beta$ is
a function of current violations.

  * •

$k,\delta$: Constants in the $\beta$ function.

  * •

$P_{D_{A}},P_{D_{B}},P_{D_{C}}$: Normal load settings for phases A, B, and C.

  * •

$n_{min},n_{max}$: Minimum and maximum multipliers for load settings.

  * •

$V_{vio_{init}}$: Initial number of voltage violations.

  * •

$f(V_{vio}-V_{vio_{init}})$: Function relating $n_{max}$ to the change in
voltage violations.

The cyber attack in our study mimics a voltage unbalance fault, subtly
designed to evade detection by standard relays. By causing a slight voltage
unbalance, it avoids triggering typical monitoring systems, highlighting the
need for more sensitive detection methods to identify such stealthy threats.
For this reason, it is essential to calculate the voltage unbalance ($U_{t}$)
for each attack to ascertain that it remains undetectable, thereby ensuring
the attack’s stealthiness.

$U_{t}=\frac{\max\left(\begin{bmatrix}|V_{a}-V_{avg}|\\\ |V_{b}-V_{avg}|\\\
|V_{c}-V_{avg}|\end{bmatrix}\right)}{V_{avg}}\times 100$ (14)

where $V_{a},V_{b},\text{ and }V_{c}$ are the voltages of phases A, B, and C,
respectively, and $V_{avg}=\frac{V_{a}+V_{b}+V_{c}}{3}$. As per ANSI C84.1
[14], voltage unbalance exceeding 3% in distribution networks can lead to
inefficient operation and potential damage to electrical equipment. The attack
vector tends to stay within allowable voltage imbalance limits, making it less
easily detectable, as shown in Algorithm 1.

Algorithm 1 Cyber Attack Vector for load nodes

1:Initialize network parameters and set initial loads
$A\textit{v}_{A,0},A\textit{v}_{B,0},A\textit{v}_{C,0}$

2:for each time step $t=1$ to $T$ do

3: Adjust load on Phase A:
$A\textit{v}_{A,t}=A\textit{v}_{A,t-1}\pm\gamma_{A}$

4: Adjust loads on Phase B:
$A\textit{v}_{B,t}=A\textit{v}_{B,t-1}\mp\gamma_{B}$

5: Adjust loads on Phase C:
$A\textit{v}_{C,t}=A\textit{v}_{C,t-1}\mp\gamma_{C}$

6: Compute node voltages: $V_{nodes,t}$

7: Check voltage unbalance factor $U_{t}$ at time $t$

8: Count voltage violations at time $t$

9:end for

10:Analyze overall attack impact based on voltage violations and unbalance
factors

##  III Game Theory for Attack Mitigation

In response to the cyber attacks initiated in section II, particularly those
targeting network stability and voltage profiles, we propose a game-theoretic
approach shown in Algorithm 2 to control network topology. This method
involves adjusting switch states within the network to maintain stability,
ensure service continuity, and minimize switching actions for practical real-
time applications.

The network consists of a set of switches $S=\\{s_{1},s_{2},\ldots,s_{n}\\}$.
Each switch $s_{i}$ can be in one of two states, represented by the strategy
set $S_{i}=\\{\text{open},\text{closed}\\}$. The current state of the network
is denoted by $V$, reflecting the voltage profile post-attack as determined by
Algorithm 1.

The key objective is to modify the switch states in a way that ensures the
radiality of the network, denoted by $R(S)$. Radiality is a constraint that
preserves the tree structure of the network, which is essential for minimizing
disruptions and maintaining operational integrity. Additionally, we aim to
minimize the number of switching actions, which is crucial for the feasibility
of real-time operations. To achieve this, we introduce a cost function $C(S)$
associated with switch operations.

The payoff function $U(S_{i})$ is defined as $f(R(S),V(S),C(S))$, which
evaluates the effectiveness of each strategy in terms of maintaining
radiality, reducing voltage violations, and minimizing the cost associated
with switch operations. The optimal switch state $s_{i}^{*}$ maximizes the
payoff, given by:

$s_{i}^{*}=\arg\max U(S_{i})$

Subject to the constraints of radiality, voltage stability, and minimal
switching cost. Given a set of nodes in an electrical network, let $V_{i}$ be
the voltage at node $i$, and $S_{ij}$ represent the state of the switch
between nodes $i$ and $j$ (1 if closed, 0 if open). The radiality of the
network can be mathematically expressed as:

$R=\prod_{i,j\in\text{Nodes}}\left(1-\delta(V_{i},0)\delta(V_{j},0)(1-S_{ij})\right)$

where $\delta(x,y)$ is the Kronecker delta function. This function equals 1 if
$x=y$ and 0 otherwise. The radiality metric $R$ will be 1 if the network
maintains radiality, and 0 if it doesn’t. The switching cost function $C(S)$
is designed to penalize frequent switching, thereby encouraging solutions that
maintain stability with minimal switch operations.

Algorithm 2 Game-Theoretic Topology Control for Attack Mitigation with
Minimized Switching

1:Start: Voltage profile from Algorithm 1 after attack

2:Initialize: Switch states $S=\\{s_{1},s_{2},\ldots,s_{n}\\}$

3:while true do

4: Observe current $V$ and network topology

5: Ensure radiality $R(S)$

6: for each switch $s_{i}$ in $S$ do

7: Define strategies $S_{i}=\\{\text{open},\text{closed}\\}$

8: Calculate payoff $U(S_{i})=f(R(S),V(S),C(S))$

9: $s_{i}^{*}\leftarrow\arg\max U(S_{i})$

10: end for

11: Update switch states $S=\\{s_{1}^{*},s_{2}^{*},\ldots,s_{n}^{*}\\}$ with
minimal changes

12: Apply $S$ to the network

13: Sleep for control interval

14:end while

##  IV Cyber-Physical Testbed Setup

We have developed a comprehensive cyber-physical testbed designed to implement
our game-theoretic algorithm for mitigating cyber attacks. This setup involves
executing Python code on a PC machine, with subsequent communication to RTDS
through a GTNETx2 card. The integration of the PC machine with RTDS is
achieved via two distinct approaches.

###  IV-1 Direct Python-to-RTDS Integration

The first approach involves directly running Python scripts on the RTDS from
the host machine, as shown in Fig. 3. This method is particularly suitable for
research purposes due to its straightforward integration. Utilizing the
Pymodbus [15] library, the host machine functions as a client, sending control
commands to the RTDS. In this setup, the GTNETx2 card acts as a server,
receiving control actions to operate network switches and monitor voltage
levels within the network. Data on these voltage levels is then relayed back
to the host machine.

Figure 3: Network diagram for research purpose

###  IV-2 RTAC-Based Integration for Industrial Application

The second approach shown in Fig. 4 represents a more industrial application.
In this method, RTAC serves as a server to the host machine. The host machine
communicates control actions to the RTAC, which then relays this data to the
RTDS, acting as a client. This communication is facilitated through a
Software-Defined Networking (SDN) infrastructure and network switches. This
approach demonstrates the practical application of RTAC in the field, enabling
the execution of advanced control actions in real-world scenarios.

Figure 4: Network diagram for deployment in field

Both methods provide valuable insights into the practical implementation of
our algorithm. While the first approach offers ease of use and quick
integration for research and development, the second approach showcases the
potential for deploying such advanced control strategies in an industrial
setting, leveraging RTAC’s capabilities.

SDN in [16] significantly enhances network scalability, centralizes control,
improves resilience, and increases visibility, making it ideal for managing
extensive control actions. The scalability of SDN allows efficient traffic
management and inspection, which is crucial for large-scale networks. Its
centralized control architecture enables dynamic and efficient packet routing,
which is essential for handling numerous control actions. These features
collectively make SDN an excellent choice for complex network management
demands. While integrating Modbus with the RTDS provides numerous advantages,
it also presents limitations that must be addressed for effective
implementation. Firstly, Modbus supports a maximum of 125 holding registers
for reading data. This restriction can be particularly challenging in large-
scale applications, as it limits the amount of data retrieved in a single
query, potentially leading to incomplete data acquisition. To circumvent this
limitation, we have implemented a systematic approach to selectively retrieve
the most crucial data, ensuring that critical information is noticed despite
the limited number of registers. Another limitation is that Modbus in RTDS
natively supports only 16-bit unsigned integers. This constraint becomes
apparent when reading voltage levels, typically represented as floating-point
numbers. We employed a technique based on the IEEE 754 Standard for Floating-
Point Arithmetic to resolve this. According to this standard, a 32-bit single-
precision floating-point number (FLOAT32) can be effectively represented using
two consecutive 16-bit registers in the Modbus data model. This mapping allows
for the accurate reading of voltage levels as exact float values. Notably,
this approach is also recommended by the RTDS support team, aligning with
industry best practices for data accuracy and representation in digital
systems. This testbed configuration allows for accurate sensor measurements in
large-scale applications, ensuring reliable data analysis and control
strategies.

##  V Numerical Testing Analysis

In our experiment shown in Fig. 5, we conduct tests on the IEEE-123 bus
system, implemented using RSCAD version FX2.1 (see Fig. 6), where the system
starts typically with zero voltage violations. We employ the Pymodbus library
for communication between the host machine and the RTDS and RTAC through
Modbus/TCP. The ACSELERATOR RTAC software facilitates the configuration of the
RTAC, including mapping of data points and verifying connectivity. Our setup
begins with sending control actions to the switches, utilizing the two
approaches discussed in the previous section. For data collection, we assume
visibility over 206 smart meters in our test setup, as shown in Fig. 7.
However, due to Modbus’s limitation of 125 tags in RTDS, we develop a script
in RSCAD that prints all voltage readings as soon as the host machine triggers
any control action. These voltage data are then saved on the same host
machine. This setup allows for the seamless integration of real-time voltage
data into our Python-based mitigation algorithm, enabling it to respond
dynamically to changes in the network’s state and effectively mitigate any
identified voltage violations. Table II succinctly categorizes the Modbus
communication within a cyber-physical system. It details three key components:
’Voltage Data’ with Holding Registers 1-206 for individual node voltages,
’Topology Control’ using Coils 1-6 for switch operations (S1-S6), and ’Load
Control’ through Holding Registers 207-215 for various setpoint adjustments.
This structure demonstrates the protocol’s integral role in managing essential
aspects of voltage monitoring, network topology, and load management.

Figure 5: Cyber-Physical testbed setup TABLE II: Modbus Data Exchange with RTAC Modbus Entity | Type | Description
---|---|---
Holding Register 1 | Voltage Data | Node 1 Voltage
Holding Register 2 | Voltage Data | Node 2 Voltage
$\vdots$ | $\vdots$ | $\vdots$
Holding Register 206 | Voltage Data | Node 206 Voltage
Coil 1 | Topology Control | Switch Control S1
Coil 2 | Topology Control | Switch Control S2
$\vdots$ | $\vdots$ | $\vdots$
Coil 6 | Topology Control | Switch Control S6
Holding Register 207 | Load Control | Setpoint 1
$\vdots$ | $\vdots$ | $\vdots$
Holding Register 215 | Load Control | Setpoint 9

We initiate the cyber attack as outlined in Algorithm 1, strategically
implementing the attack vectors to analyze their impact on the IEEE-123 bus
system, starting from zero voltage violations level. The results of our cyber
attack simulation, as detailed in Table IV, exhibit the nuanced effects of
different attack vectors on the power network’s voltage profile, highlighting
the stealthiness of the attack strategy. In Case 1, the attack vector
${Av}_{C}$ increases to 0.08 at nodes N102, N103, and N104, doubling the
normal load, while other phases reduce to a mere 0.001. This results in 26
voltage violations with an unbalance of 1.704%, indicating a significant yet
covert disruption, as the unbalance stays below typical detection thresholds.
Case 2 mirrors this strategy with attack vector ${Av}_{B}$, again leading to
26 violations but with a slightly higher unbalance of 1.888%. This subtle
increase in unbalance, while still remaining undetectable by standard systems,
hints at the attack’s evolving sophistication. Case 3 focuses on Phase C with
similar intensity, yet it results in only three violations, albeit with a
higher voltage unbalance of 2.000%. The reduced number of violations contrasts
with the increased voltage unbalance percentage, showcasing the attack’s
ability to impact network stability without triggering immediate alarms.

In Cases 4, 5, and 6, the attack intensity escalates further to 0.16,
significantly overloading one phase while reducing others to near switch-off
levels. These cases demonstrate the attack’s increasing severity: Case 4
results in 26 violations with a 2.516% unbalance, Case 5 has the same number
of violations, but a slightly higher unbalance of 2.574%, and Case 6 marks the
peak of the attack’s impact with 36 violations and a 2.999% unbalance. These
latter cases show a clear departure from the stealthiness observed in earlier
scenarios, as the magnitude of the voltage violations and unbalance soar
beyond standard limits, posing a severe and overt threat to the network’s
stability. Fig. 8 vividly illustrates the magnitude of voltage violations
resulting from the attack vector.

In the IEEE-123 bus system, switches S1 to S6 are normally closed (on),
ensuring a continuous power flow. Conversely, switches S7 and S8 are typically
open (off), allowing for flexible network management and control within the
system. The efficacy of post-attack topology control measures, implemented
through a game-theoretic approach as detailed in Table III, is clearly
demonstrated in their ability to effectively mitigate voltage violations. In
Cases 1 through 5, the strategic operation of switches, specifically S7 and
S8, successfully eliminates all voltage violations, highlighting the
approach’s precision in maintaining network stability. Notably, the game-
theoretic method excels in optimizing switch states with minimal changes - a
crucial factor for real-time operation feasibility. In most cases, this
results in maintaining S8 in the ’Off’ state while selectively operating S7,
thus minimizing the switching actions and preserving the radiality of the
system.

Case 6, however, presents a more complex scenario; despite applying similar
switch adjustments, with both S7 and S8 ’On’, the network experiences 25
instead of 36 violations. In addition, the game theory approach explores
various switch combinations for this scenario, but it ultimately determines
that the current configuration is the most optimal for voltage regulation.
This particular outcome underscores the challenge posed by varying attack
intensities and the necessity for a more nuanced control strategy. The game-
theoretic approach’s adaptability and dynamic response are essential in these
scenarios, allowing for real-time adjustments to the network’s topology. This
adaptability not only ensures the mitigation of voltage violations but also
maintains the essential radial structure of the power system. Furthermore,
Fig. 8 shows a significant reduction in voltage violations, with the impacted
nodes now comfortably within voltage standard limits.

In our study, adjusting switches S7 and S8 effectively mitigates voltage
violations from cyber attacks on smart meters. While this highlights their
importance in network stability, it is important to recognize that this
strategy’s effectiveness may vary in different scenarios.

Figure 6: IEEE-123 draft in RSCAD Figure 7: IEEE-123 Runtime in RSCAD TABLE
III: Impact of Topology Changes on Voltage Violations

Case | S1 | S2 | S3 | S4 | S5 | S6 | S7 | S8 | Violations Count
---|---|---|---|---|---|---|---|---|---
1 | On | On | On | On | On | On | On | Off | 0
2 | On | On | On | On | On | On | On | Off | 0
3 | On | On | On | On | On | On | On | Off | 0
4 | On | On | On | On | On | On | On | On | 0
5 | On | On | On | On | On | On | On | Off | 0
6 | On | On | On | On | On | On | On | On | 25

TABLE IV: Impact of Attack Vectors on Voltage Violations and Unbalance Case | Attack Vector ${Av}_{C}$ | Attack Vector ${Av}_{B}$ | Attack Vector ${Av}_{A}$ |  Violations | Voltage Unbalance (%)
---|---|---|---|---|---
|  N102 |  N103 |  N104 |  N106 |  N107 |  N99 |  N109 |  N111 |  N114 |  |
1 |  0.08 |  0.08 |  0.08 |  0.001 |  0.001 |  0.001 |  0.001 |  0.001 |  0.001 |  26 | 1.704
2 |  0.001 |  0.001 |  0.001 |  0.08 |  0.08 |  0.08 |  0.001 |  0.001 |  0.001 |  26 | 1.888
3 |  0.001 |  0.001 |  0.001 |  0.001 |  0.001 |  0.001 |  0.08 |  0.08 |  0.08 |  3 | 2.000
4 |  0.16 |  0.16 |  0.16 |  0.001 |  0.001 |  0.001 |  0.001 |  0.001 |  0.001 |  26 | 2.516
5 |  0.001 |  0.001 |  0.001 |  0.16 |  0.16 |  0.16 |  0.001 |  0.001 |  0.001 |  26 | 2.574
6 |  0.001 |  0.001 |  0.001 |  0.001 |  0.001 |  0.001 |  0.16 |  0.16 |  0.16 |  36 | 2.999

(a) Case 1

(b) Case 2

(c) Case 3

(d) Case 4

(e) Case 5

(f) Case 6

Figure 8: Mitigation of voltage violations via topology control in response to
load altering attacks

##  VI Conclusion and Future Work

This study successfully demonstrated the integration of Python with RTDS and
RTAC, creating a robust cyber-physical testbed for real-time power system
simulations. We designed and executed a stealthy cyber attack on smart meters,
emphasizing the importance of advanced detection and mitigation strategies in
modern power systems. Our game-theoretic approach for topology control proved
effective in maintaining network stability and minimizing voltage violations,
even under challenging attack scenarios.

Looking ahead, our future work will concentrate on expanding the scope and
capabilities of our testbed. This includes enhancing data handling to manage
more extensive networks and addressing the limitations of Modbus for more
comprehensive data acquisition by using other industrial protocols like DNP3
or Ethernet/IP. We also aim to explore the integration of distributed energy
resources (DERs) control alongside topology control in our cyber attack
mitigation strategies. By incorporating machine learning algorithms, we aspire
to develop more dynamic and predictive control strategies that adapt
proactively to the evolving landscape of cyber threats.

##  VII Acknowledgment

The authors are thankful to Schweitzer Engineering Lab (SEL) and RTDS for
their support in successfully implementing this project.

## References

  * [1] S. Raschka, J. Patterson, and C. Nolet, “Machine learning in python: Main developments and technology trends in data science, machine learning, and artificial intelligence,” Information, vol. 11, no. 4, p. 193, 2020.
  * [2] D. Kornienko, S. Mishina, S. Shcherbatykh, and M. Melnikov, “Principles of securing restful api web services developed with python frameworks,” in Journal of Physics: Conference Series, vol. 2094, p. 032016, IOP Publishing, 2021.
  * [3] M. Xue and C. Zhu, “The socket programming and software design for communication based on client/server,” in 2009 Pacific-Asia Conference on Circuits, Communications and Systems, pp. 775–777, IEEE, 2009.
  * [4] N. Goldenberg and A. Wool, “Accurate modeling of modbus/tcp for intrusion detection in scada systems,” international journal of critical infrastructure protection, vol. 6, no. 2, pp. 63–75, 2013.
  * [5] M. M. S. Khan, A. Palomino, J. Brugman, J. Giraldo, S. K. Kasera, and M. Parvania, “The cyberphysical power system resilience testbed: Architecture and applications,” Computer, vol. 53, no. 5, pp. 44–54, 2020.
  * [6] H. Huang, C. M. Davis, and K. R. Davis, “Real-time power system simulation with hardware devices through dnp3 in cyber-physical testbed,” in 2021 IEEE Texas Power and Energy Conference (TPEC), pp. 1–6, IEEE, 2021.
  * [7] B. Xiao, M. Starke, G. Liu, B. Ollis, P. Irminger, A. Dimitrovski, K. Prabakar, K. Dowling, and Y. Xu, “Development of hardware-in-the-loop microgrid testbed,” in 2015 IEEE energy conversion congress and exposition (ECCE), pp. 1196–1202, IEEE, 2015.
  * [8] I. A. Oyewumi, A. A. Jillepalli, P. Richardson, M. Ashrafuzzaman, B. K. Johnson, Y. Chakhchoukh, M. A. Haney, F. T. Sheldon, and D. C. de Leon, “Isaac: The idaho cps smart grid cybersecurity testbed,” in 2019 IEEE Texas Power and Energy Conference (TPEC), pp. 1–6, IEEE, 2019.
  * [9] F. Shariatzadeh, C. B. Vellaithurai, S. S. Biswas, R. Zamora, and A. K. Srivastava, “Real-time implementation of intelligent reconfiguration algorithm for microgrid,” IEEE Trans. sustainable energy, vol. 5, no. 2, pp. 598–607, 2014.
  * [10] S. S. Biswas, F. Shariatzadeh, R. Beckstrom, and A. K. Srivastava, “Real time testing and validation of smart grid devices and algorithms,” in 2013 IEEE Power & Energy Society General Meeting, pp. 1–5, IEEE, 2013.
  * [11] J. Kim and L. Tong, “On topology attack of a smart grid: Undetectable attacks and countermeasures,” IEEE Journal on Selected Areas in Communications, vol. 31, no. 7, pp. 1294–1305, 2013.
  * [12] A. Selim, J. Zhao, F. Ding, F. Miao, and S.-Y. Park, “Deep reinforcement learning for distribution system cyber attack defense with ders,” in 2023 IEEE Power & Energy Society Innovative Smart Grid Technologies Conference (ISGT), pp. 1–5, IEEE, 2023.
  * [13] Y. Liu, P. Ning, and M. K. Reiter, “False data injection attacks against state estimation in electric power grids,” ACM Transactions on Information and System Security (TISSEC), vol. 14, no. 1, pp. 1–33, 2011.
  * [14] A. Kusko and M. T. Thompson, Power quality in electrical systems, vol. 23.  McGraw-Hill New York, 2007.
  * [15] “Pymodbus documentation.” https://pymodbus.readthedocs.io/en/latest/, 2023.  Accessed: 2023-12-03.
  * [16] J. Brugman, M. Khan, S. Kasera, and M. Parvania, “Cloud based intrusion detection and prevention system for industrial control systems using software defined networking,” in 2019 Resilience Week (RWS), vol. 1, pp. 98–104, IEEE, 2019.


# Distance between various discretized fermion actions

Dian-Jun Zhao1,2, Gen Wang3,4, Fangcheng He1, Luchang Jin5,6, Peng Sun7, Yi-Bo
Yang1,2,8,9, Kuan Zhang1,2 ($\chi$QCD Collaboration) 1University of Chinese
Academy of Sciences, School of Physical Sciences, Beijing 100049, China
2CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics,
Chinese Academy of Sciences, Beijing 100190, China
3Department of Physics and Astronomy, University of Kentucky, Lexington, KY
40506, USA
4Aix-Marseille Université, Université de Toulon, CNRS, CPT, Marseille, France
5Physics Department, University of Connecticut, Storrs, CT 06269, USA
6RIKEN-BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973,
USA
7Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou, 730000,
China
8School of Fundamental Physics and Mathematical Sciences, Hangzhou Institute
for Advanced Study, UCAS, Hangzhou 310024, China
9International Centre for Theoretical Physics Asia-Pacific, Beijing/Hangzhou,
China

###### Abstract

We present the leading order mixed-action effect $\Delta_{\rm mix}\equiv
m_{\pi,{\rm vs}}^{2}-\frac{m_{\pi,{\rm vv}}^{2}+m_{\pi,{\rm ss}}^{2}}{2}$
using HISQ, clover or overlap valence fermion actions on the gauge ensembles
with kinds of sea fermion actions among a widely used lattice spacing range
$a\in[0.04,0.19]$ fm. The results suggest that $\Delta_{\rm mix}$ decreases on
the forth order of the lattice spacing on the gauge ensembles with the
dynamical chiral sea fermion, likes the Domain wall or HISQ fermion. When the
clover sea fermion action which has explicit chiral symmetry breaking is used
in the ensemble, $\Delta_{\rm mix}$ can be much larger regardless of the
valence fermion action used.

Introduction: Lattice provide an unique gauge invariant, non-perturbative
regularization method for non-ablian gauge field theories. However, the
infamous fermion doubling problem prevents a straight forward discretiztion of
the continuum Dirac fermion action being used in 4-dimensional lattice QCD
calculations. The overlap (OV) fermion Narayanan and Neuberger (1995);
Neuberger (1998); Chiu and Zenkin (1999) which satisfies the Ginsburg-Wilson
relation, would be the optimal choice for the discretized Dirac operator, but
requires a factor of ${\cal O}(100)$ cost of computational resources comparing
to the widely used Wilson-like fermions; on the other hand, the staggered
fermion and its improved versions can also provide exact chiral symmetry with
a cost much lower than the Wilson-like fermion, at the expense of mixing
between four equivalent “taste” of a given flavor.

Since the discretization error and loop corrections of different discretized
actions can be very different from each other (e.g, Ref Horsley  _et al._
(2004)), conventional lattice QCD study concentrates on given fermion action
to avoid additional complications. But using more expensive fermion action on
the ensembles with cheaper fermion action becomes popular in the past decade,
likes the works for the glue helicity Yang  _et al._ (2017), nucleon axial
charge Chang  _et al._ (2018); Jang  _et al._ (2020), hadron vacuum
polarization Borsanyi  _et al._ (2021) and so on. Thus one will be wondering
how good the above “mixed action” setup is at the finite lattice spacing $a$,
and how the “mixed action” effects vanishes in the continuum limit
$a\rightarrow 0$.

The mixed action partially quenched chiral perturbation theory (MAPQ$\chi$PT)
has been developed to quantify the mixed action effect Bar  _et al._ (2004,
2005). If we calculate the valence quark propagator using discretized fermion
action B on a gauge ensemble at pion mass $m_{\pi,{\rm ss}}$ with dynamical
light quark using action A, and tune the valence quark mass to make the
corresponding pion mass $m_{\pi,{\rm vv}}=m_{\pi,{\rm ss}}$, then the mass
$m_{\pi,{\rm vs}}$ of the pion correlator formed by one valence quark with
action B and one sea (anti-)quark with action A, can suffer from a
discretization error and be different from $m_{\pi,{\rm ss}}$. And then one
can define a low energy constant as

$\displaystyle\Delta^{\rm B/A}_{\rm mix}(m_{\pi,{\rm vv}},m_{\pi,{\rm
ss}},a)\equiv m_{\pi,{\rm vs}}^{2}-\frac{m_{\pi,{\rm vv}}^{2}+m_{\pi,{\rm
ss}}^{2}}{2},$ (1)

which is independent of $m_{\pi,{\rm vv}}$ or $m_{\pi,{\rm ss}}$ at leading
order of $\chi$PT.

Assuming $m_{\pi,{\rm vv}}=m_{\pi,{\rm ss}}\sim 300$ MeV and $a\sim$0.09 fm,
the mixed action effect will make $\delta m_{\pi}\equiv m_{\pi,{\rm
vs}}-m_{\pi,{\rm ss}}$ to be 153(59) MeV for the overlap valence quark on the
clover (CL, a typical Wilson-like fermion) sea Durr  _et al._ (2007), while
$\delta m_{\pi}$ will reduce to 30-60 MeV for the combination of the Domain
wall fermion (DW, a different implementation of the overlap fermion by
introducing an extra fifth dimension for the fermion) valence and staggered
sea Orginos and Walker-Loud (2008); Aubin  _et al._ (2008); Berkowitz  _et
al._ (2017), and it can be as small as $\sim 10$ MeV if we use the overlap
valence quark on the DW sea Lujan  _et al._ (2012). There are also studies on
the other combinations likes overlap on highly improved staggered quark (HISQ)
Follana  _et al._ (2007) at $a\simeq 0.12$ fm Basak  _et al._ (2015).

But the all above estimates are based on the imprecise calculation at
$a\sim$0.1 fm, and then can misunderstand the situation. In this work, we
calculate the mixed action effect with kinds of the valence and sea fermion
combinations, and the results suggest that $\Delta_{\rm mix}\propto a^{4}$
when the fermion action used in the sea has the chiral symmetry.

Methodology and setup: The naive fermion action below has the infamous fermion
doubling problem, and it leads to 16 fermions as opposed to 1 in 4-dimension.
The Staggered fermion introduces a redefinition of the quark field (e.g.,
$\psi_{\rm
st}(x)=\gamma_{4}^{x_{4}}\gamma_{1}^{x_{1}}\gamma_{2}^{x_{2}}\gamma_{3}^{x_{3}}\psi(x)$
at the site $x=\\{x_{1},x_{2},x_{3},x_{4}\\}$ in the MILC conversion) to
partially fix the fermion doubling problem, and further improvement likes the
HISQ action Follana  _et al._ (2007) is needed to suppress the mixing between
4 residual fermions. On the other hand, the Wilson fermion can avoid the
entire doubling problem, but it introduces an explicit chiral symmetry
breaking effect on the quark mass. Such an effect can be improved by adding a
clover term. The solution to avoid the entire fermion doubling problem without
explicit chiral symmetry breaking is the Ginsparg-Wilson action Ginsparg and
Wilson (1982), likes the overlap fermion Chiu (1999); Liu (2005). More details
of the fermion actions can be found in the Supplemental materials ref .

At the next-to-leading-order (NLO) of the partially quenched $\chi$PT, the
pion masses used in the mixed action effect definition Eq. (1) can be non-zero
when $m_{\pi,{\rm vv}}\neq m_{\pi,{\rm ss}}$, even the fermion actions used by
the valence and sea quarks are the same. Thus in this work, we majorly
concentrate on the unitary case with the valence pion mass tuned to be the
same as the sea pion mass,

$\displaystyle\Delta_{\rm mix,uni}(m_{\pi},a)\equiv\Delta_{\rm
mix}(m_{\pi},m_{\pi},a).$ (2)

$\Delta_{\rm mix,uni}$ approaches exactly zero in the continuum and then is
good reference to verify the additional discretization error in the mixed
fermion action simulation. Such a definition is different from that used in
the previous $\chi$QCD study Lujan  _et al._ (2012) but still provides
consistent result, and detailed comparison can be found in the supplemental
materials ref .

Since the HYP smearing Hasenfratz and Knechtli (2001) can make the cost of the
overlap fermion calculation to be much cheaper and the fluctuation of the
clover fermion around the physical pion mass to be smaller, we will
concentrate on the following three kinds of the valence fermion actions to
calculate $\Delta_{\rm mix}$:

1) HI: HISQ action without any additional smearing on the gauge link;

2) CL: Clover fermion with 1-step HYP smearing and tree level tadpole improved
clover coefficient $c_{sw}=\langle U_{p}\rangle^{-3/4}$ where $\langle
U_{p}\rangle$ is the vacuum expectation value of the plaquatte on the HYP
smeared configurations;

3) OV: Overlap fermion with 1-step HYP smearing and $\rho=1.5$.

Table 1: Information of the ensembles Fukaya  _et al._ (2011); Bazavov  _et al._ (2010, 2013, 2018); Blum  _et al._ (2016); Boyle  _et al._ (2016); Zhang  _et al._ (2021) used in this calculation. Action | Symbol | $L^{3}\times T$ |  $a$ (fm) |  $m_{\pi}$ (MeV) |
---|---|---|---|---|---
DW+ID | 24D | $24^{3}\times\ 64$ | 0.194 | 139 |
DW+ID | 32Df | $32^{3}\times\ 64$ | 0.143 | 143 |
DW+I | 48I | $48^{3}\times\ 96$ | 0.114 | 139 |
DW+I | 24I | $24^{3}\times\ 64$ | 0.111 | 340 |
DW+I | 64I | $64^{3}\times 128$ | 0.084 | 139 |
DW+I | 32I | $32^{3}\times\ 64$ | 0.083 | 302 |
OV0HYP+IR | JLQCD | $24^{3}\times\ 48$ | 0.112 | 290 |
HI+S(1) | a12m310 | $24^{3}\times\ 64$ | 0.121 | 310 |
HI+S(1) | a09m310 | $32^{3}\times\ 96$ | 0.088 | 310 |
HI+S(1) | a06m310 | $48^{3}\times 144$ | 0.057 | 310 |
HI+S(1) | a04m310 | $64^{3}\times 192$ | 0.043 | 310 |
CLstout+Stad | C11 | $24^{3}\times\ 72$ | 0.108 | 290 |
CLstout+Stad | C08 | $32^{3}\times\ 96$ | 0.080 | 300 |
CLstout+Stad | C06 | $48^{3}\times 144$ | 0.054 | 300 |

For the sea fermion actions, we use three kinds of the gauge ensembles to
cover the popular choices:

1) those from the RBC/UKQCD collaboration using the 2+1 flavor DW fermion
action and Iwasaki gauge action at four lattice spacings and physical pion
mass Blum  _et al._ (2016); Boyle  _et al._ (2016), while the ensembles at two
largest lattice spacing include the DSDR term Boyle  _et al._ (2016) (DW+ID,
DD for short) and those at finer two lattice spacings are not (DW+I, DW for
short);

2) those from the MILC collaboration using the 2+1+1 flavor HISQ fermion
action with one-loop Symanzik improved gauge action Hart  _et al._ (2009) at
four lattice spacings and $m_{\pi}=310$ MeV Bazavov  _et al._ (2010, 2013,
2018) (HI+S(1), HI for short);

3) those from the CLQCD collaboration using the 2+1 flavor tadpole improved
clover fermion action (with 1-step stout link smearing on the gauge link with
smearing parameter 0.125) and tadpole improved Symanzik gauge action at three
lattice spacing and $m_{\pi}\simeq 300$ MeV Zhang  _et al._ (2021)
(CLstout+Stad, CLstout for short).

We also repeat the calculation on a JLQCD ensemble Fukaya  _et al._ (2011) at
$a=0.112$ fm and $m_{\pi,{\rm ss}}=290\mathrm{MeV})$. This ensemble uses the
2+1 flavor overlap fermion ($\rho=1.3$ without any HYP smearing Hasenfratz and
Knechtli (2001) on the gauge link) and Iwasaki gauge action with additional
ratio of extra Wilson fermions and associated twisted mass bosonic spinors
Fukaya  _et al._ (2006). The brief information of all the ensembles used in
this work are collected in Table 1. More details of those ensembles can be
found in the Supplemental materials ref .

On the DW/DD ensembles, we use the existed point source DW fermion propagators
with the field sparsening compression Li  _et al._ (2021) and generate the
valence quark propagators with similar source positions; On the other
ensembles, we use the Coulomb gauge fixed wall source to take the advantage of
the $L^{3}$ enhancement of statistics. Note that the wall source using the
HISQ action will practically a grid source which picks only the even points in
each spacial direction, as the Dirac space has been mapped to the even/odd
sites and should be treated separately in the source. At the same time, we
also use the low mode substitution Li  _et al._ (2010) to suppress the
statistical uncertainty of the pion correlators using the overlap fermion on
the HISQ ensembles at small lattice spacings. The details of the unitary
Chroma Edwards and Joo (2005)+QUDA Clark  _et al._ (2010); Babich  _et al._
(2011); Clark  _et al._ (2016) interface for various fermion actions can be
found in Ref. Zhang  _et al._ (2022).

Figure 1: The mixed action effect $\Delta_{\rm mix}$ on the DW and DD
ensembles, as functions of the lattice spacing $a$. The symbol X/Y means the
case with valence fermion action X on the sea fermion action Y. The figure
also illustrated the $a^{4}$ (solid line) and $a^{2}$ (dash line) dependence
for comparison.

Results: The lattice spacing $a$ dependence of $\Delta_{\rm mix,uni}$ on the
DW/DD ensembles are shown in Fig. 1, and the symbol X/Y corresponds to the
case with valence fermion action X on the sea fermion action Y. If we take
$\Delta_{\rm mix}^{\rm
OV/DW}(a=0.114~{}\mathrm{fm})=0.00477(20)~{}\mathrm{GeV}^{2}$ and assume the
lattice spacing dependence is $a^{4}$ (black line), then the prediction at
0.084 fm will be 0.00141(6) and agree with the data 0.00141(9) perfectly. On
the other hand, the prediction with a simple $a^{2}$ dependence (dash linne)
will be 0.0026(1) and then more than 5$\sigma$ higher than the practical
result. We can also solve $n$ from the data using the $Ca^{n}$ form with
different valence fermion actions, and have $n=$ 4.0(2) (OV/DW, purple
crosses), 4.0(2) (HI/DW, red stars), and 4.0(7) (CL/DW, blue boxes),
respectively.

It is interesting to see that the similar calculations on two coarse lattice
spacings 0.141 and 0.194 fm provide similar power of the lattice spacing, as
$n=$ 4.8(2) (OV/DD, purple crosses), 4.2(7) (HI/DD, red stars), and 3.9(4)
(CL/DD, blue boxes), respectively; but they are a factor $\sim$2 smaller than
$\Delta^{X/DW}_{\rm mix}$ assuming an $a^{4}$ lattice spacing dependence, as
illustrated by the black line for the $\Delta_{\rm mix}^{\rm OV/DW}$ case.
Such a suppression of $\Delta_{\rm mix}$ would relate to the DSDR term, as the
setup of the DD and DW ensembles are the same except this term.

Figure 2: The mixed action effects of the overlap fermion (purple crosses) or
Clover fermion (blue boxes) on the HI+S(1) ensembles at four lattice spacings,
together with that of the two valence fermion actions (yellow stars). The
latter two types of the data points are shifted horizontally to improve the
visibility.

Since the simulation with physical pion mass is expensive, we turn to the
HI+S(1) ensembles Bazavov  _et al._ (2010, 2013) with $\sim 310$ MeV pion mass
to investigate the lattice spacing dependence of $\Delta_{\rm mix}$ with the
same simulation setup, and show the results in Fig. 2. More detailed
discussion on the quark mass dependence of $\Delta_{\rm mix}$ can be found in
the supplemental materials ref .

With the data at four lattice spacings, it is more obvious that the $a^{4}$
behaviors observed on the DW/DD ensembles are not accidental and can dominate
in a wide lattice spacing range. Note that the relative uncertainty at
$a=0.043$ fm is large as the mixed action effect just change the pion mass by
$\sim$1(1) MeV there, and then one can not exclude the possibility on the
dominance of the lower order lattice spacing dependence (e.g. ${\cal
O}(a^{2})$ as predicted by the lattice perturbative theory Bar  _et al._
(2004)).

In order to investigate the likeness between different fermion actions, we
define the “mixed action effect” of two valence fermion actions B and C on the
same gauge ensemble with sea fermion action A as:

$\displaystyle\bar{\Delta}^{\rm B+C/A}_{\rm mix,uni}(m_{\pi},a)\equiv
m_{\pi,{\rm BC}}^{2}$ $\displaystyle\quad\quad-\frac{m_{\pi,{\rm
BB}}^{2}+m_{\pi,{\rm CC}}^{2}}{2}|_{m_{\pi,{\rm BB}}=m_{\pi,{\rm
CC}}=m_{\pi,{\rm AA}}=m_{\pi}}.$ (3)

As shown in Fig. 2, $\bar{\Delta}^{\rm CL+OV/HI}$ also decays on $a^{4}$ while
that at 0.06 fm has large uncertainty, and satisfy the following inequalities
within the statistical uncertainty,

$\displaystyle|\Delta^{\rm B/A}_{\rm mix,uni}-\Delta^{\rm C/A}_{\rm mix,uni}|$
$\displaystyle\leq\bar{\Delta}^{\rm B+C/A}_{\rm mix,uni}\leq\Delta^{\rm
B/A}_{\rm mix,uni}+\Delta^{\rm C/A}_{\rm mix,uni}.$ (4)

Further more, based on the similar calculation on the RBC/UKQCD ensembles for
three kinds of the valence fermion actions, the following relations are also
satisfied for $\bar{\Delta}^{\rm B+C/A}_{\rm mix,uni}$,

$\displaystyle|\bar{\Delta}^{\rm B+D/A}_{\rm mix,uni}-\bar{\Delta}^{\rm
C+D/A}_{\rm mix,uni}|$ $\displaystyle\leq\bar{\Delta}^{\rm B+C/A}_{\rm
mix,uni}\leq\bar{\Delta}^{\rm B+D/A}_{\rm mix,uni}+\bar{\Delta}^{\rm
C+D/A}_{\rm mix,uni},$ (5)

where D is a valence fermion action different from B and C.

Figure 3: The mixed action effects $\Delta_{\rm mix,uni}$ for kinds of the
valence fermion actions (distances between the black star and colored data
points, dash lines) and $\Delta_{\rm mix,uni}$ for different valence fermion
actions (distances between colored data points, solid lines) on the OV0HYP+IR
(left top), DW+I (right top), HI+S(1) (left bottom), and CLstout+Stad (right
bottom) ensembles. All the distances are in unit of GeV2, and interpolated to
$a=0.11$ fm and $m_{\pi}\in[290,310]$ MeV to make a fair comparison.

After interpolated to $a\simeq 0.11$ fm and $m_{\pi}\in[290,310]$ MeV, we show
the values of $\Delta_{\rm mix,uni}$ and $\bar{\Delta}_{\rm mix,uni}$ in Fig.
3, on the OV0HYP+IR, DW+I, HI+S(1) and CLstout+Stad ensembles. The black star
in each panel is fixed to the original point, and the distances between it and
red box, blue triangle and purple crosses correspond to $\Delta_{\rm mix,uni}$
of the HISQ, clover and overlap valence fermion actions, respectively. The
distance between the colored data points correspond to $\bar{\Delta}_{\rm
mix,uni}$ of different valence fermion actions. The uncertainties of
$\Delta_{\rm mix,uni}$ and $\bar{\Delta}_{\rm mix,uni}$ are smaller than 10%
in most of the cases and then ignored in the figure.

One can see that the relations Eq. (4-5) hold in all the cases, while the
“distances” between the fermion actions can be very different. The mixed
action effect on the OV0HYP+IR ensemble (left top panel) seems to be smaller
than the DW+I case (right top panel), while can be comparable with the DW+ID
case given the factor of 2 suppression shown in Fig. 1. It is predictable as
the DSDR term provides somehow similar impact on the near zero modes of the
Dirac operator as the ratio term Fukaya  _et al._ (2006); Vranas (2006).

Comparing to the mixed action effects on the DW+I ensembles, those on the
HI+S(1) ensemble (left bottom panel) are somehow similar. But when we consider
the valence CL or OV action, $\Delta_{\rm mix,uni}$ on the HI+S(1) ensemble
are still larger than those on the DW+I ensemble since the HI action is
somehow “far away” from either CL, OV or DW actions. At the same time, the
impact of the gauge action used by DW+I or HI+S(1) would be minor, and
requires further study using different fermion actions and the same gauge
action.

In contrary to the above three cases, the mixed action effects on the
CLstout+Stad ensemble (right bottom panel) are much larger, even for the case
of the HYP smeared clover valence fermion action on stout smeared clover sea
action . The lattice spacing dependence of $\Delta_{\rm mix,uni}$ on the
CLstout+Stad ensembles looks like $a^{2}$ at $a\sim 0.1$ fm but is closer to
$a^{4}$ when $a$ becomes smaller, and can be found in the supplemental
materials ref . The origin of the huge mixed action effect would relate to the
additive chiral symmetry breaking of the clover fermion action, and more
investigation is essential to figure out a good solution to suppress it.

Summary: Based on the calculation of the mixed action pion mass using three
kinds of the fermion actions, including overlap fermion, HISQ staggered
fermion, and clover Wilson fermion, on kinds of the dynamical ensembles, we
found that the leading mixed action effect would be $a^{4}$, and are
particularly small when the sea fermion action has the chiral symmetry, much
smaller than the naive theoretical estimate at small lattice spacings. We note
that a similar $a^{4}$ behavior was observed in the taste mixing effect of the
HISQ fermion which is also found to be roughly $a^{4}$ (or
$\alpha_{s}^{n}a^{2}$ with large $n$) in Ref. Borsanyi  _et al._ (2021); Aubin
_et al._ (2022). At the same time, the ensembles using the clover fermion
action have much larger mixed action effect for all the three kinds of the
valence fermion actions, and then the mixed action calculation on those
ensembles can be more challenging. The simulations with the exactly the same
gauge action but different fermion action, should be very valuable to answer
whether changing the gauge action can be helpful to suppress the mixed action
effect with similar cost.

Among the ensembles we used in this work, the ensemble using the fermion
action satisfying the Ginsburg-Wilson relation and Iwasaki gauge action can
provide the smallest mixed action effect, and the ratio or DSDR term can
further suppress it with the price of the topological charge fluctuation
suppression Fukaya  _et al._ (2011); Boyle  _et al._ (2016); The HISQ fermion
and one-loop Symanzik improved gauge action ensembles have a larger mixed
action effects for either CL or OV but can still be an acceptable choice at
small lattice spacing thanks to the $a^{4}$ behavior. We noted that the
$a^{4}$ scaling of the mixed action effects observed in this work cannot be
explained with the partially quenched chiral perturbation theory
(MAPQ$\chi$PT), which predict a $a^{2}$ scaling. The mechanism of the $a^{4}$
scaling of $\Delta_{\rm mix}$ and whether the mixed action effects can display
a $a^{2}$ scaling in some other observables need further studies.

## Acknowledgment

We thank the CLQCD, JLQCD, RBC/UKQCD, and MILC collaborations for providing us
their gauge configurations with dynamical fermions. The calculations were
performed using the Chroma software suite Edwards and Joo (2005) with QUDA
Clark  _et al._ (2010); Babich  _et al._ (2011); Clark  _et al._ (2016) and
the GWU code Alexandru  _et al._ (2012, 2011) through HIP programming model Bi
_et al._ (2020). The domain wall propagators are generated using Grid, GPT,
and Qlattice. The numerical calculation has been done on CAS Xiaodao-1
computing environment, and HPC Cluster of ITP-CAS. This research used
resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge
National Laboratory, which is supported by the Office of Science of the U.S.
Department of Energy under Contract No. DE-AC05-00OR22725. This work used
Stampede time under the Extreme Science and Engineering Discovery Environment
(XSEDE), which is supported by National Science Foundation Grant No.
ACI-1053575. We also thank the National Energy Research Scientific Computing
Center (NERSC) for providing HPC resources that have contributed to the
research results reported within this paper. We acknowledge the facilities of
the USQCD Collaboration used for this research in part, which are funded by
the Office of Science of the U.S. Department of Energy. Y. Yang is supported
in part by the Strategic Priority Research Program of Chinese Academy of
Sciences, Grant No. XDB34030303 and XDPB15. P. Sun is supported in part by the
National Natural Science Foundation of China under Grant No. 11975127. Y. Yang
and P. Sun are also supported in part by a NSFC-DFG joint grant under grant
No. 12061131006 and SCHA 458/22. G. Wang is supported by the U.S. DOE Grant
No. DE-SC0013065, DOE Grant No. DE-AC05-06OR23177 which is within the
framework of the TMD Topical Collaboration and the French National Research
Agency under the contract ANR-20-CE31-0016. L.C.J. acknowledges support by DOE
Office of Science Early Career Award DE-SC0021147 and DOE grant DE-SC0010339.

## References

  * Narayanan and Neuberger (1995) R. Narayanan and H. Neuberger, Nucl. Phys. B443, 305 (1995), arXiv:hep-th/9411108 [hep-th] .
  * Neuberger (1998) H. Neuberger, Phys. Lett. B417, 141 (1998), arXiv:hep-lat/9707022 [hep-lat] .
  * Chiu and Zenkin (1999) T.-W. Chiu and S. V. Zenkin, Phys. Rev. D59, 074501 (1999), arXiv:hep-lat/9806019 [hep-lat] .
  * Horsley  _et al._ (2004) R. Horsley, H. Perlt, P. E. L. Rakow, G. Schierholz, and A. Schiller (QCDSF), Nucl. Phys. B693, 3 (2004), [Erratum: Nucl. Phys.B713,601(2005)], arXiv:hep-lat/0404007 [hep-lat] .
  * Yang  _et al._ (2017) Y.-B. Yang, R. S. Sufian, A. Alexandru, T. Draper, M. J. Glatzmaier, K.-F. Liu, and Y. Zhao, Phys. Rev. Lett. 118, 102001 (2017), arXiv:1609.05937 [hep-ph] .
  * Chang  _et al._ (2018) C. C. Chang, A. N. Nicholson, E. Rinaldi, E. Berkowitz, N. Garron, D. A. Brantley, H. Monge-Camacho, C. J. Monahan, C. Bouchard, M. A. Clark, and et al., Nature 558, 91 (2018).
  * Jang  _et al._ (2020) Y.-C. Jang, R. Gupta, B. Yoon, and T. Bhattacharya, Phys. Rev. Lett. 124, 072002 (2020), arXiv:1905.06470 [hep-lat] .
  * Borsanyi  _et al._ (2021) S. Borsanyi _et al._ , Nature 593, 51 (2021), arXiv:2002.12347 [hep-lat] .
  * Bar  _et al._ (2004) O. Bar, G. Rupak, and N. Shoresh, Phys. Rev. D 70, 034508 (2004), arXiv:hep-lat/0306021 .
  * Bar  _et al._ (2005) O. Bar, C. Bernard, G. Rupak, and N. Shoresh, Phys. Rev. D72, 054502 (2005), arXiv:hep-lat/0503009 [hep-lat] .
  * Durr  _et al._ (2007) S. Durr, Z. Fodor, C. hoelbling, S. D. katz, S. Krieg, T. Kurth, L. Lellouch, T. Lippert, K. K. Szabo, and G. Vulvert,  _Proceedings, Lattice 2007:_ , PoS LATTICE2007, 115 (2007), arXiv:0710.4769 [hep-lat] .
  * Orginos and Walker-Loud (2008) K. Orginos and A. Walker-Loud, Phys. Rev. D77, 094505 (2008), arXiv:0705.0572 [hep-lat] .
  * Aubin  _et al._ (2008) C. Aubin, J. Laiho, and R. S. Van de Water, Phys. Rev. D77, 114501 (2008), arXiv:0803.0129 [hep-lat] .
  * Berkowitz  _et al._ (2017) E. Berkowitz _et al._ , Phys. Rev. D96, 054513 (2017), arXiv:1701.07559 [hep-lat] .
  * Lujan  _et al._ (2012) M. Lujan, A. Alexandru, Y. Chen, T. Draper, W. Freeman, M. Gong, F. X. Lee, A. Li, K. F. Liu, and N. Mathur, Phys. Rev. D86, 014501 (2012), arXiv:1204.6256 [hep-lat] .
  * Follana  _et al._ (2007) E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong (HPQCD, UKQCD), Phys. Rev. D 75, 054502 (2007), arXiv:hep-lat/0610092 .
  * Basak  _et al._ (2015) S. Basak, S. Datta, N. Mathur, A. T. Lytle, P. Majumdar, and M. Padmanath (ILGTI), PoS LATTICE2014, 083 (2015), arXiv:1412.7248 [hep-lat] .
  * Ginsparg and Wilson (1982) P. H. Ginsparg and K. G. Wilson, Phys. Rev. D25, 2649 (1982).
  * Chiu (1999) T.-W. Chiu, Phys. Rev. D60, 034503 (1999), arXiv:hep-lat/9810052 [hep-lat] .
  * Liu (2005) K.-F. Liu, Int. J. Mod. Phys. A20, 7241 (2005), arXiv:hep-lat/0206002 [hep-lat] .
  * (21) Supplementary materials .
  * Hasenfratz and Knechtli (2001) A. Hasenfratz and F. Knechtli, Phys. Rev. D64, 034504 (2001), arXiv:hep-lat/0103029 [hep-lat] .
  * Fukaya  _et al._ (2011) H. Fukaya, S. Aoki, T. Chiu, S. Hashimoto, T. Kaneko, J. Noaki, T. Onogi, and N. Yamada (JLQCD, TWQCD), Phys. Rev. D 83, 074501 (2011), arXiv:1012.4052 [hep-lat] .
  * Bazavov  _et al._ (2010) A. Bazavov _et al._ (MILC), Phys. Rev. D 82, 074501 (2010), arXiv:1004.0342 [hep-lat] .
  * Bazavov  _et al._ (2013) A. Bazavov _et al._ (MILC), Phys. Rev. D 87, 054505 (2013), arXiv:1212.4768 [hep-lat] .
  * Bazavov  _et al._ (2018) A. Bazavov _et al._ , Phys. Rev. D 98, 074512 (2018), arXiv:1712.09262 [hep-lat] .
  * Blum  _et al._ (2016) T. Blum _et al._ (RBC, UKQCD), Phys. Rev. D93, 074505 (2016), arXiv:1411.7017 [hep-lat] .
  * Boyle  _et al._ (2016) P. Boyle _et al._ , Phys. Rev. D 93, 054502 (2016), arXiv:1511.01950 [hep-lat] .
  * Zhang  _et al._ (2021) Q.-A. Zhang, J. Hua, F. Huang, R. Li, Y. Li, C.-D. Lu, P. Sun, W. Sun, W. Wang, and Y.-B. Yang, (2021), arXiv:2103.07064 [hep-lat] .
  * Hart  _et al._ (2009) A. Hart, G. M. von Hippel, and R. R. Horgan (HPQCD), Phys. Rev. D 79, 074008 (2009), arXiv:0812.0503 [hep-lat] .
  * Fukaya  _et al._ (2006) H. Fukaya, S. Hashimoto, K.-I. Ishikawa, T. Kaneko, H. Matsufuru, T. Onogi, and N. Yamada (JLQCD), Phys. Rev. D74, 094505 (2006), arXiv:hep-lat/0607020 [hep-lat] .
  * Li  _et al._ (2021) Y. Li, S.-C. Xia, X. Feng, L.-C. Jin, and C. Liu, Phys. Rev. D 103, 014514 (2021), arXiv:2009.01029 [hep-lat] .
  * Li  _et al._ (2010) A. Li _et al._ ($\chi$QCD), Phys. Rev. D82, 114501 (2010), arXiv:1005.5424 [hep-lat] .
  * Edwards and Joo (2005) R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD), Nucl. Phys. B Proc. Suppl. 140, 832 (2005), arXiv:hep-lat/0409003 .
  * Clark  _et al._ (2010) M. A. Clark, R. Babich, K. Barros, R. C. Brower, and C. Rebbi, Comput. Phys. Commun. 181, 1517 (2010), arXiv:0911.3191 [hep-lat] .
  * Babich  _et al._ (2011) R. Babich, M. A. Clark, B. Joo, G. Shi, R. C. Brower, and S. Gottlieb, in  _Proceeding, SC11_ (2011) arXiv:1109.2935 [hep-lat] .
  * Clark  _et al._ (2016) M. A. Clark, B. Joó, A. Strelchenko, M. Cheng, A. Gambhir, and R. Brower, (2016), arXiv:1612.07873 [hep-lat] .
  * Zhang  _et al._ (2022) K. Zhang, W. Sun, Y.-B. Yang, and R.-Q. Zhang, in  _38th International Symposium on Lattice Field Theory_ (2022) arXiv:2201.09004 [hep-lat] .
  * Vranas (2006) P. M. Vranas, Phys. Rev. D 74, 034512 (2006), arXiv:hep-lat/0606014 .
  * Aubin  _et al._ (2022) C. Aubin, T. Blum, M. Golterman, and S. Peris, (2022), arXiv:2204.12256 [hep-lat] .
  * Alexandru  _et al._ (2012) A. Alexandru, C. Pelissier, B. Gamari, and F. Lee, J. Comput. Phys. 231, 1866 (2012), arXiv:1103.5103 [hep-lat] .
  * Alexandru  _et al._ (2011) A. Alexandru, M. Lujan, C. Pelissier, B. Gamari, and F. X. Lee, in  _Proceedings, SAAHPC’11_ (2011) pp. 123–130, arXiv:1106.4964 [hep-lat] .
  * Bi  _et al._ (2020) Y.-J. Bi, Y. Xiao, W.-Y. Guo, M. Gong, P. Sun, S. Xu, and Y.-B. Yang,  _Proceedings, Lattice 2019_ , PoS LATTICE2019, 286 (2020), arXiv:2001.05706 [hep-lat] .
  * Luscher and Weisz (1985) M. Luscher and P. Weisz, Commun. Math. Phys. 97, 59 (1985), [Erratum: Commun.Math.Phys. 98, 433 (1985)].
  * Orginos and Toussaint (1999) K. Orginos and D. Toussaint (MILC), Phys. Rev. D 59, 014501 (1999), arXiv:hep-lat/9805009 .
  * Sharpe (1997) S. R. Sharpe, Phys. Rev. D 56, 7052 (1997), [Erratum: Phys.Rev.D 62, 099901 (2000)], arXiv:hep-lat/9707018 .

## Supplemental materials

###  .1 Note on the discretized fermion and gauge actions

The naive fermion action below has the infamous fermion doubling problem. In
4-dimension, it leads to 16 fermions as opposed to 1,

$\displaystyle S_{\rm naive}(m)=\sum_{x,y}\bar{\psi}(x)D_{\rm
naive}(x,y,m)\psi(y),$ (6)

with $D_{\rm
naive}(x,y,m)=\frac{1}{2}\sum_{\mu=1,...,4,\eta=\pm}\big{[}\eta\gamma_{\mu}U_{\mu}(x,x+\eta\hat{n}_{\mu}a)\delta_{y,x+\eta\hat{n}_{\mu}a}\big{]}-ma\delta_{y,x}$.
The Staggered fermion

$\displaystyle S_{\rm st}(m)=\sum_{x,y}\bar{\psi}_{\rm st}(x)D_{\rm
st}(x,y,m)\psi_{\rm st}(y),$ (7)

introduces a redefinition of the quark field (e.g., $\psi_{\rm
st}(x)=\gamma_{4}^{x_{4}}\gamma_{1}^{x_{1}}\gamma_{2}^{x_{2}}\gamma_{3}^{x_{3}}\psi(x)$
at the site $x=\\{x_{1},x_{2},x_{3},x_{4}\\}$ in the MILC conversion) and
rotate $D_{\rm naive}$ into $D_{\rm
st}(x,y,m)=\frac{1}{2}\sum_{\mu=1,...,4,\eta=\pm}\big{[}\eta\eta_{\mu}U_{\mu}(x,x+\eta\hat{n}_{\mu}a)\delta_{y,x+\eta\hat{n}_{\mu}a}\big{]}-ma\delta_{y,x}$,
where
$\\{\eta_{1},\eta_{2},\eta_{3},\eta_{4}\\}=\\{(-1)^{x_{4}},(-1)^{x_{1}+x_{4}},(-1)^{x_{1}+x_{2}+x_{4}},1\\}$
in MILC conversion. Staggered fermion does not completely fixed the fermion
doubling problem, there are 4 residual fermions, usually named as tastes of
fermion. Further improvement likes the HISQ action Follana  _et al._ (2007) is
needed to suppress the mixing between the tastes.

On the other hand, the Wilson fermion

$\displaystyle S_{\rm w}(m)=\sum_{x,y}\bar{\psi}(x)D_{\rm w}(x,y,m)\psi(y)$
(8)

with $D_{\rm w}(x,y,m)=D_{\rm
naive}(x,y,m)+\frac{1}{2}\sum_{\mu=1,...,4,\eta=\pm}U_{\mu}(x,x+\eta\hat{n}_{\mu}a)\delta_{y,x+\eta\hat{n}_{\mu}a}-4\delta_{y,x}$
can avoid the entire doubling problem. But it introduces an explicit chiral
symmetry breaking effect on the quark mass. This can be improved by adding a
clover term:

$\displaystyle S_{\rm cl}(m)=S_{\rm w}(m)+c_{\rm
sw}\sum_{x}\bar{\psi}(x)\sigma_{\mu\nu}F^{\mu\nu}\psi(x),$ (9)

where $c_{\rm sw}=1+{\cal O}(\alpha_{s})$ is a tunable parameter.

The solution to avoid the entire fermion doubling problem without explicit
chiral symmetry breaking is the Ginsparg-Wilson action Ginsparg and Wilson
(1982), likes the overlap fermion Chiu (1999); Liu (2005),

$\displaystyle S_{\rm c}(m)=\sum_{x}\bar{\psi}(x)[(1-\frac{m}{2\rho})D_{\rm
ov}(x,y,m)+\delta_{x,y}m]\psi(y),$ (10)

where $D_{\rm ov}=\rho(1+\frac{D_{\rm w}(-\rho)}{\sqrt{D^{\dagger}_{\rm
w}(-\rho)D_{\rm w}(-\rho)}})$ with a positive $\rho$ such as $\rho=1.5$.

The gauge actions in the ensembles we used in this work can be considered as
particular cases of the improved Luescher-Weisz action Luscher and Weisz
(1985),

$\displaystyle S=\sum_{x}\mathrm{Re}\big{[}c_{0}U_{\rm plaquattes}+c_{1}U_{\rm
Rectangles}+c_{2}U_{\rm parallelograms}\big{]},$ (11)

where $c_{0,1,2}$ are coefficients to be determined perturbatively or non-
perturbatively, $U_{\rm plaquattes}$ and $U_{\rm Rectangles}$ are the standard
$1\times 1$ and $2\times 2$ Wilson loop, and $U_{\rm parallelograms}$ is the
“parallelograms” Wilson loops with 6 links defined in Ref. Luscher and Weisz
(1985).

Table 2: Information of the ensembles Bazavov  _et al._ (2010, 2013, 2018); Blum  _et al._ (2016); Boyle  _et al._ (2016); Zhang  _et al._ (2021); Fukaya  _et al._ (2011) used in this section. The symbol I corresponds to the standard Iwasaki action, ID/IR means the Iwasaki action and extra DSDR or JLQCD ratio ratio term are used in the simulation, and S(1)/Stad is used for the Symanzik gauge action with fully 1-loop improvement or simple tree-level tadpole improvement, respectively. Action | Symbol | $L^{3}\times T$ | $c_{0}$ | $c_{1}$ | $c_{2}$ | $u_{0}$ |  $a$ (fm) |  $m_{\pi,{\rm ss}}$ (MeV)
---|---|---|---|---|---|---|---|---
DW+ID | 24DH | $24^{3}\times\ 64$ | 5.957 | -0.541 | 0 | 0.827 | 0.194 | 337
DW+ID | 24D | $24^{3}\times\ 64$ | 5.957 | -0.541 | 0 | 0.827 | 0.194 | 139
DW+ID | 32D | $32^{3}\times\ 64$ | 5.957 | -0.541 | 0 | 0.827 | 0.194 | 139
DW+ID | 32Df | $32^{3}\times\ 64$ | 6.384 | -0.579 | 0 | 0.846 | 0.143 | 143
DW+ID | 32DfH | $32^{3}\times\ 64$ | 6.384 | -0.579 | 0 | 0.846 | 0.143 | 250
DW+I | 48I | $48^{3}\times\ 96$ | 7.770 | -0.705 | 0 | 0.875 | 0.114 | 139
DW+I | 24I | $24^{3}\times\ 64$ | 7.770 | -0.705 | 0 | 0.876 | 0.111 | 340
DW+I | 24IH | $24^{3}\times\ 64$ | 7.770 | -0.705 | 0 | 0.876 | 0.111 | 432
DW+I | 24IH2 | $24^{3}\times\ 64$ | 7.770 | -0.705 | 0 | 0.876 | 0.111 | 576
DW+I | 64I | $64^{3}\times 128$ | 8.208 | -0.745 | 0 | 0.886 | 0.084 | 139
DW+I | 32I | $32^{3}\times\ 64$ | 8.208 | -0.745 | 0 | 0.886 | 0.083 | 302
DW+I | 32IH | $32^{3}\times\ 64$ | 8.208 | -0.745 | 0 | 0.886 | 0.083 | 412
OV0HYP+IR | JLQCD | $24^{3}\times\ 48$ | 8.390 | -0.761 | 0 | 0.884 | 0.112 | 290
HI+S(1) | a12m310 | $24^{3}\times\ 64$ | 6.000 | -0.298 | 0.022 | 0.864 | 0.121 | 310
HI+S(1) | a12m130 | $48^{3}\times\ 64$ | 6.000 | -0.298 | 0.022 | 0.864 | 0.121 | 130
HI+S(1) | a09m310 | $32^{3}\times\ 96$ | 6.300 | -0.281 | 0.021 | 0.874 | 0.088 | 310
HI+S(1) | a09m130 | $64^{3}\times\ 96$ | 6.300 | -0.281 | 0.021 | 0.874 | 0.088 | 130
HI+S(1) | a06m310 | $48^{3}\times 144$ | 6.720 | -0.263 | 0.020 | 0.886 | 0.057 | 310
HI+S(1) | a04m310 | $64^{3}\times 192$ | 7.000 | -0.254 | 0.019 | 0.892 | 0.043 | 310
CLstout+Stad | C11 | $24^{3}\times\ 72$ | 6.200 | -0.440 | 0 | 0.855 | 0.108 | 290
CLstout+Stad | C08 | $32^{3}\times\ 96$ | 6.410 | -0.430 | 0 | 0.863 | 0.080 | 300
CLstout+Stad | C06 | $48^{3}\times 144$ | 6.720 | -0.424 | 0 | 0.873 | 0.055 | 300

On the DW+I Blum  _et al._ (2016); Boyle  _et al._ (2016) and Fukaya  _et al._
(2011) ensembles, the Iwasaki gauge action uses the fixed coefficients
$c_{0,1,2}$ at arbitrary lattice spacing,

$\displaystyle c_{0}=(1-8c)\beta,\ c_{1}=c\beta,\ c_{2}=0,$ (12)

where $c=-0.331$, $\beta=6/g^{2}$ and $g$ is the bare coupling constant. On
the OV0HYP+IR ensemble, extra Wilson fermions and associated twisted mass
bosonic spinors are introduced to to generate a weight

$\displaystyle{\rm det}[H_{W}(-\rho)]/{\rm det}[H_{W}^{2}(-\rho)+m_{t}^{2}]$
(13)

with $m_{t}\sim 0.2$ in the functional integrals, to suppress the near-zero
mode of $H_{W}(-\rho)$ Fukaya  _et al._ (2006). Similarly, the DW+ID ensemble
uses an improved version of the above weight as the dislocation suppressing
determinant ratio (DSDR) Vranas (2006),

$\displaystyle{\rm det}[H_{W}^{2}(-\rho)+\epsilon_{f}^{2}]/{\rm
det}[H_{W}^{2}(-\rho)+\epsilon_{b}^{2}]$ (14)

with $\epsilon_{f}=0.02$ and $\epsilon_{b}=0.5$. As shown in Table 2, the
extra weight can make the $c_{0}$ of the OV0HYP+IR ensemble at 0.11 fm to be
$\sim 8$% larger than that of the DW+I ensemble at similar lattice spacing.

In the other hand, the one-loop improved Symanzik gauge action in the HI+S(1)
ensembles Bazavov  _et al._ (2010, 2013) uses the lattice spacing dependent
coefficients Hart  _et al._ (2009) with $N_{f}=4$,

$\displaystyle c_{0}=\beta^{\prime},\
c_{1}=-\frac{\beta^{\prime}}{20u_{0}^{2}}[1-(0.6264-1.1746N_{f})]\mathrm{ln}u_{0},\
c_{2}=\frac{\beta^{\prime}}{u_{0}^{2}}(0.0433-0.0156N_{f})\mathrm{ln}u_{0},$
(15)

where $u_{0}=\langle\frac{1}{18}U_{\rm plaquattes}\rangle^{1/4}$ is the
tadpole improvement factor which can be determined self-consistently during
the HMC, and the $\mathrm{ln}u_{0}$ dependence reflects the 1-loop correction
through the approximation $\alpha_{s}\simeq-\frac{4\mathrm{ln}u_{0}}{3.0684}$
Orginos and Toussaint (1999). Note that
$\beta^{\prime}=\frac{5}{3u_{0}^{4}}\beta$ also includes the additional
factors in front of the plaquatte term and then differs from $\beta$ by a
factor $\sim 3$.

At the same time, the gauge action used by the CLstout+Stad ensembles Zhang
_et al._ (2021) is the tree level Symanzik gauge action without the
$\mathrm{ln}u_{0}$ terms,

$\displaystyle c_{0}=\beta^{\prime},\
c_{1}=-\frac{\beta^{\prime}}{20u_{0}^{2}},\ c_{2}=0.$ (16)

Based on the information shown in Table 2, the similar Symanzik gauge actions
used by the HI+S(1) and CLstout+Stad ensembles with the same $c_{0}$ provides
quite close lattice spacing, even though they use different sea fermion
actions and different $c_{1,2}$. But the value of $u_{0}$ seems to be quite
sensitive to those differences.

###  .2 Comparison with the previous $\chi$QCD $\Delta_{\rm mix}$ calculation

At the next-to-leading-order (NLO) of the partially quenched $\chi$PT, the
pion masses used in the mixed action effect definition

$\displaystyle\Delta_{\rm mix}(m_{\pi,{\rm vv}},m_{\pi,{\rm ss}},a)$
$\displaystyle\equiv m_{\pi,{\rm vs}}^{2}-\frac{m_{\pi,{\rm
vv}}^{2}+m_{\pi,{\rm ss}}^{2}}{2},$ (17)

can be expressed as Sharpe (1997),

$\displaystyle m^{2}_{\pi,{\rm vv}}$ $\displaystyle=\Lambda_{\chi}^{2}2y_{\rm
v}\big{\\{}1+\frac{2}{N_{f}}[(2y_{\rm v}-y_{\rm s})\mathrm{ln}(2y_{\rm
v})+(y_{\rm v}-y_{\rm s})]+2y_{\rm v}(2\alpha_{8}-\alpha_{5})+2y_{\rm
s}N_{f}(2\alpha_{6}-\alpha_{4})+{\cal O}(y^{2})\big{\\}},$ (18) $\displaystyle
m^{2}_{\pi,{\rm vs}}$ $\displaystyle=\Lambda_{\chi}^{2}(y_{\rm v}+y_{\rm
s})\big{\\{}1+\frac{2}{N_{f}}y_{\rm v}\mathrm{ln}(2y_{\rm v})+(y_{\rm
v}+y_{\rm s})(2\alpha_{8}-\alpha_{5})+2y_{\rm
s}N_{f}(2\alpha_{6}-\alpha_{4})+{\cal O}(y^{2})\big{\\}},$ (19) $\displaystyle
m^{2}_{\pi,{\rm ss}}$ $\displaystyle=\Lambda_{\chi}^{2}2y_{\rm
s}\big{\\{}1+\frac{2}{N_{f}}y_{\rm s}\mathrm{ln}(2y_{\rm s})+2y_{\rm
s}[2\alpha_{8}-\alpha_{5}+N_{f}(2\alpha_{6}-\alpha_{4})]+{\cal
O}(y^{2})\big{\\}},$ (20)

where $\Lambda_{\chi}=4\pi F$ is the intrinsic scale of $\chi$PT with $F$
being the pion decay constant in the chiral limit, $y=\frac{\Sigma
m}{F^{2}\Lambda_{\chi}^{2}}$ is the dimensionless expansion parameter with
$\Sigma$ and $m$ being the chiral condensate and quark mass respectively, and
$\alpha_{i}=128\pi^{2}L_{i}$ is the NLO low energy constants of $\chi$PT. Thus
$\Delta_{\rm mix}(m_{\pi,{\rm vv}},m_{\pi,{\rm ss}},a)$ can be non-zero when
$m_{\pi,{\rm vv}}\neq m_{\pi,{\rm ss}}$, even the fermion actions used by the
valence and sea quarks are the same.

Figure 4: $\Delta_{\rm mix}(m_{\pi,{\rm
vv}},290~{}\mathrm{MeV},0.11~{}\mathrm{fm})$ on the OV0HYP+IR ensemble. The
$\Delta_{\rm mix}$ using the exactly the same action for both the valence and
sea quarks (black stars) is consistent with zero, while using the overlap
fermion action with HYP smeared link and $\rho=1.5$ (purple crosses) can make
$\Delta_{\rm mix}$ to be obviously non-zero while still small. The NLO
PQ$\chi$PT prediction of $\Delta_{\rm mix}$ is also shown as a black dash line
for comparison.

To quantify this effect, we carried out an unitary action calculation on the
OV0HYP+IR ensemble Fukaya  _et al._ (2011) from JLQCD collobaration at
$a=0.112$ fm and $m_{\pi,{\rm ss}}=290\mathrm{MeV})$. As in the Fig. 4, the
black stars show the partially quenched (PQ) results with the unitary OV0HYP
fermion action as that used in the ensemble, and the black dash curve shows
the NLO $\chi$PT prediction of $\Delta_{\rm mix}(m_{\pi,{\rm vv}},m_{\pi,{\rm
ss}}=290\mathrm{MeV})$ using the low energy constants from the literature.
They agree with each other well when $m_{\pi,{\rm vv}}$ is close to
$m_{\pi,{\rm ss}}$, and it is understandable that the lattice results can be
somehow different from the $\chi$PT prediction due to the finite volume effect
when $m_{\pi,{\rm vv}}$ is small (or the higher order $\chi$PT correction when
$m_{\pi,{\rm vv}}>400~{}\mathrm{MeV}$).

We also calculate $\Delta_{\rm mix}$ using the overlap fermion with HYP
smeared gauge link and $\rho=1.5$ which is the default $\chi$QCD setup, and
show the result (purple crosses) in Fig. 4. In such a case, the mixed action
effect at the unitary pion mass $\Delta_{\rm mix,uni}(m_{\pi,{\rm
ss}})\equiv\Delta_{\rm mix}(m_{\pi,{\rm ss}},m_{\pi,{\rm ss}})=0.0008(1)$ is
nonzero, and depends on the valence pion mass obviously.

In this work, we majorly concentrate on the unitary case with the valence pion
mass tuned to be the same as the sea pion mass,

$\displaystyle\Delta_{\rm mix,uni}(m_{\pi},a)\equiv\Delta_{\rm
mix}(m_{\pi},m_{\pi},a).$ (21)

$\Delta_{\rm mix,uni}$ approaches exactly zero in the continuum and then is
good reference to verify the additional discretization error in the mixed
fermion action simulation. Such a definition is different from that used the
previous $\chi$QCD study Lujan  _et al._ (2012), in which the mixed action
effect is defined by the following parameterization,

$\displaystyle m_{vs}^{2}-\frac{1}{2}m^{2}_{ss}=Bm_{v}+\tilde{\Delta}_{\rm
mix},$ (22)

and calculated on two DW+I ensembles with $\sim$300 MeV pion mass Blum  _et
al._ (2016). Based on above definition, Ref. Lujan  _et al._ (2012) gave
$\tilde{\Delta}^{\rm OV/DW}_{\rm mix}(0.083~{}\mathrm{fm})=0.004(2)$ and
$\tilde{\Delta}^{\rm OV/DW}_{\rm
mix}(0.111~{}\mathrm{fm})=0.009(2)~{}\mathrm{GeV}^{2}$ (based on Table V of
Ref. Lujan  _et al._ (2012) with Method II). If we assume a simple $Ca^{n}$
form and solve $n$, the data suggests an $a^{2.4(1.2)}$ behavior with large
uncertainty on the power.

Figure 5: The mixed action effect of the overlap fermion on the RBC/UKQCD
Domain wall fermion ensembles Blum  _et al._ (2016), as functions of the
valence pion mass $m_{\pi,\rm{vv}}$ with two lattice spacings and with/w.o.
HYP smearing. The lines show the fits in the range of
$m_{\pi,\rm{vv}}\in[0.05,0.15]~{}\mathrm{GeV}^{2}$.

We repeated the calculation on 24I and 32I ensembles with several valence
quark masses, and show the result as functions of valence pion mass in Fig. 5.
Based on the leading order $\chi$PT, $\tilde{\Delta}_{\rm mix}$ corresponds to
the $\Delta_{\rm mix}$ in the chiral limit of $m_{\pi,\rm{vv}}$. As shown in
the figures, the negative $m_{\pi,\rm{vv}}$ dependence of $\Delta_{\rm mix}$
(purple boxes for 24I and red cross for 32I) makes the $\tilde{\Delta}_{\rm
mix}$ (illustrated as the intercepts at $m_{\pi,\rm{vv}}=0$ at the leading
order of PQ$\chi$PT) to be larger than $\Delta_{\rm mix,uni}$, and consistent
with the previous study within the uncertainty. But with much more precise
data, similar parameterization suggests an $a^{4.6(0.4)}$ behavior which means
that the mixed action effect can be significantly smaller than what predicted
by a naive $a^{2}$ discretization error at small lattice spacings.

The results of $\Delta_{\rm mix}$ using the overlap valence fermion without
HYP smearing are also presented in Fig. 5 as the blue dots (24I) and green
triangles (32I). As shown in the figure, $\Delta_{\rm mix}$ without HYP
smearing is $\sim$ 20-30% than the HYP smeared cases at all the
$m_{\pi,\rm{vv}}$, and then have a similar $a^{4}$ lattice spacing dependence.

###  .3 Quark mass dependence and More about the clover fermion sea case

Figure 6: The mixed action effect $\Delta_{\rm mix}$ on the DW+I ensembles, as
functions of the lattice spacing $a$. The left panel shows the physical pion
mass case and the right panel shows that using 300 MeV pion mass. The figures
also illustrate the $a^{4}$ (solid line) and $a^{2}$ (dash line) dependence
for comparison.

Figure 7: The mixed action effects of the overlap fermion (purple crosses) or
Clover fermion (blue boxes) on the HI+S(1) ensembles using the HISQ fermion at
four lattice spacings, together with that of the two valence fermion actions
(yellow stars). The left panel shows the physical pion mass case and the right
panel shows that using 300 MeV pion mass. The latter two types of the data
points are shifted horizontally to improve the visibility.

In the main text, we show the lattice spacing dependence of $\Delta_{\rm mix}$
on the DW+I ensemble at physical pion mass, and also those on the HI+S(1) at
$m_{\pi}\simeq$ 300 MeV. In Fig. 6 and Fig. 7, we show the $\Delta_{\rm mix}$
on those two ensembles with both the physical pion mass (left panels) and
$m_{\pi}\simeq$ 300 MeV (right panels). Even though the value of $\Delta_{\rm
mix}$ can depend on the pion mass at given lattice spacing, but the lattice
spacing dependence are always $a^{4}$ in all the cases.

Figure 8: The mixed action effects of the overlap fermion (OV, purple
crosses), HISQ fermion (HI, red boxes) or HYP smeared clover fermion (CL, blue
triangles) on the Clover fermion ensembles at three lattice spacings,

For completeness, we also report the calculation on the CLstout+Stad ensembles
Zhang  _et al._ (2021) in the lattice spacing range $a\in[0.054,0.108]$ fm.
Unlike the DW+I or HI+S(1) ensembles using the chiral fermion, $\Delta_{\rm
mix,uni}$ at $a=0.080$ and 0.108 fm look like a $a^{2}$ behavior , but the
lattice spacing dependence becomes closer to $a^{4}$ from $a=0.080$ fm to
0.054 fm. At the same time, the mixed action effect of the overlap fermion is
larger than that of the HISQ fermion, which is also different from the
behavior on the DW+I or HI+S(1) ensembles which use the chiral sea fermion
actions. Since the CLstout+Stad ensembles and HI+S(1) ensembles using similar
Symanzik gauge action with minor difference on the 1-loop correction, the huge
$\Delta_{\rm mix,uni}$ would majorly related to the clover fermion action in
the ensembles, while more accurate comparison with the exactly the same gauge
action is worth further study in the future.