# Deflection of Coronal Mass Ejections in Unipolar Ambient Magnetic Fields

Michal Ben-Nun Predictive Science Inc., 9990 Mesa Rim Road, Suite 170, San
Diego, CA 92121, USA  Tibor Török Predictive Science Inc., 9990 Mesa Rim Road,
Suite 170, San Diego, CA 92121, USA  Erika Palmerio Predictive Science Inc.,
9990 Mesa Rim Road, Suite 170, San Diego, CA 92121, USA  Cooper Downs
Predictive Science Inc., 9990 Mesa Rim Road, Suite 170, San Diego, CA 92121,
USA  Viacheslav S. Titov Predictive Science Inc., 9990 Mesa Rim Road, Suite
170, San Diego, CA 92121, USA  Mark G. Linton U.S. Naval Research Laboratory,
4555 Overlook Avenue, SW Washington, DC 20375, USA  Ronald M. Caplan
Predictive Science Inc., 9990 Mesa Rim Road, Suite 170, San Diego, CA 92121,
USA  Roberto Lionello Predictive Science Inc., 9990 Mesa Rim Road, Suite 170,
San Diego, CA 92121, USA

###### Abstract

The trajectories of coronal mass ejections (CMEs) are often seen to
substantially deviate from a purely radial propagation direction. Such
deviations occur predominantly in the corona and have been attributed to
“channeling” or deflection of the eruptive flux by asymmetric ambient magnetic
fields. Here, we investigate an additional mechanism that does not require any
asymmetry of the pre-eruptive ambient field. Using magnetohydrodynamic
numerical simulations, we show that the trajectory of CMEs through the solar
corona can significantly deviate from a radial direction when propagation
takes place in a unipolar radial field. We demonstrate that the deviation is
most prominent below ${\sim}15\,R_{\odot}$ and can be attributed to an
“effective ${\bf I}\times{\bf B}$ force” that arises from the intrusion of a
magnetic flux rope with a net axial electric current into a unipolar
background field. These results are important for predictions of CME
trajectories in the context of space weather forecasts, as well as for
reaching a deeper understanding of the fundamental physics underlying CME
interactions with the ambient fields in the extended solar corona.

##  1 Introduction

Coronal mass ejections (CMEs) are huge expulsions of magnetized plasma from
the Sun’s atmosphere into interplanetary space, and the main driver of space-
weather effects at Earth (e.g., Baker & Lanzerotti, 2016; Temmer, 2021). CMEs,
and the often associated flares and prominence eruptions, are the largest
transient energy-release events in the solar system and constitute different
manifestations of a sudden and violent reconfiguration of the coronal magnetic
field (e.g., Forbes, 2000). While there is a general agreement that CMEs and
flares are powered by the free magnetic energy stored in current-carrying
structures in the low corona, there exists an ongoing debate on the physical
mechanisms that initiate and drive these events (e.g., Aulanier, 2014; Green
et al., 2018). As CMEs expand and then propagate through the corona and the
inner heliosphere, their interaction with the surrounding magnetic field and
solar wind, or with other CMEs, can significantly change their properties
(e.g., Lugaz et al., 2017; Manchester et al., 2017). The compression regions
and shock waves created by sufficiently fast CMEs (e.g., Gopalswamy, 2006)
play an important role for generating gradual solar energetic particle (SEP)
events (e.g., Kahler, 1992; Desai & Giacalone, 2016). Upon their arrival at
Earth, CMEs can trigger geomagnetic storms as they interact with the
terrestrial magnetosphere (e.g., Pulkkinen, 2007; Zhang et al., 2021). The
resulting geoeffectiveness depends primarily on the CME’s magnetic
configuration, velocity, and ram pressure (e.g., Srivastava & Venkatakrishnan,
2004; Siscoe & Schwenn, 2006).

Whether or not a CME hits Earth depends first and foremost on the location of
its source region on the Sun, but to a large extent also on its trajectory.
While most CMEs propagate away from the Sun in a more or less radial
direction, their trajectories can sometimes deviate from a radial propagation
by several tens of degrees (e.g., Isavnin et al., 2014), which can lead to
false space-weather alerts. A striking example is the January 7, 2014 event
(e.g., Möstl et al., 2015). This fast CME (radial speed ${\sim}2500$ km s-1)
originated in an active region (AR) close to the disk center and was
accompanied by an X1.2-class flare. Based on its speed and origin, many
forecasters predicted a transit time to Earth of ${\sim}36$ hours and a strong
geomagnetic storm. However, the CME arrived significantly later (after about
49 hours), and it delivered only a glancing blow to Earth’s magnetosphere,
resulting in very minor geomagnetic activity (see Mays et al., 2015).
Therefore, in order to improve our capabilities of predicting space-weather
effects associated with CMEs, one important element is to understand the
circumstances and mechanisms that determine CME trajectories.

The deviation of CMEs from a radial trajectory has been often attributed to
their deflection (a change of trajectory) at the typically stronger magnetic
fields (and Alfvén velocities) of coronal holes (CHs) (e.g., Gopalswamy et
al., 2009; Panasenco et al., 2013). Indeed, CMEs that are launched at high
latitudes, close to a polar CH, are often seen to deflect toward the
heliospheric current sheet (HCS) (e.g., Filippov et al., 2001; Kilpua et al.,
2009), and this also has been found in numerical simulations (e.g., Zuccarello
et al., 2012; Bemporad et al., 2012; Talpeanu et al., 2022). The amount of
deflection seems to depend on several parameters, such as the CH area (width)
and the CME speed (e.g., Xie et al., 2009; Mohamed et al., 2012; Wang et al.,
2020, 2022), and such dependencies have recently been studied systematically
in a series of numerical simulations (Sahade et al., 2020, 2021).

Simple analytical calculations of the restoring forces resulting from the
compression of the ambient magnetic field by an expanding CME have shown that
CMEs should be deflected toward regions of lower magnetic energy density (Shen
et al., 2011), and this idea was successfully tested for a sample of eight
CMEs by Gui et al. (2011), although Sieyra et al. (2020) found the opposite
behavior in four of the 13 cases they studied. Using the forces acting on the
surface of a CME draped by its ambient field, Kay et al. (2013, 2015)
developed a semi-analytical tool (ForeCAT) that models the deflection (and
rotation) of CMEs during their propagation.

Non-radial CME trajectories have been also attributed to the presence of
asymmetric magnetic fields in the CME source region, which lead to an
immediate “channeling” or “asymmetric expansion” of the erupting flux (e.g.,
Liewer et al., 2015; Wang et al., 2015; Sahade et al., 2022). Indeed, a pre-
eruptive magnetic flux rope (MFR) that is embedded in an asymmetric field
adopts a tilted configuration (see, e.g., Fig. 4(c) in Titov et al., 2014). To
first order, it can be expected that the erupting MFR will rise in the
direction given by the tilt angle, as in the simulation shown in Fig. 8 in
Török et al. (2013). Such initial channeling was likely the main contributor
to the strongly non-radial trajectory of the January 7, 2014 CME, since the
source-region field was highly asymmetric (Möstl et al., 2015). Of course,
CMEs that experience an initial non-radial rise may still be subjected to
deflection as they travel away from their source region, sometimes resulting
in a reversal of the propagation direction (e.g., Sahade et al., 2023).

These two mechanisms (channeling and deflection) appear to be the main
contributors to the deviation of CME trajectories from a radial direction, as
such deviations seem to take place predominantly in the corona, often below
5–10 $R_{\odot}$ (Isavnin et al., 2014; Kay & Opher, 2015). In interplanetary
space, deviations from a radial propagation direction occur primarily due to
the interaction with other CMEs (e.g., Xiong et al., 2009; Lugaz et al.,
2012), but also via (the presumably relatively weak effect of) flux pile-up
due to the presence of the Parker spiral (e.g., Gosling et al., 1987; Wang et
al., 2004; Zhuang et al., 2019).

In the investigation presented here we focus on an additional CME deflection
mechanism, which, to the best of our knowledge, has been discussed so far only
by Lugaz et al. (2011). The mechanism is similar to the one proposed by
Isenberg & Forbes (2007) to explain the rotation of erupting MFRs about their
rise direction in the presence of an external shear field (oriented in the
same direction as the initial MFR axis). The interaction of the axial MFR
current, ${\bf I}_{\mathrm{FR}}$, with the shear field, ${\bf
B}_{\mathrm{EX}}$, produces a Lorentz force, ${\bf F}={\bf
I}_{\mathrm{FR}}\times{\bf B}_{\mathrm{EX}}$, which points in opposite
(horizontal) directions in the two MFR legs. The forces are balanced as long
as the MFR remains in a stable magnetic equilibrium, but start to rotate the
MFR’s top part as soon as an eruption commences. This effect was modeled and
studied systematically in Kliem et al. (2012) and Zhou et al. (2023).

In principle, the same mechanism can be applied to CMEs that propagate in a
unipolar radial field (see Lugaz et al., 2011). This will be the case, for
instance, if the CME is launched below a pseudostreamer (e.g., Török et al.,
2011) or outside the streamer belt (e.g., Liu, 2007). Here the CME takes the
role of the erupting MFR, while the radial field takes the role of the shear
field. For a CME whose legs are connected to the solar surface, we then expect
that its front part will be deflected out of the plane of propagation (since
in the CME legs the axial MFR current will be approximately aligned with the
radial field). We note that this may be analogous to the rotation of
spheromaks in an ambient field, caused by the “spheromak tilting instability”
(see Asvestari et al. 2022 and references therein.)

An essential assumption made in Isenberg & Forbes (2007) is that the
background shear field can freely penetrate the MFR body (see their Fig.
11(b)) and so produce the above Lorentz force, just as in the case of a
current-carrying conductor placed into a vacuum magnetic field. In a highly
conducting solar plasma, however, the interaction between the MFR and the
ambient field has a more intricate and mediated character. Intruding into a
surrounding magnetized plasma, the erupting MFR preserves its internal field
structure (as long as its axial field strength remains larger than the
strength of the ambient field) by forming shielding current layers at its
surface. For an MFR with a net axial (toroidal) current that intrudes into a
unipolar radial field, this leads to a magnetic-pressure difference at the two
flanks of the MFR (see § 4) whose resultant force might, in principle, differ
from the ${\bf I}_{\mathrm{FR}}\times{\bf B}_{\mathrm{EX}}$ force conjectured
by Isenberg & Forbes (2007). However, in the framework of an idealized 2D
model, where the MFR is considered as a solid perfect conductor, Yeh (1983)
showed that this resultant force is formally equal to ${\bf
I}_{\mathrm{FR}}\times{\bf B}_{\mathrm{EX}}$ if the original (unperturbed)
ambient field is uniform, i.e., if the currents producing it are sufficiently
far away from the MFR current, so that magnetic-pressure gradients in that
field can be neglected (otherwise, additional force terms come into play; see
Eq. (15) in Yeh, 1983). It is important to note that the resultant force,
which we refer to as the “effective ${\bf I}\times{\bf B}$ force” henceforth
(to emphasize the formal resemblance to ${\bf I}_{\mathrm{FR}}\times{\bf
B}_{\mathrm{EX}}$), does not require any asymmetry of the original background
field. Consequently, a deflection of CMEs propagating in a radial field by
this force should always be expected, in addition to the channeling/deflection
caused by net restoring forces resulting from the presence of asymmetry in the
original background field, as described above.

We note that, in addition to the approach by Yeh (1983), the surrounding
magnetized plasma should be considered as a perfect conductor as well. This
implies the formation of an additional current layer that tries to shield the
background medium from the field produced by the MFR current, thereby
decreasing the effective ${\bf I}\times{\bf B}$ force. However, as we shall
see in § 4, this effect is relatively weak, at least during the initial
propagation of the CME, where most of its deflection takes place.

The goal of this paper is to investigate the deflection of CMEs by the
effective ${\bf I}\times{\bf B}$ force in a systematic manner, using
magnetohydrodynamic (MHD) simulations. We first consider, as a “proof of
concept”, a very idealized configuration in which the initial MFR field and
background field are superimposed (§ 2). We then employ a more realistic
setup, in which CMEs are launched from bipolar ARs low in the corona before
they encounter open-field regions, so that no initial superposition of the MFR
and the with the large-scale, radial background field is present (§ 3). The
results of this second investigation are presented in § 4, which is followed
by a discussion and our conclusions in § 5.

Figure 1:  Simplified “proof-of-concept” numerical experiment. (a) Freely
expanding toroidal MFR in a Cartesian box (side view, perpendicular to the
torus axis). Shown are the initial state (left) and the MFR during its
expansion (right). Rainbow-colored (purple) field lines show the MFR core (the
field produced by the MFR current). The symbols denote the direction of the
axial current. Red/blue colors show the horizontal component of the Lorentz
force density, $\bf{j}\times\bf{B}$, in two MFR cross-sections. (b) Same
initial MFR as in (a), now with a uniform vertical magnetic field, $B_{z}>0$,
superimposed. This leads to a rotation of the MFR during its expansion.

##  2 Simplified “Proof-Of-Concept” Experiment

Before we discuss our CME simulations in §3 and §4, we briefly describe here a
simplified “proof-of-concept” experiment that we used to test the effect of a
uniform ambient field on the evolution of an expanding MFR for a case where
the MFR field and the ambient field are superimposed, as in Isenberg & Forbes
(2007). This experiment is shown in Figure 1.

We consider the simplest possible configuration, namely a freely expanding
current ring (a toroidal MFR) in a homogeneous background field. The
simulation is performed in a Cartesian domain with the MHD code described in
Török & Kliem (2003). We use the $\beta=0$ MHD approximation, where thermal
pressure and gravity are neglected and the system is driven merely by inertial
and Lorentz forces. The current ring is modeled by modifying the
Titov–Démoulin MFR model (TD; Titov & Démoulin, 1999), yielding thus the so-
called “modified Titov–Démoulin” model (TDm; Titov et al., 2014). In contrast
to the TD model, which describes a partial current ring (a line-tied coronal
MFR), we consider here a full current ring, and we remove the external fields
that are present in the TD model. This means that the MFR is not in
equilibrium and starts to expand right away at the beginning of the
simulation. We place the MFR here such that the symmetry axis of the torus is
parallel to the $x$ axis.

Figure 1(a) shows a case without any external field. The left panel depicts
the initial configuration. The inset shows that the horizontal Lorentz force
densities along the symmetry axis of the torus are evenly distributed within
the MFR. Consequently, as the MFR expands (right panel), it fully remains in
its initial plane. Figure 1(b) shows a different experiment, where we add a
uniform vertical field to the system. In the case shown, the strength of that
field is 10% of the maximum field strength within the MFR. We can see in the
left panel that the initial horizontal Lorentz force densities are now
oppositely directed in the top and bottom parts of the MFR, and that they
point into the direction we expect from crossing the axial MFR current with
the upward directed background field. Indeed, as the MFR expands, we now see a
significant rotation (right panel), as expected from the considerations in
Isenberg & Forbes (2007).

Figure 2:  Global background magnetic configuration for our CME runs. (a)
Initial configuration (global dipole). The sphere shows the solar surface,
colored by the radial magnetic field, $B_{r}$. Only half of the simulation
domain (out to about $15R_{\odot}$) is shown here. Field lines are colored by
$B_{r}/|{\bf B}|$. (b) Configuration after MHD relaxation, during which the
solar wind opens up flux and the heliospheric current sheet (shown here in the
plane $\phi=\pi$ by the quantity $|{\bf j}|/|{\bf B}|$) is formed. (c) Zoom
into the streamer base, with the $(r,\theta)$ numerical mesh depicted in the
plane $\phi=\pi$ (only every third grid point is shown in each direction). The
mesh is finer at the Equator and at 45°North, where an AR will be inserted.

While this experiment demonstrates that the presence of a unipolar ambient
field can indeed lead to a significant rotation of an expanding MFR, we have
to keep in mind two strong simplifications. First, in our experiment the
initial MFR and the vertical background field are superimposed , leading to
strong sideways directed forces within the MFR. In a real CME, the pre-
eruptive configuration (e.g., a line-tied MFR) is embedded in a closed,
locally bipolar ambient field, which is often referred to as the “strapping
field”. Before the erupting MFR can encounter a unipolar radial field, it
first has to overcome the strapping field, which may happen by pushing it to
the side (e.g., Amari et al., 1996), via “breakout” reconnection above the MFR
(e.g., Antiochos et al., 1999; Lynch et al., 2008), or via “tether-cutting”
reconnection below the MFR that transfers strapping-field flux into flux of
the MFR (e.g., Moore et al., 2001). Even then, as described in § 1 and
illustrated in Figure 6 below, the background radial field cannot freely
penetrate the MFR’s magnetic field, due to the formation of shielding
currents. Second, during the eruption the MFR will strongly expand, so its
current and Lorentz force densities will be less concentrated than in our
simplified experiment, where we inserted a compact MFR directly into a
unipolar background field. Thus, in order to see whether a significant
deflection still can take place after the initial MFR expansion and without a
penetration of the MFR by the background field, we have to consider a model
that includes the initial eruption that produces a CME. The model we are using
for this task is described in the next section.

Figure 3:  Initial configuration for Runs 1, 3, 5, and 7. The sphere shows
$B_{r}$ at the solar surface; north/south is at the top/bottom. In Runs 1 and
3, the AR is placed below the streamer belt and HCS (cf. Figure 2(b)), and in
Runs 5 and 7 within an area of unipolar open field at 45°north. Green (closed)
and magenta (open) field lines are shown, together with electric currents
visualized by $|{\bf j}|/|{\bf B}|$) in the plane $\phi=\pi$ that cuts through
the AR center. The insets show the AR magnetogram and the MFR (orange field
lines) in a top-down view, where east/west is on the left/right. The MFR axial
current (black arrows) is oppositely directed in the configurations on the
left and right, respectively, because of the opposite orientations of the AR
main polarities. All MFRs shown are right-handed (positive helicity).

##  3 Numerical Setup

Our simulations are similar to those described in Lionello et al. (2013). They
are performed with the Magnetohydrodynamic Algorithm Outside a Sphere (MAS)
code, which solves the 3D time-dependent resistive MHD equations on a non-
uniform spherical mesh (see Török et al., 2018, and references therein).

Here we use MAS in the polytropic approximation (Linker et al., 1999) with a
fixed polytropic index of $\gamma=1.05$. For simplicity, and to better isolate
the impact of the effective ${\bf I}\times{\bf B}$ force, we do not consider
additional waves that are able to accelerate a fast solar wind. This results
in an overall slow solar wind profile that is roughly uniform in latitude. Our
computational domain extends from $1\,R_{\odot}$ to $30\,R_{\odot}$ and is
discretized on a non-uniform mesh in all three spherical coordinates
($r,\theta,\phi$), using $(405\times 368\times 318)$ points. The radial
resolution of the mesh varies between $0.002\,R_{\odot}$ at the solar surface
to $0.413\,R_{\odot}$ at the outer boundary.

As described below, we use two different locations where we place an AR,
namely either at the equator ($0^{\circ}$ latitude) or at $45^{\circ}$ north.
To resolve the AR, we construct a finer mesh around these two latitudes. The
minimum and maximum resolution of the mesh along longitudes (north–south;
$\theta$) is $0.113\degree$ and $3.98\degree$, respectively, and along
latitudes (east–west; $\phi$) it is $0.114\degree$ and $7.67\degree$,
respectively. The resistivity and viscosity are chosen to be uniform and set
such that the corresponding diffusion times are much larger than the Alfvén
time of ${\approx}\,24$ minutes (by factors $10^{5}$ and $10^{3}$,
respectively).

The background configuration of the corona starts from a global dipole field
that is subjected to an MHD relaxation, during which the solar wind partially
opens up the magnetic field. The resulting configuration (Figure 2) resembles
solar minimum conditions, and at the solar surface the typical strength of the
radial field is a few G. Following the relaxation we insert into this
configuration a bipolar AR that contains an MFR in stable equilibrium, using
the TDm model and two sub-photospheric magnetic charges to create the
strapping field, as in the original TD model. Without loss of generality, we
place the AR center at $\phi=\pi$. The maximum strength of the field in the AR
is about 400 G.

The insertion is done via superposition of the AR field and the background
field. Since the AR field is two orders of magnitude larger than the
background field, the strapping field spatially divides the MFR from the
large-scale radial background field, into which the MFR intrudes only after
its eruption. Note, also, that in the TDm model the MFR and the strapping
field are superimposed in such a way that the initial net force acting on the
MFR approximately vanishes (see Titov et al. 2014 for details). Both aspects
are different from the simple superposition of the MFR and background field in
the configuration described in § 2. For simplicity, since we are only
interested in the CME trajectory here, we do not modify the plasma parameters
when inserting the AR.

We consider eight different configurations, which differ by the location of
the inserted AR (equator or $45^{\circ}$ north, i.e., under the HCS or in a
unipolar ambient field), its magnetic orientation (north–south or
south–north), and the handedness of the MFR (right-handed or left-handed). In
all eight configurations, the MFR axis is east–west oriented. The odd-numbered
runs (Runs 1, 3, 5, and 7) and even-numbered runs (Runs 2, 4, 6, and 8) differ
only in the direction of the MFR’s axial current with respect to its axial
field (the handedness), yielding a right-handed MFR for odd-numbered runs and
a left-handed MFR for even-numbered runs (corresponding to positive and
negative helicity, respectively). The initial configurations of the four runs
with right-handed MFRs are shown in Figure 3. In Run 1 (top left panel), the
AR is located below the HCS and it has the same magnetic orientation
(north–south) as the global background field, resulting in a bipolar source
region (see the closed green field lines). In Run 3 (top right panel), the AR
orientation is reversed, resulting in a quadrupolar source region. In Run 5
(bottom left panel), the AR is inserted at $45^{\circ}$ north with the same
magnetic orientation as in Run 1, resulting in a bipolar source region in
which the closed field at larger heights is tilted towards the equator.
Finally, in Run 7 (bottom right panel) the AR orientation is reversed again,
resulting in a pseudostreamer-like source region, where the MFR is located in
the southern lobe of the pseudostreamer.

We note that for these simulations we have specifically chosen to insert each
combination of AR and MFR into the exact same background MHD solution to
initiate the CME. Alternatively, we could have run eight separate MHD
relaxations using the surface flux-distribution that results from the
superposition of the AR, MFR, and global dipole for each case. Then, the
current-carrying part of the MFR field could be embedded into the background
configuration to initiate the CME without modifying $B_{r}$ at the surface
(e.g., Török et al., 2018; Titov et al., 2022). While the latter approach is
desired for case-study modeling, doing so would yield a slight displacement of
the HCS away from the equatorial plane in the relaxed background solution,
which would formally be different for each case. The latitudinal displacement
of the HCS would be small for Runs 1–4, but on the order of 8–14∘ for Runs 5–8
(as estimated from PFSS extrapolations using the POT3D code; Caplan et al.,
2021). This displacement would also be oppositely directed for Runs 5/6
(northward) and 7/8 (southward), respectively. Therefore, because our goal is
to study the effective ${\bf I}\times{\bf B}$ force in a background that is
both unipolar and as similar as possible from run to run, we opted for the
simpler setup, where the respective AR-MFR configurations are inserted into
the same dipolar global background.

For each of the eight configurations we perform a number of short test runs to
determine the maximum MFR current that still yields a stable configuration.
Once this value is determined, we trigger an eruption in each production run
by setting the MFR current slightly (about $1\%$) above the stable-equilibrium
value. This ensures that the eruption starts slowly and then accelerates, as
observed for most CMEs (e.g., Vršnak, 2001). Each simulation is run for 50
Alfvén times (about 20 hours). Figure 4 shows, as an example, a snapshot of
the erupting MFR (or CME) for Run 7.

Figure 4:  Erupting MFR after 3 Alfvén times (72 minutes) in Run 7. Field
lines are colored by the normalized transverse magnetic field
$(B_{\theta}^{2}+B_{\phi}^{2})^{1/2}/|\bf B|$ (denoted here as Btang/B).
Electric currents are shown by $|{\bf j}|/|{\bf B}|$ in the $\phi=\pi$ plane,
outlining the HCS and the shock wave produced by the CME.  Figure 5:
Distribution of tracer particles in and above the initial MFR (here for Run 7;
cf. corresponding inset in Figure 3). For clarity, only selected tracers are
displayed (as small spheres). Shown are 19 of the 1000 arc (yellow) and axis
(magenta) tracers, and three of the 49 cross-sections (nos. 8, 25, and 42; in
cyan) with 129 tracers each. Field lines originating at the tracers are shown
for the arc and the axis. The transparent gray iso-surface shows
$j=0.025\,j_{\mathrm{max}}$, outlining the edge of the MFR.  Figure 6:  Sketch
illustrating the effective ${\bf I}\times{\bf B}$ force introduced in § 1. (a)
Unipolar background field. (b) The intrusion of a perfectly conducting MFR
i.e., $I=0$, with a fully balancing “direct” and “return” axial current, not
shown) leads to the formation of shielding surface currents (SC; shown in red
and blue, where red depicts currents directed into the plane and blue depicts
currents directed towards the viewer) and a deformation of the background
field. (c) An MFR with a direct current ($I>0$) additionally produces a
magnetic-pressure difference in the ambient field, yielding a sideways
deflection force. (d) Most general case with a return current (RC; light blue)
and a direct current (DC), where the RC partially balances the DC. See text
for further details.

At the beginning of each simulation, we distribute “tracer particles” (fluid
elements) within and slightly above the MFR, and follow their trajectories
over the course of the simulation. To do so, we interpolate, at each time-step
and for each tracer, the velocities from the adjacent grid points to the
tracer location, and advect the tracers according to the resulting velocity
vector and current time step.

Figure 5 illustrates how the tracers are distributed initially. A first group
of tracers is distributed evenly along the MFR axis. A second group consists
of “projections” of each axis tracer along the radial direction, placed at a
distance of 1.5 times the MFR minor radius away from the axis. This “arc” of
tracers represents the closed arcade field located just above the MFR.
Finally, a third group is distributed within evenly spaced cross-sections of
the MFR (perpendicular to the axis). We place $1,000$ axis tracer particles,
$1,000$ arc tracer particles and $129$ tracer particles at each of the MFR
cross-sections (there are a total of $49$ such cross-sections). For our
analysis described in the next section, we only include tracer particles that
propagated in the radial direction at least 10% above their initial position.
We found that on average (for Runs 5–8, which are the runs that exhibit a
deviation from a radial propagation; see below) $635$ axis and $436$ arc
tracer particles fulfilled this condition, and all of these particles erupted
as part of the CME. Therefore, we included five central cross-sections (within
the vicinity of the MFR apex) in the analysis. This resulted in a similar
number of $645$ cross-section tracer particles, which all erupted as well.

##  4 Results

###  4.1 Effective ${\bf I}\times{\bf B}$ Force and Current Distribution

Before discussing our simulation results, let us first consider in more detail
the effective ${\bf I}\times{\bf B}$ force introduced in § 1. For our
considerations, we assume the following properties of a CME propagating in a
unipolar radial background field. First, the CME-MFR is sufficiently elongated
for its magnetic structure to be approximately invariant along the axial
direction. Second, over the course of its evolution, the MFR tends to preserve
its integrity and approximately circular cross-sectional shape. Third, all
currents remain concentrated within the MFR and its immediate neighborhood, so
that the ambient magnetic field remains practically current-free. Fourth,
adjacent magnetic field lines are largely tangential to the current-carrying
volume. We note that these properties, which in large part are present in our
simulations, allow one to determine the forces acting on the CME-MFR by
integrating the Maxwell stress tensor over the surface of the current-carrying
volume (rather than evaluating the Lorentz forces). We leave a quantitative
analysis of the forces to a later study—here we merely employ the above
properties to discuss the nature of the effective ${\bf I}\times{\bf B}$ force
in a qualitative manner, as illustrated in Figure 6.

Figure 7:  Distribution of electric currents in the $\phi=\pi$ plane for Runs
1, 3, 5, and 7 after 5 Alfvén times ($\approx$ 2 hours), visualized by
$j_{p}\,r^{2}$ (in code units), where $j_{p}$ is the current density along the
$\phi$ direction (which points into the plane) and $r$ is the (sun-centered)
radial coordinate. Hence, red colors show currents directed into the plane,
while blue colors show colors directed towards the viewer. “HCS”, “DC”, “RC”,
and “SC” stand for the heliospheric current sheet, direct (axial) current,
return current, and shielding current, respectively.

Consider, first, a unipolar background magnetic field as shown in panel (a).
The intrusion of an MFR (in the simplest case, a perfect solid conductor,
i.e., a fully current-neutralized MFR) into such a field implies the formation
of oppositely directed surface shielding currents (SC) that prevent the
background field from penetrating the MFR and force it to “flow” around it
(panel (b)). However, if the MFR carries a direct current (DC), the poloidal
(azimuthal) magnetic field produced by that current additionally contributes
to the ambient field, increasing the tangential component of the ambient field
on one side of the MFR and decreasing it on the other (panel (c)). This
results in a sustained difference in magnetic pressure on the two flanks of
the MFR, i.e., in a force that pushes the MFR sideways, with the oppositely
directed field lines on the right-hand side facilitating magnetic reconnection
(see § 5). We refer to this sideward motion as “deflection” henceforth, but
note that we do not mean here a change of direction due to a collision with
some other structure, such as a CH. The direction of the deflection depends on
the respective orientations of the initial background field and the MFR axial
current.

We note that the case shown in panel (c) has been considered earlier by Yeh
(1983), who demonstrated that, for a uniform background field ${\bf B_{0}}$
and MFR axial current $I$, the corresponding net force exactly coincides with
${\bf I}\times{\bf B_{0}}$, as if $B_{0}$ would freely penetrate the MFR body
as assumed by Isenberg & Forbes (2007) (see § 1). Our simulations described
below suggest a more general case, which is sketched in panel (d) and includes
a return current (RC) within the MFR. The presence of such a current can be
understood if one considers that the magnetized plasma surrounding the MFR is
an almost perfect fluid conductor, trying to prevent the azimuthal MFR field
from penetrating the external volume by generating an RC. Note that the
presence of an RC weakens the deflection force compared to the case shown in
panel (c), but that force should still be present as long as the RC remains
smaller than the direct current (DC) that flows in the opposite direction in
the core of the MFR.

Figure 7 shows the distributions of electric currents that develop during the
eruption for Runs 1, 3, 5, and 7 (i.e., runs with a right-handed MFR). It can
be seen that the DC, RC, and SC sketched in Figure 6(d) are present in each
case. There are additional currents present in the system, predominantly at
the wave/shock-front and in the wake of the CME, but those do not seem to be
relevant for the CME deflection considered here. The sign of the DC in the MFR
core (corresponding to $I_{\mathrm{axial}}$ in Figure 3) is determined by the
magnetic orientation of the AR (see insets in Figure 3). It is, in all cases,
surrounded by an RC of opposite sign near the periphery of the CME-MFR. The
RC, in turn, is enclosed by the SC, whose sign is determined by the
orientation of the ambient field. In Runs 1 and 3, where that field is
oppositely directed on the two sides of the HCS, the SC has the same
(positive) sign at both flanks of the CME. In Runs 5 and 7, on the other hand,
the displaced radial field has the same sign at both flanks of the CME, so the
SC has opposite signs at the respective CME flanks. The overall current
distributions in the simulations are in line with the theoretical
considerations outlined above. Hence, we expect to find a deflection of the
CME by the effective ${\bf I}\times{\bf B}$ force in those runs where the CME
propagates in a unipolar ambient field.

Figure 8 shows the CMEs for Runs 1, 3, 5, and 7, about ten hours after the CME
is launched and before its front reaches the outer simulation boundary at
$30\,R_{\odot}$. We can see that the CMEs in Runs 1 and 3, which are launched
under the HCS, do not exhibit any signatures of deviation from a radial
propagation. This is expected, since the effective $\bf{I}\times\bf{B}$ force
points into opposite directions on the two sides of the HCS (due to the
oppositely-directed radial background fields), so that the overall net effect
is zero.

In contrast, Runs 5 and 7, where the CME propagates into a purely unipolar (in
this case, positive) radial field, exhibit a significant deviation from a
radial propagation direction (indicated by the dashed black line). The
respective directions of the deflection (southward for Run 5, northward for
Run 7) are as expected from the effective $\bf{I}\times\bf{B}$ force for the
respective orientations of the MFR axial currents (see Figure 3).

Figure 8:  CMEs for Runs 1, 3, 5, and 7 after 25 Alfvén times ($\approx$10
hours), visualized by $|{\bf j}|/|{\bf B}|$ (in code units) in the $\phi=\pi$
plane. The dashed black line marks the radial direction above the AR. The
radial extension of the simulation domain is $(1-30)R_{\odot}$. Runs 1 and 3
propagate along the HCS, while Runs 5 and 7 propagate in a unipolar radial
ambient field (cf. Figure 3).

As we describe below, our even-numbered runs behave essentially in the same
way as the odd-numbered ones, except for the opposite sense of rotation of the
erupting MFR about its rise direction. Therefore, we focus our analysis of CME
deflection on Runs 5 and 7. Quantifying the deflection is not straightforward
because the CME-MFRs do not expand in a simple, self-similar manner. We
therefore consider the bulk angular CME deflection and employ two independent
methods for estimating its amount.

###  4.2 Tracer Particles

First, we employ our tracer particles to estimate the angular CME deflection.
The first three panels in Figure 9 depict, for comparison, the angular
deviations in spherical space for the three groups of tracers described in § 3
(for each group, we follow the tracers’ positions until the first tracer hits
the outer boundary at 30 $R_{\odot}$). Shown are the east–west (E–W; along
constant latitudes) and north–south (N–S; along constant longitudes)
deviations of individual tracers (light-blue curves) over time for Run 5 and
Run 7 (top and bottom rows, respectively). The rightmost panel shows the full
set of tracers included in our analysis.

We associate the bulk angular deflection of the CME with the angular
deflection of the center of mass of the full set of tracers. To obtain the
latter, we compute their mean position in Cartesian coordinates, and convert
that position back to spherical coordinates. The red curves in the rightmost
panel of Figure 9 show the resulting evolution of the center of mass over the
course of the simulation. For completeness, we also show the corresponding
values for our three tracer groups separately.

For Run 5, we thus find a total (final) southward and eastward deflection of
17°and 5°, respectively, while for Run 7 we find a total northward and
westward deflection of 13°and 9°, respectively. The bottom panel in Figure 11
below shows the evolution of the CME deflection in the N–S direction over time
for both runs. In Runs 1 and 3, where there is no deflection, the
corresponding angular tracer deviations (not shown here) merely represent the
lateral expansion of the CME, with all mean values being zero.

Figure 9:  Trajectories of tracer particles in the $(\theta,\phi)$ plane (with
respect to their initial positions in that plane) along the N–S
($-\Delta\theta$) and E–W ($\Delta\phi$) directions for Run 5 (top) and Run 7
(bottom). The mean trajectories, obtained from a center-of-mass calculation
(see text for details) are shown in red. The first three panels show,
respectively, tracers starting from the initial MFR axis, the arc above this
axis, and five MFR cross-sections around the initial axis apex; the rightmost
panel shows all tracers combined. The annotated values are the mean angular
deviations at the time when the first tracer of the respective group reaches
the outer boundary of the simulation domain.

We note that the selected tracers represent the evolving CME volume fairly
well, but not completely. Specifically, they lag behind the visible CME
leading edge, which has already left the simulation domain when the first
tracers arrive at $r=30\,R_{\odot}$. This discrepancy is because of the so-
called stand-off distance between the actual CME-MFR (as covered by the
tracers) and the “observed” CME leading edge, i.e., the shock/wave front (see
Ma et al. 2011, Gopalswamy et al. 2012, and Figure 10 below). We also note
that the arc tracers, which are initially located entirely above the MFR,
largely become part of the CME-MFR, presumably due to reconnection occurring
between the expanding MFR and adjacent flux systems (e.g., Aulanier & Dudík,
2019). This justifies the inclusion of those tracers for our deflection
analysis.

Due to the asymmetric distribution of the closed-field regions above the AR in
Runs 5 and 7 (see green field lines in Figure 3), some deflection along the
N-S direction already occurs before the MFR reaches the purely radial field at
a height of less than 2 $R_{\odot}$. However, as can be seen in Figure 11
(bottom), this initial phase amounts only to a relatively small fraction of
the total deflection (it takes place way before the first respective data
point in the figure). At even lower heights, on the scale of the AR, the field
is almost perfectly symmetric, since the AR field is two orders of magnitude
larger than the ambient field. Hence, during the very early phase of the MFR
eruption the rise direction should be almost perfectly radial.

While the deflections in N–S (or S–N) direction are expected from the
respective direction of the effective $\bf{I}\times\bf{B}$ force, the cause of
the (mostly significantly smaller) deflections in E–W (or W–E) direction is
less obvious. They most likely result from the bulk rotation of the top part
of the erupting MFR about its rise direction (see Figure 4), which locally
changes the orientation of the axial MFR current and, hence, the direction of
the effective $\bf{I}\times\bf{B}$ force. Indeed, the eastward (westward)
deflection in Run 5 (Run 7) is in line with what one would expect from the
sense of rotation (clockwise for right-handed MFRs; e.g., Green et al., 2007;
Lynch et al., 2009; Kliem et al., 2012) and the resulting change of direction
of the effective $\bf{I}\times\bf{B}$ force. Moreover, for our Runs 6 and 8
(which differ from Runs 5 and 7 only in the MFR handedness) we find the same
southward and northward, but exactly opposite westward and eastward
deflection, respectively, which concurs with the opposite sense of rotation
(counter-clockwise) in those runs. Such rotations about the rise direction are
expected for erupting MFRs even in the absence of an external shear field, in
which case they are caused by the conversion of MFR twist into writhe (e.g.,
Török et al., 2010).

Figure 10:  Example of GCS fitting results for (a–f) Run 5 and (g–l) Run 7.
Panels (a–c) and (g–h) show base-difference synthetic coronagraph images
(where the frame at $t=0.0$ hours is subtracted from the frame at $t=8.8$
hours) from the simulated L5, L1, and L4 points, while panels (d–f) and (j–l)
show the same images with the GCS wireframe projection (in red) overlaid.

###  4.3 GCS Modeling

The second method we use for determining the CME deflection is based on a
technique that is widely employed to characterize CME propagation through the
solar corona for real events, namely forward-modeling via the Graduated
Cylindrical Shell (GCS; Thernisien, 2011) model. In brief, the GCS model
consists of a parametrized (hollow) shell, with a half-toroidal main body and
two conical legs connected back to the Sun (reminiscent of an MFR morphology).
Forward-modeling with the GCS technique is performed by visually fitting the
shell to simultaneous coronagraph images, while tuning its six free parameters
($\theta$, $\phi$, nose height, axial tilt, center-to-leg half-width, and
aspect ratio) until the wireframe’s projection onto each plane of sky best
matches the observational data. This technique has been employed to evaluate
CME deflection through the corona in a number of studies (e.g., Gui et al.,
2011; Kay et al., 2017). In order to apply the GCS model to our simulated
CMEs, synthetic coronagraph images need to be produced. Theoretically, there
are no constraints on where to place synthetic observers for the generation of
such images but, in order to mimic a likely real-event situation, we choose to
set three viewing perspectives at a heliocentric distance of 1 au and along
the solar equatorial plane. One observer is placed at $\phi=180^{\circ}$, i.e.
aligned in longitude with the source AR, and the remaining two are placed at
$\phi=120^{\circ}$ and $\phi=240^{\circ}$, simulating the Sun–Earth L1, L5,
and L4 points, respectively. SOHO and the twin STEREO spacecraft (all three
equipped with coronagraphs) were in this approximate configuration in October
2009; furthermore, the L4 and L5 points have gained considerable attention in
recent years because of their potential benefits for space-weather forecasting
(e.g., Vourlidas, 2015; Bemporad, 2021; Posner et al., 2021). Once the
synthetic observers have been set, we produce corresponding coronagraph images
(i.e., total K-corona brightness maps; see Mikić & Linker 1996 and Howard &
Tappin 2009 for a detailed description of the procedure) throughout the
simulation’s temporal (${\sim}20$ hours) and spatial ($[1$–$30]\,R_{\odot}$)
domains.

We perform GCS fittings of the CME for both Run 5 and Run 7, at regular
intervals (in time) until the CME’s nose/apex exits the coronagraph’s field of
view in at least one data set, for a total of 10 fits. An example of such
fits, performed at $t=8.8$ hours, is shown in Figure 10. It is evident from
these snapshots that the CMEs in Run 5 and Run 7 have already deflected
considerably southward and northward, respectively, at this point in the
simulation (cf. Figure 8). By $t=9.8$ hours and at a GCS estimated apex height
of ${\sim}29\,R_{\odot}$, just before the CME leading edge reaches the outer
boundary, we find a total deflection in latitude of $17^{\circ}$ southward for
Run 5 and $15^{\circ}$ northward for Run 7 with the GCS method.

As for the E–W direction, our GCS fitting results do not exhibit significant
deflections, with values within ${\pm}2^{\circ}$ of the AR longitude
($\phi=180^{\circ}$). This is different from the tracer-particles method (see
§ 4.2), and may be due to the larger uncertainties that are known to be
associated with estimates of $\phi$ in comparison to $\theta$—Thernisien et
al. (2009) performed a sensitivity analysis of the GCS model parameters,
finding a mean deviation of ${\pm}1.8^{\circ}$ in $\theta$ and
${\pm}4.3^{\circ}$ in $\phi$—indicating that a deflection of just a few
degrees in $\phi$ is unlikely to be recovered reliably via the GCS technique.

###  4.4 Method Comparison & Deflection Evolution

The evolution of the CME deflection across the whole simulation domain is
shown as a function of height for both runs and methods in the top panel of
Figure 11. For comparison, we plot in the bottom panel the complete temporal
evolution as obtained from the tracer-particles method. Here we consider only
the dominant deflection along the N–S direction. The plots show that, for both
runs, the strongest deflection occurs within a few solar radii above the
surface, which is in line with previous observational studies (e.g., Isavnin
et al., 2014; Kay & Opher, 2015). Subsequently, the added deflection
successively decreases, and there are no more significant changes beyond
$\approx 20\,R_{\odot}$. For Run 5 (blue) the GCS and tracer-particles methods
yield consistent deflection angles, except for the first data point, where the
value obtained from the tracer-particle method is slightly outside of the
estimated GCS error bars. For Run 7 the correspondence is not quite as good:
starting with the fourth data point, the GCS method consistently yields a
somewhat larger deflection angle, with a difference of 2° for the final data
point.

Figure 11:  Top: Heliocentric latitude at the position of the CME leading-edge
height (as determined from the GCS model apex) for Run 5 (blue) and Run 7
(orange). The circles and squares denote the latitudes estimated from the GCS
fitting method and the tracer-particles method, respectively, at corresponding
times. The green cross marks the location of the source AR at $45^{\circ}$
north. For the GCS fitting method, the error bars indicate the estimates for
the average deviations in latitude (${\pm}1.8^{\circ}$) and height
(${\pm}0.48\,R_{\odot}$), found by Thernisien et al. (2009). Bottom:
Heliocentric latitude but as a function of time from the tracer-particles
method. The blue and orange lines use all the tracer-particle data (up to 14
hours), while the squares correspond to the data points shown in the top
panel. The green vertical line marks the arrival time of the CME leading edge
at the outer boundary of the simulation domain.  Figure 12:  Contours of
$B^{2}$ (in code units) for Runs 5 (top) and 7 (bottom) at 2 (left) and 5
(right) $R_{\odot}$ for the total magnetic field (including the AR) at the
beginning of the simulation. The inner sphere shows the radial magnetic field
(in gray-scale) at the solar surface. The view is centered on the AR at
45°north. The small magenta (green) ball marks the position of the solar north
pole (its radial projection onto the respective outer sphere).

By inspecting the morphology of the CMEs in both runs in synthetic white-light
data (Figure 10), we see that the CME in Run 5 is characterized by a
relatively symmetric overall structure, while the CME in Run 7 displays
stronger asymmetries, with its main body “bending” northward with respect to
the underlying signature of the current sheet (see also Figure 8). The GCS
shell, however, does not allow for such deformations, making the fitting
procedure more difficult in the case of CMEs with irregular morphologies.
Additionally, since the GCS forward-modeling method is performed “manually”—or
“by eye”—it follows that there are intrinsic uncertainties related to the
specific user’s fitting choices. The effects of the user’s subjectivity on GCS
results have been explored and discussed in Verbeke et al. (2022). Taking
these considerations into account, we deem the results obtained with the
tracer-particles method to be more robust.

Looking again at Figure 11, we see that the respective amounts of deflection
differ significantly between Run 5 and Run 7. To understand this, we have to
consider the influence of the global ambient field on the CME trajectory.
Figure 12 shows the initial (i.e., after the AR has been inserted)
distribution of the (normalized) magnetic energy density, $B^{2}$, of the
total magnetic field at two representative coronal heights. For Run 5 (top
panels), $B^{2}$ decreases monotonically from the north pole to the equator at
both heights, i.e., the net restoring force acting due to the compression of
the ambient field by the CME (Shen et al., 2011) always points southward,
supporting the southward directed effective $\bf{I}\times\bf{B}$ force at all
times. For Run 7, however, the effective $\bf{I}\times\bf{B}$ force is
northward directed, so the net restoring force acts against it, at least at
coronal heights where the global field dominates (bottom right panel). For
lower heights (bottom left panel), the situation is more complicated than in
Run 5, due to the opposite magnetic orientation of the AR. Here $B^{2}$ has a
local minimum at a latitude higher than that of the AR, so initially the
restoring force should support the northward deflection by the effective
$\bf{I}\times\bf{B}$ force. Overall, however, the counteracting effect of the
net restoring force dominates, leading to a smaller total deflection than in
Run 5.

##  5 Summary and Discussion

We used MHD simulations to systematically investigate the deflection of CMEs
in a unipolar coronal magnetic field, focusing here on the deflection caused
by the “effective $\mathbf{I}\times\mathbf{B}$ force” described in §§ 1 and 4.
Our “proof of concept” simulation described in § 2 showed that a freely
expanding MFR remains in its initial plane, but the addition (superposition)
of a uniform vertical magnetic field leads to a rotation of the MFR during its
expansion, as suggested by Isenberg & Forbes (2007). However, this simulation
ignores one important condition for real CMEs, namely that the ambient field
cannot freely penetrate the CME-MFR. Moreover, the axial MFR current remains
much more compact than in a real CME, where the MFR undergoes a strong
expansion.

We addressed these limitations by constructing a global coronal configuration
that resembles solar-minimum conditions, consisting of a purely unipolar field
in each hemisphere and an HCS in the equatorial plane. Into this background
configuration we inserted a bipolar AR that contained an MFR. We created eight
configurations that differed in the location of the AR (under the HCS or the
radial unipolar field at 45°,N), its magnetic orientation (N–S or S–N), and
the MFR handedness (right- or left-handed). For each configuration we
initiated an eruption of the MFR and quantified the CME deflection using two
independent methods, namely a set of tracer particles (fluid elements) and GCS
modeling.

We found that CMEs launched under the HCS did not exhibit any deviation from a
radial propagation, due to the cancellation of deflection forces across the
HCS. On the other hand, CMEs launched into the (in our case positive) unipolar
radial field exhibited a significant southward or northward angular
deflection, depending on the AR orientation, i.e., the direction of the axial
MFR current. While the (relatively modest) asymmetry of the closed field
overlying the initial MFR likely leads to some channeling (i.e., additional
deflection) of the MFR in the early phase of the eruption, our results
indicate that the bulk of the CME’s deflection is accrued during its
propagation in the open radial field.

Both the GCS technique and the tracer-particles method estimate the total
southward deflection (Run 5) to be $17^{\circ}$. For the total northward
deflection (Run 7) the GCS technique estimate was $15^{\circ}$, while the
tracer-particles method gave a somewhat lower estimate of
$\mathbf{13^{\circ}}$. We attributed the latter difference to limitations of
the GCS method for asymmetric CME geometries. As discussed in § 4.4, the fact
that the northward deflection is smaller than the southward one is consistent
with the northward-pointing magnetic pressure gradient of the large-scale
coronal background field. The resulting, southward-directed net restoring
force (due to the compression of that field by the CME) increases the
deflection of the southward moving CME, whereas it counteracts the deflection
of the northward moving one. The substantial total northward deflection in Run
7 suggests that the effective ${\bf I}\times{\bf B}$ force is significantly
stronger than the net restoring force for the configurations we considered
here. Overall, our results are in line with those found by Lugaz et al.
(2011), who reported a CME deflection by the effective ${\bf I}\times{\bf B}$
force of $10^{\circ}$ (at $5.0\,R_{\odot}$) and $11.5^{\circ}$ (at
$8.3\,R_{\odot}$) in an MHD simulation of an observed event.

For the CMEs that propagate in the unipolar radial ambient field we also found
a (smaller) deflection in the E–W (or W–E) direction, which we attributed to
changes of the orientation of the axial current due to the bulk rotation of
the erupting MFR about its rise direction. For MFRs with the same initial
orientation of the axial current but opposite handedness (i.e., opposite sense
of rotation), we found the exactly same amount, but opposite direction of that
deflection. This indicates that the initial MFR twist, a main factor in
determining the amount of rotation, influences the partition of the deflection
into longitudinal and latitudinal components. The CMEs launched under the HCS,
while undergoing a rotation as well, do not exhibit any eastward or westward
deflection due to force cancellation across the HCS.

We note that the expansion and propagation of the CMEs in our simulations is
accompanied by continuous reconnection between the CME-MFR and the ambient
field. This occurs, for example, between the axial field of the MFR legs and
the ambient field in the form of “interchange reconnection” (e.g., Reinard &
Fisk, 2004; Lugaz et al., 2011; Crooker & Owens, 2012). For our deflected
CMEs, which propagate in a unipolar radial ambient field, interchange
reconnection is dominant in the CME-MFR leg in which the MFR axial field is
largely antiparallel to the ambient field. Over time, this leads to the
“opening” of a considerable fraction of the MFR flux. This partially changes
the path of the axial MFR current and should, therefore, have some effect on
the CME deflection by the effective ${\bf I}\times{\bf B}$ force. We expect
reconnection also to take place between the azimuthal MFR flux and the ambient
field via “breakthrough reconnection” (Titov et al., 2022), which is analogous
to “breakout reconnection” (Antiochos et al., 1999), except that it can occur
in bipolar magnetic configurations. As described in Titov et al. (2022),
breakthrough reconnection takes place at quasi-separatrix layers (QSLs; Priest
& Démoulin, 1995; Demoulin et al., 1996). In our cases of CME propagation in a
unipolar ambient field, this QSL is expected to form on the flank of the CME-
MFR, where the MFR azimuthal field is directed opposite to the ambient field,
i.e., on the right-hand side of the MFR in Figure 6(c,d). Reconnection on this
side of the MFR then “peels away” adjacent ambient field, allowing the
compressed flux on the other side of the MFR to push the CME more efficiently.
Quantifying the exact locations of reconnection, the respective amounts of
magnetic fluxes transferred, and the role of breakthrough reconnection for the
CME deflection, are nontrivial tasks, which we leave for future work.

As shown in Figure 11, the deflection continuously decelerates as the CME
moves away from the Sun, until there is no more notable increase beyond
$\approx 20\,R_{\odot}$. The reason for this saturation of the CME deflection
is not obvious. It may take place, for instance, because other forces such as
the aerodynamic drag force or the magnetic tension force in the MFR legs start
to counterbalance the effective ${\bf I}\times{\bf B}$ force and/or because
the strength of the MFR return current with respect to the direct current
increases over time until it fully neutralizes the latter, so that the
effective ${\bf I}\times{\bf B}$ force vanishes. Such effects are not
straightforward to quantify, as they require a detailed assessment of the
current-evolution and the various forces involved, so we leave also this
question for a future investigation.

Follow-up work will also explore several model parameters that have not been
varied here. For example, we expect the ratio of the strength of the AR/MFR
field to the strength of the large-scale radial field (assumed to be
${\sim}100$ in the pre-eruptive configuration in this study) to have an effect
on the total deflection. Additionally, the orientation of the source-region
polarity inversion line, i.e, of the initial MFR (in our case strictly E–W or
W–E) may play a role, as well as the amount by which the erupting MFR rotates
about its rise direction, which can be controlled to some extent by, e.g., the
initial MFR diameter (Kliem et al., 2012). Finally, future work will examine
more closely the early evolution (0–5 hours) of the erupting MFR, as this
interval corresponds to the range of coronal heights (${<}15\,R_{\odot}$)
where most of the deflection was observed to take place.

The simulations presented in this work demonstrate the significance of the
effective ${\bf I}\times{\bf B}$ force for the development of non-radial CME
trajectories, in addition to initial flux channeling or deflection by
asymmetric source-region or ambient magnetic fields. The effective ${\bf
I}\times{\bf B}$ force should act most efficiently on CMEs propagating in a
largely unipolar radial ambient field. That is, it seems most relevant for
CMEs launched below a pseudostreamer or outside the streamer belt, and
presumably during periods of lower solar activity, when the HCS is less
complex and larger coherent areas of unipolar field should exist at low and
middle latitudes. Its relevance could be tested, for instance, with future
Parker Solar Probe (Fox et al., 2016) in-situ observations of CMEs below
${\sim}20\,R_{\odot}$.

Furthermore, our results could be applied to predictions of CME trajectories
throughout the extended solar corona, with direct implications for the so-
called “CME hit/miss problem” in space weather forecasting (e.g., Wold et al.,
2018; Verbeke et al., 2019). This could be achieved, for example, by
incorporating the effective $\mathbf{I}\times\mathbf{B}$ force (or, more
generally, the concept of shielding currents that develop in response to CME-
MFRs intruding non-uniform ambient fields, as outlined in Yeh 1983) into
(semi-)empirical models used for such predictions.

## Acknowledgments

We thank the anonymous referee for their patience and very detailed and
helpful comments that significantly improved the quality of this article. The
work described here was supported by NASA’s HSR, HSO-Connect, and HTMS
programs (award No. 80NSSC20K1317, 80NSSC20K1285, and 80NSSC20K1274,
respectively) and by NSF’s PREEVENTS and Solar–Terrestrial programs (award No.
ICER-1854790 and AGS-1923377, respectively). T.T. was additionally supported
by Naval Research Laboratory (NRL) contract N0017319C2003 to Predictive
Science Inc., a subcontract of NASA/LWS grant 80HQTR19T0084 to NRL. M.G.L.
acknowledges support by the latter contract and by the Office of Naval
Research.

## References

  * Amari et al. (1996) Amari, T., Luciani, J. F., Aly, J. J., & Tagger, M. 1996, ApJ, 466, L39
  * Antiochos et al. (1999) Antiochos, S. K., DeVore, C. R., & Klimchuk, J. A. 1999, ApJ, 510, 485
  * Asvestari et al. (2022) Asvestari, E., Rindlisbacher, T., Pomoell, J., & Kilpua, E. K. J. 2022, ApJ, 926, 87
  * Aulanier (2014) Aulanier, G. 2014, in Nature of Prominences and their Role in Space Weather, ed. B. Schmieder, J.-M. Malherbe, & S. T. Wu, Vol. 300, 184
  * Aulanier & Dudík (2019) Aulanier, G., & Dudík, J. 2019, A&A, 621, A72
  * Baker & Lanzerotti (2016) Baker, D. N., & Lanzerotti, L. J. 2016, , 84, 166
  * Bemporad (2021) Bemporad, A. 2021, , 8, 627576
  * Bemporad et al. (2012) Bemporad, A., Zuccarello, F. P., Jacobs, C., Mierla, M., & Poedts, S. 2012, Sol. Phys., 281, 223
  * Caplan et al. (2021) Caplan, R. M., Downs, C., Linker, J. A., & Mikic, Z. 2021, ApJ, 915, 44
  * Crooker & Owens (2012) Crooker, N. U., & Owens, M. J. 2012, Space Sci. Rev., 172, 201
  * Demoulin et al. (1996) Demoulin, P., Henoux, J. C., Priest, E. R., & Mandrini, C. H. 1996, A&A, 308, 643
  * Desai & Giacalone (2016) Desai, M., & Giacalone, J. 2016, , 13, 3
  * Filippov et al. (2001) Filippov, B. P., Gopalswamy, N., & Lozhechkin, A. V. 2001, Sol. Phys., 203, 119
  * Forbes (2000) Forbes, T. G. 2000, J. Geophys. Res., 105, 23153
  * Fox et al. (2016) Fox, N. J., Velli, M. C., Bale, S. D., et al. 2016, Space Sci. Rev., 204, 7
  * Gopalswamy (2006) Gopalswamy, N. 2006, , 27, 243
  * Gopalswamy et al. (2009) Gopalswamy, N., Mäkelä, P., Xie, H., Akiyama, S., & Yashiro, S. 2009, J. Geophys. Res., 114, A00A22
  * Gopalswamy et al. (2012) Gopalswamy, N., Nitta, N., Akiyama, S., Mäkelä, P., & Yashiro, S. 2012, ApJ, 744, 72
  * Gosling et al. (1987) Gosling, J. T., Thomsen, M. F., Bame, S. J., & Zwickl, R. D. 1987, J. Geophys. Res., 92, 12399
  * Green et al. (2007) Green, L. M., Kliem, B., Török, T., van Driel-Gesztelyi, L., & Attrill, G. D. R. 2007, Sol. Phys., 246, 365
  * Green et al. (2018) Green, L. M., Török, T., Vršnak, B., Manchester, W., & Veronig, A. 2018, Space Sci. Rev., 214, 46
  * Gui et al. (2011) Gui, B., Shen, C., Wang, Y., et al. 2011, Sol. Phys., 271, 111
  * Howard & Tappin (2009) Howard, T. A., & Tappin, S. J. 2009, Space Sci. Rev., 147, 31
  * Isavnin et al. (2014) Isavnin, A., Vourlidas, A., & Kilpua, E. K. J. 2014, Sol. Phys., 289, 2141
  * Isenberg & Forbes (2007) Isenberg, P. A., & Forbes, T. G. 2007, ApJ, 670, 1453
  * Kahler (1992) Kahler, S. W. 1992, ARA&A, 30, 113
  * Kay et al. (2017) Kay, C., Gopalswamy, N., Xie, H., & Yashiro, S. 2017, Sol. Phys., 292, 78
  * Kay & Opher (2015) Kay, C., & Opher, M. 2015, ApJ, 811, L36
  * Kay et al. (2013) Kay, C., Opher, M., & Evans, R. M. 2013, ApJ, 775, 5
  * Kay et al. (2015) Kay, C., Opher, M., & Evans, R. M. 2015, ApJ, 805, 168
  * Kilpua et al. (2009) Kilpua, E. K. J., Pomoell, J., Vourlidas, A., et al. 2009, , 27, 4491
  * Kliem et al. (2012) Kliem, B., Török, T., & Thompson, W. T. 2012, Sol. Phys., 281, 137
  * Liewer et al. (2015) Liewer, P., Panasenco, O., Vourlidas, A., & Colaninno, R. 2015, Sol. Phys., 290, 3343
  * Linker et al. (1999) Linker, J. A., Mikić, Z., Biesecker, D. A., et al. 1999, J. Geophys. Res., 104, 9809
  * Lionello et al. (2013) Lionello, R., Downs, C., Linker, J. A., et al. 2013, ApJ, 777, 76
  * Liu (2007) Liu, Y. 2007, ApJ, 654, L171
  * Lugaz et al. (2011) Lugaz, N., Downs, C., Shibata, K., et al. 2011, ApJ, 738, 127
  * Lugaz et al. (2012) Lugaz, N., Farrugia, C. J., Davies, J. A., et al. 2012, ApJ, 759, 68
  * Lugaz et al. (2017) Lugaz, N., Temmer, M., Wang, Y., & Farrugia, C. J. 2017, Sol. Phys., 292, 64
  * Lynch et al. (2008) Lynch, B. J., Antiochos, S. K., DeVore, C. R., Luhmann, J. G., & Zurbuchen, T. H. 2008, ApJ, 683, 1192
  * Lynch et al. (2009) Lynch, B. J., Antiochos, S. K., Li, Y., Luhmann, J. G., & DeVore, C. R. 2009, ApJ, 697, 1918
  * Ma et al. (2011) Ma, S., Raymond, J. C., Golub, L., et al. 2011, ApJ, 738, 160
  * Manchester et al. (2017) Manchester, W., Kilpua, E. K. J., Liu, Y. D., et al. 2017, Space Sci. Rev., 212, 1159
  * Mays et al. (2015) Mays, M. L., Thompson, B. J., Jian, L. K., et al. 2015, ApJ, 812, 145
  * Mikić & Linker (1996) Mikić, Z., & Linker, J. A. 1996, in American Institute of Physics Conference Series, Vol. 382, Proceedings of the eigth International solar wind Conference: Solar wind eight, ed. D. Winterhalter, J. T. Gosling, S. R. Habbal, W. S. Kurth, & M. Neugebauer, 104
  * Mohamed et al. (2012) Mohamed, A. A., Gopalswamy, N., Yashiro, S., et al. 2012, , 117, A01103
  * Moore et al. (2001) Moore, R. L., Sterling, A. C., Hudson, H. S., & Lemen, J. R. 2001, ApJ, 552, 833
  * Möstl et al. (2015) Möstl, C., Rollett, T., Frahm, R. A., et al. 2015, , 6, 7135
  * Panasenco et al. (2013) Panasenco, O., Martin, S. F., Velli, M., & Vourlidas, A. 2013, Sol. Phys., 287, 391
  * Posner et al. (2021) Posner, A., Arge, C. N., Staub, J., et al. 2021, , 19, e2021SW002777
  * Priest & Démoulin (1995) Priest, E. R., & Démoulin, P. 1995, J. Geophys. Res., 100, 23443
  * Pulkkinen (2007) Pulkkinen, T. 2007, , 4, 1
  * Reinard & Fisk (2004) Reinard, A. A., & Fisk, L. A. 2004, ApJ, 608, 533
  * Sahade et al. (2021) Sahade, A., Cécere, M., Costa, A., & Cremades, H. 2021, A&A, 652, A111
  * Sahade et al. (2020) Sahade, A., Cécere, M., & Krause, G. 2020, ApJ, 896, 53
  * Sahade et al. (2022) Sahade, A., Cécere, M., Sieyra, M. V., et al. 2022, A&A, 662, A113
  * Sahade et al. (2023) Sahade, A., Vourlidas, A., Balmaceda, L. A., & Cécere, M. 2023, ApJ, 953, 150
  * Shen et al. (2011) Shen, C., Wang, Y., Gui, B., Ye, P., & Wang, S. 2011, Sol. Phys., 269, 389
  * Sieyra et al. (2020) Sieyra, M. V., Cécere, M., Cremades, H., et al. 2020, Sol. Phys., 295, 126
  * Siscoe & Schwenn (2006) Siscoe, G., & Schwenn, R. 2006, Space Sci. Rev., 123, 453
  * Srivastava & Venkatakrishnan (2004) Srivastava, N., & Venkatakrishnan, P. 2004, J. Geophys. Res., 109, A10103
  * Talpeanu et al. (2022) Talpeanu, D. C., Poedts, S., D’Huys, E., Mierla, M., & Richardson, I. G. 2022, A&A, 663, A32
  * Temmer (2021) Temmer, M. 2021, , 18, 4
  * Thernisien (2011) Thernisien, A. 2011, ApJS, 194, 33
  * Thernisien et al. (2009) Thernisien, A., Vourlidas, A., & Howard, R. A. 2009, Sol. Phys., 256, 111
  * Titov & Démoulin (1999) Titov, V. S., & Démoulin, P. 1999, A&A, 351, 707
  * Titov et al. (2022) Titov, V. S., Downs, C., Török, T., & Linker, J. A. 2022, ApJ, 936, 121
  * Titov et al. (2014) Titov, V. S., Török, T., Mikic, Z., & Linker, J. A. 2014, ApJ, 790, 163
  * Török et al. (2010) Török, T., Berger, M. A., & Kliem, B. 2010, A&A, 516, A49
  * Török & Kliem (2003) Török, T., & Kliem, B. 2003, A&A, 406, 1043
  * Török et al. (2011) Török, T., Panasenco, O., Titov, V. S., et al. 2011, ApJ, 739, L63
  * Török et al. (2013) Török, T., Temmer, M., Valori, G., et al. 2013, Sol. Phys., 286, 453
  * Török et al. (2018) Török, T., Downs, C., Linker, J. A., et al. 2018, ApJ, 856, 75
  * Verbeke et al. (2019) Verbeke, C., Mays, M. L., Temmer, M., et al. 2019, , 17, 6
  * Verbeke et al. (2022) Verbeke, C., Mays, M. L., Kay, C., et al. 2022, , in press
  * Vourlidas (2015) Vourlidas, A. 2015, , 13, 197
  * Vršnak (2001) Vršnak, B. 2001, J. Geophys. Res., 106, 25249
  * Wang et al. (2020) Wang, J., Hoeksema, J. T., & Liu, S. 2020, , 125, e2019JA027530
  * Wang et al. (2022) Wang, J., Liu, S., & Luo, B. 2022, , in press
  * Wang et al. (2015) Wang, R., Liu, Y. D., Dai, X., et al. 2015, ApJ, 814, 80
  * Wang et al. (2004) Wang, Y., Shen, C., Wang, S., & Ye, P. 2004, Sol. Phys., 222, 329
  * Wold et al. (2018) Wold, A. M., Mays, M. L., Taktakishvili, A., et al. 2018, , 8, A17
  * Xie et al. (2009) Xie, H., St. Cyr, O. C., Gopalswamy, N., et al. 2009, Sol. Phys., 259, 143
  * Xiong et al. (2009) Xiong, M., Zheng, H., & Wang, S. 2009, J. Geophys. Res., 114, A11101
  * Yeh (1983) Yeh, T. 1983, ApJ, 264, 630
  * Zhang et al. (2021) Zhang, J., Temmer, M., Gopalswamy, N., et al. 2021, , 8, 56
  * Zhou et al. (2023) Zhou, Z., Jiang, C., Yu, X., et al. 2023, , 11, 1119637
  * Zhuang et al. (2019) Zhuang, B., Wang, Y., Hu, Y., et al. 2019, ApJ, 876, 73
  * Zuccarello et al. (2012) Zuccarello, F. P., Bemporad, A., Jacobs, C., et al. 2012, ApJ, 744, 66


††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Intrinsic Nonlinear Planar Hall Effect

Yue-Xin Huang  Research Laboratory for Quantum Materials, Singapore University
of Technology and Design, Singapore 487372, Singapore  Xiaolong Feng  Research
Laboratory for Quantum Materials, Singapore University of Technology and
Design, Singapore 487372, Singapore  Hui Wang  Research Laboratory for Quantum
Materials, Singapore University of Technology and Design, Singapore 487372,
Singapore  Cong Xiao<EMAIL_ADDRESS>Department of Physics, The University of
Hong Kong, Hong Kong, China  HKU-UCAS Joint Institute of Theoretical and
Computational Physics at Hong Kong, China  Shengyuan A. Yang  Research
Laboratory for Quantum Materials, Singapore University of Technology and
Design, Singapore 487372, Singapore

###### Abstract

We propose an intrinsic nonlinear planar Hall effect, which is of band
geometric origin, independent of scattering, and scales with the second order
of electric field and first order of magnetic field. We show that this effect
is less symmetry constrained compared to other nonlinear transport effects and
is supported in a large class of nonmagnetic polar and chiral crystals. Its
characteristic angular dependence provides an effective way to control the
nonlinear output. Combined with first-principles calculations, we evaluate
this effect in the Janus monolayer MoSSe and report experimentally measurable
results. Our work reveals an intrinsic transport effect, which offers a new
tool for material characterization and a new mechanism for nonlinear device
application.

_Intrinsic_ transport properties, which are independent of scattering and
solely manifest the band structure geometry, have been a focus in condensed
matter physics research. As prominent examples, the intrinsic contributions to
the anomalous Hall effect [1, 2, 3], spin Hall effect [4, 5, 6], and anomalous
Nernst effect [7] are connected to the Berry curvature of the bands [8] and
found great success in explaining experimental observations.

In this work, we reveal a new intrinsic transport in the _nonlinear_ planar
Hall effect (NPHE). The planar Hall effect (PHE) refers to the setup where the
applied magnetic $B$ field is inside the transport plane formed by the
transverse response current $j_{\text{H}}$ and the driving $E$ field. This
setup dictates that the response is distinct from the ordinary Hall effect, as
the Lorentz force is not in action here. The linear PHE was initially studied
in magnetic materials [9, 10, 11], as a useful tool for probing the
magnetization reversal. Later, it received interest in nonmagnetic materials,
especially topological materials such as Weyl semimetals [12, 13, 14, 15] and
topological insulators [16, 17, 18]. More recently, NPHE, the nonlinear
version of the effect with $j_{\text{H}}\sim E^{2}B$, was reported in
experiment [19, 20]. Notably, He  _et al._ [19] observed the NPHE at the two-
dimensional (2D) surface of a topological insulator and connected it to the
spin-momentum locking of the topological surface states. A few theories were
put forward to explain the effect [19, 21, 22, 23]. However, it is important
to note that all the NPHE mechanisms proposed so far are of _extrinsic_
nature, i.e., they depend on scattering, which can be readily seen as their
contributions to the the Hall current scale as $\sim\tau^{2}$ in the electron
relaxation time $\tau$.

Here, we propose the _intrinsic_ NPHE, which is completely determined by the
band structure therefore represents an intrinsic material property. We show
that the intrinsic NPHE tensor is expressed in terms of the band geometric
quantities including the Berry-connection polarizability (BCP) [24, 25] and
its spin susceptibility. We clarify the symmetry property of the effect and
find that unlike the linear PHE or nonlinear anomalous Hall effect which are
strongly constrained by symmetry, the intrinsic NPHE can exist in a wide range
of non-centrosymmetric crystals. Moreover, the response exhibits a
characteristic angular dependence in the applied fields, providing a
convenient switch for the nonlinear output signal in applications. As a unique
advantage, intrinsic effects allow a quantitative evaluation. We combine our
theory with first-principles calculations and demonstrate a sizable intrinsic
NPHE in a famous polar 2D material, the monolayer MoSSe. Our work discovers a
new intrinsic transport property, which offers a tool for characterizing
fundamental band geometry of materials and holds great potential for nonlinear
device applications.

Origin of intrinsic NPHE. Consider a 2D nonmagnetic system occupying the
$x$-$y$ plane. First, note that since the current ($E$ field) is odd (even)
under time reversal operation, the intrinsic transport cannot occur without
the applied $B$ field. Second, in the 2D planar Hall setup, the $B$ field
couples only to the electron spin via Zeeman coupling but not to its orbital
motion. Hence, the effect of $B$ field here can be captured as a perturbation
to the band structure, which causes a spin splitting and transforms the
original nonmagnetic band structure to a kind of “magnetic” band structure.

Starting from this perturbed band structure, the intrinsic NPHE of the
original system is equivalent to an intrinsic nonlinear anomalous Hall effect
of the new (perturbed) system. It can be most readily derived in the framework
of semiclassical theory [26, 8]. Like the intrinsic anomalous Hall effect, in
this framework, it is connected to the so-called anomalous velocity of
electrons $\bm{v}^{A}=\bm{E}\times\bm{\Omega}$ (set $e=\hbar=1$) [27, 28],
which is evidently transverse to the $E$ field. Here, $\bm{\Omega}$ is the
Berry curvature, which points out of plane for a 2D system. Note that to study
the nonlinear response at the $E^{2}$ order, one has to include in
$\bm{\Omega}$ the field correction of the Berry curvature $\bm{\Omega}^{E}$
that is linear in $E$. Such a field correction has been captured in the
recently developed extended semiclassical theory [24, 29]. The result shows
that $\bm{\Omega}^{E}=\nabla_{\bm{k}}\times\bm{\mathcal{A}}^{E}$, where
${\mathcal{A}}^{E}_{a}=\tilde{G}_{ab}E_{b}$ is the field corrected Berry
connection, $\tilde{G}_{ab}$ is a gauge-invariant quantity known as the BCP
tensor, the indices $a,b$ label the Cartesian components and repeated indices
are summed over.

The intrinsic nonlinear Hall current is then obtained as
$\bm{j}^{\text{int}}=\int
f_{0}(\tilde{\varepsilon})\bm{E}\times\bm{\Omega}^{E}$, with $f_{0}$ the
equilibrium Fermi distribution function. In this result, the material’s
information is encoded in the band energy $\tilde{\varepsilon}$ and the BCP
$\tilde{G}_{ab}$, but note that here they are defined with respect to the
$B$-perturbed band structure (which is why we add a tilde in these symbols).
By straightforward calculations, we can expand them to the $O(B^{1})$ order
and express all quantities using the original _unperturbed_ eigenstates
$|u_{n}(\bm{k})\rangle$ and band energies $\varepsilon_{n}(\bm{k})$. For
instance, we have

$\tilde{G}_{ab}=G_{ab}+\Lambda_{abc}B_{c},$ (1)

where $G_{ab}$ is the BCP for the unperturbed system, and
$\Lambda_{abc}=\partial_{B_{c}}\tilde{G}_{ab}|_{B=0}$ can be interpreted as
the spin susceptibility of BCP. Explicitly, for a band with index $n$,

$G_{ab}^{n}(\bm{k})=2\mathrm{{\operatorname{Re}}}\sum_{m\neq
n}\frac{v_{a}^{nm}v_{b}^{mn}}{(\varepsilon_{n}-\varepsilon_{m})^{3}},$ (2)

and

$\displaystyle\Lambda_{abc}^{n}(\bm{k})=2\mathrm{{\operatorname{Re}}}\sum_{m\neq
n}\bigg{[}\frac{3v_{a}^{nm}v_{b}^{mn}\left(\mathcal{M}_{c}^{n}-\mathcal{M}_{c}^{m}\right)}{(\varepsilon_{n}-\varepsilon_{m})^{4}}-\sum_{\ell\neq
n}\frac{\left(v_{a}^{\ell m}v_{b}^{mn}+v_{b}^{\ell
m}v_{a}^{mn}\right)\mathcal{M}_{c}^{n\ell}}{\left(\varepsilon_{n}-\varepsilon_{\ell}\right)(\varepsilon_{n}-\varepsilon_{m})^{3}}-\sum_{\ell\neq
m}\frac{\left(v_{a}^{\ell n}v_{b}^{nm}+v_{b}^{\ell
n}v_{a}^{nm}\right)\mathcal{M}_{c}^{m\ell}}{\left(\varepsilon_{m}-\varepsilon_{\ell}\right)(\varepsilon_{n}-\varepsilon_{m})^{3}}\bigg{]},$
(3)

where $v_{a}^{nm}=\langle u_{n}|\hat{v}_{a}\mathbf{|}u_{m}\rangle$ is the
velocity matrix element, $\bm{\mathcal{M}}^{mn}=-g\mu_{B}\bm{s}^{mn}$ is the
spin magnetic moment matrix element, with $g$ the $g$-factor, $\mu_{B}$ the
Bohr magneton, and $\bm{s}^{mn}$ the spin matrix element; and
$\bm{\mathcal{M}}^{n}\equiv\bm{\mathcal{M}}^{nn}$ denotes the diagonal
element. Meanwhile,
$\tilde{\varepsilon}_{n}(\bm{k})={\varepsilon}_{n}-\bm{\mathcal{M}}^{n}\cdot\bm{B}$.
Substituting these expansions into the nonlinear Hall current and keeping
terms up to $O(E^{2}B)$, one obtains the intrinsic NPHE current
$\bm{j}^{\text{int}}$ (details in the Supplemental Material [30]).

The result is cast into an intrinsic NPHE conductivity $\chi^{\text{int}}$
defined by

$j_{a}^{\text{int}}=\chi_{abcd}^{\text{int}}E_{b}E_{c}B_{d},$ (4)

where all the indices $\in\\{x,y\\}$. Then, we find

$\chi_{abcd}^{\text{int}}=\sum_{n}\int\frac{d^{2}\bm{k}}{\left(2\pi\right)^{2}}f_{0}^{\prime}\
\alpha_{abcd}^{n}\left(\bm{k}\right),$ (5)

with

$\alpha_{abcd}^{n}=(v_{a}^{n}\Lambda_{bcd}^{n}-v_{b}^{n}\Lambda_{acd}^{n})+\left(\partial_{a}G_{bc}^{n}-\partial_{b}G_{ac}^{n}\right)\mathcal{M}_{d}^{n}.$
(6)

Evidently, $\chi^{\text{int}}$ is antisymmetric in its first two indices,
which ensures $j_{a}^{\text{int}}E_{a}=0$, so $j_{a}^{\text{int}}$ is truly a
dissipationless Hall current. The expressions in (5, 6) confirm that the
effect is a Fermi surface property determined completely by the material’s
band structure, manifesting its intrinsic nature. In Eq. (6), the first term
of $\alpha$ is an antisymmetrized product of the band velocity $v_{a}^{n}$ and
the spin susceptibility of BCP; whereas the second term represents a product
of the k-space dipole of BCP and the spin magnetic moment.

As emphasized, the intrinsic NPHE $\chi^{\text{int}}$ is independent of
scattering. This is distinct from previously discussed extrinsic contributions
which explicitly involve scattering and scale as $\tau^{2}$ [19, 22, 23] (note
that contribution $\sim\tau E^{2}B$ is forbidden for nonmagnetic systems by
the time reversal symmetry). In practice, the different contributions can be
readily separated by their different $\tau$ scaling.

Symmetry property and angular dependence. Crystal symmetries constrain the
form of the $\chi^{\text{int}}$ tensor. Clearly, $\chi^{\text{int}}$ would
take its simplest form in a coordinate system adapted to the crystal axis [31,
32]. For example, for a rectangular lattice, one would prefer to choose $x$
and $y$ along the $a$ and $b$ axis for the calculation of $\chi^{\text{int}}$.
In real experiments, the applied $E$ and $B$ fields are in general not aligned
with the crystal axis, and one typically rotates these fields in the plane. To
analyze such general cases, we specify $E$ ($B$) field by its polar angle
$\theta$ ($\varphi$) from the $x$ direction (which is already chosen as some
crystal axis), i.e., $(E_{x},E_{y})=E\left(\cos\theta,\sin\theta\right)$ and
$(B_{x},B_{y})=B\left(\cos\varphi,\sin\varphi\right)$, then the intrinsic NPHE
current is in the direction $\hat{z}\times\bm{E}$ with a magnitude

$j^{\text{int}}=\chi^{\text{int}}_{\text{H}}(\theta,\varphi)E^{2}B,$ (7)

where the angle-dependent scalar coefficient $\chi^{\text{int}}_{\text{H}}$ is
a combination of the $\chi^{\text{int}}$ tensor elements:

$\displaystyle\chi^{\text{int}}_{\text{H}}(\theta,\varphi)=$
$\displaystyle(\chi_{yxxx}^{\text{int}}\cos\varphi+\chi_{yxxy}^{\text{int}}\sin\varphi)\cos\theta$
$\displaystyle-(\chi_{xyyx}^{\text{int}}\cos\varphi+\chi_{xyyy}^{\text{int}}\sin\varphi)\sin\theta.$
(8)

Note that since $\chi^{\text{int}}_{abcd}$ is enforced to be antisymmetric in
$a$ and $b$ indices, there are at most four independent tensor elements, which
we choose here as $\chi_{yxxx}^{\text{int}}$, $\chi_{yxxy}^{\text{int}}$,
$\chi_{xyyx}^{\text{int}}$ and $\chi_{xyyy}^{\text{int}}$.

Table 1: Constraint on the $\chi^{\text{int}}$ tensor from point group symmetries pertaining to 2D materials. “$\checkmark$” (“$\times$”) means that the element is symmetry allowed (forbidden). Because $\chi^{\text{int}}$ is time reversal ($\mathcal{T}$) even, symmetry operations $\mathcal{R}$ and $\mathcal{RT}$ impose the same constraint. |  $\mathcal{P}$, $S_{6}$, $\sigma_{z}$ |  $C_{2}^{z}$ |  $S_{4}^{z}$ |  $C_{3,4,6}^{z}$ |  $\ C_{2}^{x},$ $C_{2}^{y}$ |  $\sigma_{x},$ $\sigma_{y}$
---|---|---|---|---|---|---
$\chi_{xyyx}^{\text{int}}$ | $\times$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\times$
$\chi_{yxxx}^{\text{int}}$ | $\times$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\times$ | $\checkmark$
$\chi_{xyyy}^{\text{int}}$ | $\times$ | $\checkmark$ | $\chi_{yxxx}^{\text{int}}$ | $-\chi_{yxxx}^{\text{int}}$ | $\times$ | $\checkmark$
$\chi_{yxxy}^{\text{int}}$ | $\times$ | $\checkmark$ | $-\chi_{xyyx}^{\text{int}}$ | $\chi_{xyyx}^{\text{int}}$ | $\checkmark$ | $\times$

Next, we analyze the symmetry property of the $\chi^{\text{int}}$ tensor.
There are two quick observations. First, the inversion symmetry must be
broken, because the $B$ field is even under inversion whereas the current and
the $E$ field are odd. Second, the system must not have a horizontal mirror
$\sigma_{z}$, because it would flip the $B$ field but keep other vectors in
(4) unchanged. Then, out of the 18 noncentrosymmetric point groups pertaining
to 2D systems, 16 allow the intrinsic NPHE. These include $S_{4}$, $D_{2d}$,
all polar point groups and chiral point groups, i.e., $C_{n}$, $C_{nv}$
($n=1,2,3,4,6$) and $D_{n}$ ($n=2,3,4,6$) [33]. This indicates that the
intrinsic NPHE can be supported in broad crystal systems. In comparison, the
constraint on the linear PHE (with $j_{\text{H}}\sim EB$) is much more
stringent: among these groups, it is only allowed in $C_{1}$, $C_{1v}$, and
$C_{2}$ (in-plane polar axis). Similarly, the nonlinear anomalous Hall effect
by Berry curvature dipole is suppressed in 2D by any symmetry other than a
single mirror line [34].

In Table 1, we list the constraints on the $\chi^{\text{int}}$ elements by
symmetries in the 16 allowed groups. It follows from these results that the
intrinsic NPHE response exhibits rich and interesting angular dependence as
one rotates the $E$ or $B$ fields. For example, polar groups $C_{nv}$ with
$n>2$ support only one independent tensor element
$\chi_{xyyy}^{\text{int}}=-\chi_{yxxx}^{\text{int}}$, and
$\chi_{\text{H}}^{\text{int}}$ is reduced to

$\chi_{\text{H}}^{\text{int}}=\chi_{yxxx}^{\text{int}}\cos(\varphi-\theta),$
(9)

which depends only on the relative angle between the ${E}$ and ${B}$ fields
through a simple cosine relation. The response is maximized when the two
fields are parallel or anti-parallel, but vanishes when they are
perpendicular. In contrast, materials with dihedral chiral point groups
$D_{n}$ with $n>2$ host only
$\chi_{xyyx}^{\text{int}}=\chi_{yxxy}^{\text{int}}$, then

$\chi_{\text{H}}^{\text{int}}=\chi_{yxxy}^{\text{int}}\sin(\varphi-\theta),$
(10)

which is maximized when the two fields are perpendicular, but vanishes when
they become collinear. These cases are illustrated in Fig. 1.

Figure 1: (a) Schematic of the measurement setup for NPHE. (b, c) For a
material with point group ${C_{nv}}$ or ${D_{n}}$ ($n>2$), the intrinsic NPHE
response exhibits a simple cosine or sine dependence on the angle
$(\varphi-\theta)$ between $E$ and $B$ fields, when $B$ field is rotated in
the plane.

Our above analysis demonstrates the intrinsic NPHE as a promising tool for
characterizing a large class of materials not accessible by other nonlinear
transport measurements. The clarified angular dependence offers a convenient
knob to control the nonlinear response and to switch it on and off by simply
rotating the $B$ (or $E$) field (as in Fig. 1(b,c)). This could be useful for
the design of new devices for nonlinear rectification or frequency doubling.

A model study. To illustrate the features of the effect and the underlying
geometric quantities, we first apply our theory to a tight-binding model on a
honeycomb lattice with $C_{3v}$ symmetry. The model reads

$\displaystyle H=$ $\displaystyle-t\sum_{\langle
ij\rangle,\alpha}c_{i\alpha}^{\dagger}c_{j\alpha}+\frac{\Delta}{2}\sum_{i,\alpha}\xi_{i}c^{\dagger}_{i\alpha}c_{i\alpha}$
$\displaystyle+it_{R}\sum_{\langle
ij\rangle,\alpha\beta}(\bm{s}\times{\bm{d}}_{ij})^{z}_{\alpha\beta}c^{\dagger}_{i\alpha}c_{j\beta},$
(11)

where the first term is the nearest neighbor hopping, the second term is a
staggered potential with $\xi_{i}=\pm 1$ on the two sublattices, the last term
is a Rashba spin-orbit coupling (SOC), $\alpha$ and $\beta$ are the spin
labels, and ${\bm{d}}_{ij}$ is a unit vector pointing from site $j$ to $i$.
This model is similar to that for silicene or germanene on a substrate (some
intrinsic SOC terms inessential to our discussion are neglected) [35, 36]. The
staggered potential violates the inversion symmetry and the Rashba SOC breaks
the horizontal mirror. They together lower the system symmetry from $D_{6h}$
to $C_{3v}$, which supports the intrinsic NPHE.

Figure 2: (a) Low-energy bands for the honeycomb lattice model (11), with two
valleys at $K$ and $K^{\prime}$. The inset shows the honeycomb lattice, which
has two sublattices. (b) Calculated response coefficient
$\chi_{yxxx}^{\text{int}}$ versus Fermi energy $\mu$. (c-f) show the
distribution of (c) $\Lambda_{yxx}$, (d) $\Lambda_{xxx}$, (e) $G_{xy}$ and (f)
$\alpha_{yxxx}$ for the upper valence band in $k$-space. In the calculation,
we take $t=1$ eV, $t_{R}=10$ meV, and $\Delta=100$ meV.

Figure 2(a) shows the energy spectrum of the model, which exhibits a binary
valley structure at $K$ and $K^{\prime}$. The splitting of the conduction
(valence) band is due to the SOC. For the coordinate setup in Fig. 2(a), only
one tensor element $\chi_{yxxx}^{\text{int}}$ is relevant. In Figs. 2(c)-(f),
we plot the $k$-space distribution of relevant band geometric quantities,
including $\Lambda_{yxx}$, $\Lambda_{xxx}$, $G_{xy}$ and $\alpha_{yxxx}$, for
the upper valence band. One observes that these quantities are concentrated
around the two valleys, where the energy difference between the bands is
small. This can be understood by noting that similar to the Berry curvature,
all these geometric quantities encodes the interband coherence, so they are
peaked mainly at small-gap regions. For the model here, the intrinsic NPHE
should become large when the Fermi surface crosses the hot spots around the
band edge, which is confirmed by the result in Fig. 2(b).

This example shows that to have a sizable intrinsic NPHE, the system should
have its Fermi level near certain band near degeneracies, where the band
geometric quantities are pronounced. In addition, a sizable SOC is also a
necessary condition. Generally, without SOC, the band structure perturbed by
the $B$ field would consist of two decoupled spin channels. Each spin channel
can be regarded as a spinless system with an effective time reversal symmetry
[37, 38], which then suppresses the intrinsic NPHE. This point is also
verified in our model calculation.

Figure 3: (a) Top and (b) side views of MoSSe monolayer. (c) Brillouin zone of
MoSSe. (d) Calculated band structure along the high-symmetry paths. (e)
Enlarged view around the valence band degeneracy point at $\Gamma$, as
indicated by the arrow in (d). (f) Calculated intrinsic NPHE response
$\chi_{yxxx}^{\mathrm{int}}$, which shows a peak for Fermi level around the
point in (e).

Application to 2D MoSSe. As discussed, the advantage of an intrinsic effect is
that it can be evaluated from first principles to yield quantitative
predictions in concrete materials. Here, guided by the symmetry constraint and
lessons from the model study, we investigate the 2D Janus monolayer MoSSe,
which is a polar material synthesized in 2017 [39, 40] and under active
research since then [41, 42, 43].

The structure of monolayer MoSSe is shown in Figs. 3(a, b), which is similar
to H-MoS2 but with a whole layer of S replaced by Se. The crystal has space
group $P3m1$ and polar point group $C_{3v}$. It follows that the intrinsic
NPHE response obeys Eq. (9), with a single independent tensor element
$\chi_{yxxx}^{\text{int}}$. We evaluate $\chi_{yxxx}^{\text{int}}$ by
combining our theory with first-principles calculations (details are presented
in [30]). The calculated band structure in Fig. 3(d) matches well with
previous results [40]. The system is a semiconductor with a direct gap $\sim
1.47$ eV at $K$ and $K^{\prime}$. The degeneracy of the top two valence bands
is lifted by a sizable SOC with a large splitting $\sim 168$ meV at $K$. More
important for our discussion here is a Rashba-type SOC due to the the broken
horizontal mirror, a characteristic of the Janus structure. This manifests as
a Rashba-type splitting at the $\Gamma$ point [44], as indicated by the arrow
in Fig. 3(d) and more clearly in the enlarged Fig. 3(e). We find that this
band degeneracy point at $\Gamma$ gives a sizable contribution to the
intrinsic NPHE.

The calculated $\chi_{yxxx}^{\text{int}}$ as a function of chemical potential
$\mu$ for the $p$-doped case is shown in Fig. 3(f). Indeed, pronounced
response is found when $\mu$ is located around the Rashba-type band degeneracy
point at $\Gamma$. From the calculation, under an in-plane $B$ field of 1 T
the intrinsic nonlinear Hall conductivity can reach $\sim 3\times 10^{-4}$
$(\Omega\text{V})^{-1}$ if converted to the conventional 3D unit (by dividing
the layer thickness of 0.32 nm), when the Fermi level is about 80 meV below
the valence band maximum. This value is comparable to those calculated in
CuMnAs [45] and Mn2Au [46], and is two orders of magnitude larger than the
experimental result measured in CuMnAs [47], so it should be readily
detectable. In practice, the carrier doping can be conveniently controlled for
2D materials by electric gating [48, 49]. Therefore, our predicted NPHE and
its special angular dependence (9) (see Fig. 1(b)) can be directly tested in
experiment.

Discussion. We have proposed a new intrinsic transport effect along with a new
intrinsic material property, i.e., the intrinsic NPHE conductivity
$\chi^{\text{int}}$. It offers a promising characterization tool for a large
class of materials, especially the polar and chiral crystals, most of which do
not support the linear PHE nor the nonlinear anomalous Hall effect from Berry
curvature dipole. Actually, the monolayer MoSSe studied here represents such
an example.

Our subject here is on the intrinsic NPHE. As mentioned, there are also other
extrinsic contributions. In experiment, they can be separated according to
their different scaling with $\tau$ (e.g., by plotting against the
longitudinal conductivity) [31, 32]. In addition, we point out that the
extrinsic contribution is in general not purely a Hall response, namely, the
response tensor $\chi^{\text{ext}}_{abcd}$ is not guaranteed to be
antisymmetric in $a$ and $b$. It follows that their angular dependence will
generally differ from that of intrinsic NPHE.

## References

  * Jungwirth  _et al._ [2002] T. Jungwirth, Q. Niu, and A. H. MacDonald, Phys. Rev. Lett. 88, 207208 (2002).
  * Onoda and Nagaosa [2002] M. Onoda and N. Nagaosa, J. Phys. Soc. Jpn. 71, 19 (2002).
  * Nagaosa  _et al._ [2010] N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, and N. P. Ong, Rev. Mod. Phys. 82, 1539 (2010).
  * Murakami  _et al._ [2003] S. Murakami, N. Nagaosa, and S.-C. Zhang, Science 301, 1348 (2003).
  * Sinova  _et al._ [2004] J. Sinova, D. Culcer, Q. Niu, N. A. Sinitsyn, T. Jungwirth, and A. H. MacDonald, Phys. Rev. Lett. 92, 126603 (2004).
  * Sinova  _et al._ [2015] J. Sinova, S. O. Valenzuela, J. Wunderlich, C. H. Back, and T. Jungwirth, Rev. Mod. Phys. 87, 1213 (2015).
  * Xiao  _et al._ [2006] D. Xiao, Y. Yao, Z. Fang, and Q. Niu, Phys. Rev. Lett. 97, 026603 (2006).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * Tang  _et al._ [2003] H. X. Tang, R. K. Kawakami, D. D. Awschalom, and M. L. Roukes, Phys. Rev. Lett. 90, 107201 (2003).
  * Bowen  _et al._ [2005] M. Bowen, K.-J. Friedland, J. Herfort, H.-P. Schönherr, and K. H. Ploog, Phys. Rev. B 71, 172401 (2005).
  * Seemann  _et al._ [2011] K. M. Seemann, F. Freimuth, H. Zhang, S. Blügel, Y. Mokrousov, D. E. Bürgler, and C. M. Schneider, Phys. Rev. Lett. 107, 086603 (2011).
  * Nandy  _et al._ [2017] S. Nandy, G. Sharma, A. Taraphder, and S. Tewari, Phys. Rev. Lett. 119, 176804 (2017).
  * Kumar  _et al._ [2018] N. Kumar, S. N. Guin, C. Felser, and C. Shekhar, Phys. Rev. B 98, 041103(R) (2018).
  * Deng  _et al._ [2019] M.-X. Deng, H.-J. Duan, W. Luo, W. Y. Deng, R.-Q. Wang, and L. Sheng, Phys. Rev. B 99, 165146 (2019).
  * Ma  _et al._ [2019a] D. Ma, H. Jiang, H. Liu, and X. C. Xie, Phys. Rev. B 99, 115121 (2019a).
  * Breunig  _et al._ [2017] O. Breunig, Z. Wang, A. A. Taskin, J. Lux, A. Rosch, and Y. Ando, Nat. Commun. 8, 15545 (2017).
  * Wu  _et al._ [2018a] B. Wu, X.-C. Pan, W. Wu, F. Fei, B. Chen, Q. Liu, H. Bu, L. Cao, F. Song, and B. Wang, Appl. Phys. Lett. 113, 011902 (2018a).
  * Rakhmilevich  _et al._ [2018] D. Rakhmilevich, F. Wang, W. Zhao, M. H. W. Chan, J. S. Moodera, C. Liu, and C.-Z. Chang, Phys. Rev. B 98, 094404 (2018).
  * He  _et al._ [2019] P. He, S. S.-L. Zhang, D. Zhu, S. Shi, O. G. Heinonen, G. Vignale, and H. Yang, Phys. Rev. Lett. 123, 016801 (2019).
  * Esin  _et al._ [2021] V. D. Esin, A. V. Timonina, N. N. Kolesnikov, and E. V. Deviatov, J. Exp. Theor. Phys. 133, 792 (2021).
  * Zheng  _et al._ [2020] S.-H. Zheng, H.-J. Duan, J.-K. Wang, J.-Y. Li, M.-X. Deng, and R.-Q. Wang, Phys. Rev. B 101, 041408(R) (2020).
  * Rao  _et al._ [2021] W. Rao, Y.-L. Zhou, Y.-j. Wu, H.-J. Duan, M.-X. Deng, and R.-Q. Wang, Phys. Rev. B 103, 155415 (2021).
  * Battilomo  _et al._ [2021] R. Battilomo, N. Scopigno, and C. Ortix, Phys. Rev. Res. 3, L012006 (2021).
  * Gao  _et al._ [2014] Y. Gao, S. A. Yang, and Q. Niu, Phys. Rev. Lett. 112, 166601 (2014).
  * Liu  _et al._ [2022] H. Liu, J. Zhao, Y.-X. Huang, X. Feng, C. Xiao, W. Wu, S. Lai, W.-b. Gao, and S. A. Yang, Phys. Rev. B 105, 045118 (2022).
  * Sundaram and Niu [1999] G. Sundaram and Q. Niu, Phys. Rev. B 59, 14915 (1999).
  * Karplus and Luttinger [1954] R. Karplus and J. M. Luttinger, Phys. Rev. 95, 1154 (1954).
  * Chang and Niu [1995] M.-C. Chang and Q. Niu, Phys. Rev. Lett. 75, 1348 (1995).
  * Gao  _et al._ [2015] Y. Gao, S. A. Yang, and Q. Niu, Phys. Rev. B 91, 214405 (2015).
  * [30] See Supplemental Material for theoretical and calculation details.
  * Kang  _et al._ [2019] K. Kang, T. Li, E. Sohn, J. Shan, and K. F. Mak, Nat. Mater. 18, 324 (2019).
  * Lai  _et al._ [2021] S. Lai, H. Liu, Z. Zhang, J. Zhao, X. Feng, N. Wang, C. Tang, Y. Liu, K. S. Novoselov, S. A. Yang, and W.-b. Gao, Nat. Nanotechnol. 16, 869 (2021).
  * [33] Here, for $C_{2v}$, the twofold axis is normal to the plane.
  * Sodemann and Fu [2015] I. Sodemann and L. Fu, Phys. Rev. Lett. 115, 216806 (2015).
  * Liu  _et al._ [2011a] C.-C. Liu, H. Jiang, and Y. Yao, Phys. Rev. B 84, 195430 (2011a).
  * Liu  _et al._ [2011b] C.-C. Liu, W. Feng, and Y. Yao, Phys. Rev. Lett. 107, 076802 (2011b).
  * Vanderbilt [2018] D. Vanderbilt,  _Berry phases in electronic structure theory_ (Cambridge University Press, Cambridge, 2018).
  * Wu  _et al._ [2018b] W. Wu, Y. Liu, S. Li, C. Zhong, Z.-M. Yu, X.-L. Sheng, Y. X. Zhao, and S. A. Yang, Phys. Rev. B 97, 115125 (2018b).
  * Zhang  _et al._ [2017] J. Zhang, S. Jia, I. Kholmanov, L. Dong, D. Er, W. Chen, H. Guo, Z. Jin, V. B. Shenoy, L. Shi, _et al._ , ACS Nano 11, 8192 (2017).
  * Lu  _et al._ [2017] A.-Y. Lu, H. Zhu, J. Xiao, C.-P. Chuu, Y. Han, M.-H. Chiu, C.-C. Cheng, C.-W. Yang, K.-H. Wei, Y. Yang, _et al._ , Nat. Nanotechnol. 12, 744 (2017).
  * Dong  _et al._ [2017] L. Dong, J. Lou, and V. B. Shenoy, ACS Nano 11, 8242 (2017).
  * Guo  _et al._ [2021] Y. Guo, Y. Lin, K. Xie, B. Yuan, J. Zhu, P.-C. Shen, A.-Y. Lu, C. Su, E. Shi, K. Zhang, _et al._ , Proc. Natl. Acad. Sci. U.S.A. 118, e2106124118 (2021).
  * Zheng  _et al._ [2021] T. Zheng, Y.-C. Lin, Y. Yu, P. Valencia-Acuna, A. A. Puretzky, R. Torsi, C. Liu, I. N. Ivanov, G. Duscher, D. B. Geohegan, _et al._ , Nano Lett. 21, 931 (2021).
  * Yu  _et al._ [2021] S.-B. Yu, M. Zhou, D. Zhang, and K. Chang, Phys. Rev. B 104, 075435 (2021).
  * Wang  _et al._ [2021] C. Wang, Y. Gao, and D. Xiao, Phys. Rev. Lett. 127, 277201 (2021).
  * Liu  _et al._ [2021] H. Liu, J. Zhao, Y.-X. Huang, W. Wu, X.-L. Sheng, C. Xiao, and S. A. Yang, Phys. Rev. Lett. 127, 277202 (2021).
  * Godinho  _et al._ [2018] J. Godinho, H. Reichlová, D. Kriegner, V. Novák, K. Olejník, Z. Kašpar, Z. Šobáň, P. Wadley, R. P. Campion, R. M. Otxoa, _et al._ , Nature communications 9, 1 (2018).
  * Chen  _et al._ [2010] J. Chen, H. J. Qin, F. Yang, J. Liu, T. Guan, F. M. Qu, G. H. Zhang, J. R. Shi, X. C. Xie, C. L. Yang, K. H. Wu, Y. Q. Li, and L. Lu, Phys. Rev. Lett. 105, 176602 (2010).
  * Ma  _et al._ [2019b] Q. Ma, S.-Y. Xu, H. Shen, D. MacNeill, V. Fatemi, T.-R. Chang, A. M. Mier Valdivia, S. Wu, Z. Du, C.-H. Hsu, S. Fang, Q. D. Gibson, K. Watanabe, T. Taniguchi, R. J. Cava, E. Kaxiras, H.-Z. Lu, H. Lin, L. Fu, N. Gedik, and P. Jarillo-Herrero, Nature 565, 337 (2019b).