# Sequential Small Coronal Mass Ejections Observed In situ and in White-Light
Images by Parker Solar Probe

Brian E. Wood11affiliation: Naval Research Laboratory, Space Science Division,
Washington, DC 20375, USA<EMAIL_ADDRESS>, Phillip
Hess11affiliation: Naval Research Laboratory, Space Science Division,
Washington, DC 20375, USA<EMAIL_ADDRESS>, Yu Chen22affiliation:
Center for Space Plasma and Aeronomic Research (CSPAR), The University of
Alabama in Huntsville, Huntsville, AL 35805, USA , Qiang Hu22affiliation:
Center for Space Plasma and Aeronomic Research (CSPAR), The University of
Alabama in Huntsville, Huntsville, AL 35805, USA 33affiliation: Department of
Space Science, The University of Alabama in Huntsville, Huntsville, AL 35805,
USA

###### Abstract

We reconstruct the morphology and kinematics of a series of small transients
that erupt from the Sun on 2021 April 24 using observations primarily from
Parker Solar Probe (PSP). These sequential small coronal mass ejections (CMEs)
may be the product of continuous reconnection at a current sheet, a
macroscopic example of the more microscopic reconnection activity that has
been proposed to accelerate the solar wind more generally. These particular
CMEs are of interest because they are the first CMEs to hit PSP and be
simultaneously imaged by it, using the Wide-field Imager for Solar Probe
(WISPR) instrument. Based on imaging from WISPR and STEREO-A, we identify and
model six discrete transients, and determine that it is the second of them
(CME2) that first hits PSP, although PSP later more obliquely encounters the
third transient as well. Signatures of these encounters are seen in the PSP in
situ data. Within these data, we identify six candidate magnetic flux ropes
(MFRs), all but one of which are associated with the second transient. The
five CME2 MFRs have orientations roughly consistent with PSP encountering the
right sides of roughly E-W oriented MFRs, which are sloping back towards the
Sun.

Sun: coronal mass ejections (CMEs) — solar wind — interplanetary medium

##  1 Introduction

Parker Solar Probe (PSP) provides an opportunity to study coronal mass
ejections (CMEs) closer to the Sun than ever before. However, the fraction of
time PSP spends close to the Sun in its highly elliptical orbit is relatively
small, and the odds of being hit by a CME during a particular perihelion
passage are therefore relatively low. Thus, there are so far a limited number
of events available for study during the 15 orbits that PSP has made around
the Sun as of 2023 April. The most well-studied event dates back to the very
first perihelion passage, where on 2018 November 11 shortly after perihelion
(2018 November 6) PSP encountered a CME while at a distance of 53 R⊙ from the
Sun (Korreck et al., 2020; Nieves-Chinchilla et al., 2020; Good et al., 2020).
Other events observed in situ by PSP farther from the Sun have been analyzed
by Winslow et al. (2021), Hu et al. (2022), and Möstl et al. (2022).

The PSP mission is primarily designed to study the solar wind and CMEs using a
suite of on board particle and field detectors. However, PSP also possesses a
white light imaging instrument, the Wide-field Imager for Solar Probe (WISPR)
(Vourlidas et al., 2016). The two heliospheric imagers that consistitute
PSP/WISPR observe the solar wind in the ram direction of PSP’s orbit, and are
typically taking data only for about a month surrounding each perihelion. A
number of CMEs have been observed by WISPR during these periods, dating back
to the first orbit, with a small circular CME observed on 2018 November 1
(Hess et al., 2020; Rouillard et al., 2020), followed by a “streamer blob”
event seen very near perihelion on 2018 November 5 (Wood et al., 2020). The
transients studied so far have tended to be rather small and slow, but these
images have still allowed CME structures to be seen in much greater detail
than ever before (Liewer et al., 2020, 2021; Wood et al., 2021; Braga et al.,
2022; Howard et al., 2022).

Given that WISPR observes in the direction of the spacecraft motion, WISPR
sees transients passing in front of the spacecraft in its orbit, which are
typically not hitting PSP. We report here on the first CME that both directly
hits PSP and is imaged by WISPR while doing so, offering the first opportunity
for WISPR and PSP’s in situ instruments to be used to study the same CME.
Unfortunately, this CME encounter is a complicated one, with WISPR observing
many fronts passing through the field of view on 2021 April 24-25, just before
PSP’s eighth perihelion on April 29, with the spacecraft 46 R⊙ from the Sun.
Necessary clarity is provided by coronagraphic observations from STEREO-A,
showing that the eruptive activity is more accurately described as a series of
small CMEs erupting off the east limb, as viewed from STEREO-A’s perspective.
A central task of our analysis will be to utilize the WISPR and STEREO-A
images to determine which of the small transients are actually hitting PSP, a
necessary step toward properly interpreting the in situ signatures of the
activity observed by PSP. There is also a prominence eruption associated with
this activity, observed in EUV images from both STEREO-A and Solar Orbiter,
but we refer the reader to Niembro et al. (2023) for a detailed discussion of
the EUV data.

##  2 Observations

Figure 1: (a) On the left is an image of a series of three small CMEs
(numbered 1-3) from UT 06:53:30 on 2021 April 24, taken by the COR2-A
coronagraph on board STEREO-A. To the right of the image is a synthetic image
based on the 3-D morphological reconstruction of the CMEs described in Section
3. The real and synthetic images are displayed in a running-difference format,
and artificial noise is added to the synthetic image for aesthetic purposes.
(b) Another COR2-A image, from UT 22:53:30 on 2021 April 24, showing three
later small CMEs (numbered 4-6), with the synthetic model image again shown
next to it. (A movie version of this figure is available in the online version
of this article.)

Although the principal stimuli of our analysis are the images and in situ
observations from PSP on 2021 April 24-25, we first provide context for these
data by presenting coronagraphic observations at that time from 1 au, in
Figure 1. These images are from the COR2-A coronagraph on the STEREO-A
spacecraft, which observes at angular distances from Sun-center of
$0.7^{\circ}-4.2^{\circ}$ ($2.5-15.6$ R⊙) (Howard et al., 2008). The images
are displayed in a running-difference format, subtracting the previous image
to emphasize the dynamic CME fronts. On 2021 April 24, COR2-A observes a
series of small CMEs erupting off the east limb of the Sun. (The COR2-A
observations of the activity are covered more thoroughly in the movie version
of Figure 1.) This sequence of transients is not easy to separate into
distinct eruptions, but in Figure 1 we identify six separate little CMEs.
Ideally, it would be preferable if it could be determined if there are truly
separate magnetic flux rope (MFR) structures within these separate transients,
to ensure these are in fact separate eruptions, but as will be discussed in
more detail below, the available imaging data do not allow this to be clearly
determined. This kind of sustained activity is not uncommon in coronagraphic
data, and has been frequently observed by, for example, the Large Angle and
Spectrometric COronagraph (LASCO) instrument on board the SOlar and
Heliospheric Observatory (SOHO), which has been continously observing the Sun
since 1996 (Brueckner et al., 1995). These CMEs are small enough that they
could be considered in the category of “streamer blob” events (Sheeley et al.,
1997; Wang et al., 1998).

Figure 2: The positions of Earth, PSP, and STEREO-A in the ecliptic plane on
2021 April 24 (in HEE coordinates). At PSP’s position, the green and red lines
indicate the fields of view of the WISPR-I and WISPR-O detectors. At
STEREO-A’s position, the blue line indicates the field of view of COR2-A. The
dashed lines indicate the central trajectories of six small CMEs observed by
all these imagers.

Solar transients this small would normally not garner much attention, but in
this instance the eruptions are directed right at PSP, making this our first
opportunity to see what this kind of activity looks like up close. Figure 2
shows the locations of STEREO-A and PSP relative to Earth on 2021 April 24. It
also explicitly shows the fields of view of the COR2-A coronagraph, and of
PSP/WISPR’s two heliospheric imagers: WISPR-I, imaging at
$13^{\circ}-53^{\circ}$ from the Sun, and WISPR-O, imaging at
$50^{\circ}-108^{\circ}$ (Vourlidas et al., 2016). Dashed lines in the figure
indicate the central trajectories inferred for the six CMEs, based on the
stereoscopic analysis that will be described in Section 3. The CMEs all pass
very close to PSP, which is almost directly behind the Sun relative to Earth,
explaining why SOHO/LASCO operating near Earth is unable to see any of this
activity. As for STEREO-A imagery, the CMEs are observed by COR2-A, as shown
in Figure 1, but they are too faint to be discerned closer to the Sun by the
COR1-A coronagraph; and the heliospheric imagers, HI1-A and HI2-A, are
observing the wrong side of the Sun to view them. Thus, from STEREO-A we only
have the COR2-A data. As noted in Section 1, there is a prominence eruption
associated with this activity apparent in EUV imagery, but we refer the reader
to Niembro et al. (2023) for a detailed analysis of those data, along with a
discussion of potential signatures of the prominence in the PSP in situ data.
We do note here that we see no clear evidence of the prominence in the white
light images, which are naturally observing farther from the Sun than the EUV
images.

Figure 3: (a) On the left is a PSP/WISPR-I image of a series of fronts
(numbered 1-3) observed on UT 23:27:18 on 2021 April 24, associated with the
three CMEs observed by COR2-A in Figure 1(a). To the right of the image is a
synthetic image based on the 3-D morphological reconstruction of the CMEs
described in Section 3. The real and synthetic images are displayed in a
average-difference format, and artificial noise is added to the synthetic
image for aesthetic purposes. (b) Another WISPR-I image, from UT 10:57:15 on
2021 April 25, showing three later fronts (numbered 4-6), corresponding to the
COR2-A CMEs seen in Figure 1(b), with the synthetic model image again shown
next to it. (A movie version of this figure is available in the online version
of this article, which also includes imagery from WISPR-O as well as WISPR-I.)

In the WISPR images from 2021 April 24-25, a confusing series of fronts is
observed passing through the field of view, as illustrated by two WISPR-I
images in Figure 3. The movie version of this figure shows this more
completely, and shows the WISPR-O imagery as well. The images are shown in an
average-difference format, with an average WISPR-I image subtracted from the
sequence. This activity is clearly associated with the sequence of transients
seen by COR2-A in Figure 1, but connecting the WISPR-I fronts with the
individual small COR2-A CMEs is nontrivial. The numbers in Figure 3 indicate
what fronts we ultimately associate with each CME in Figure 1. We will be
referring to these CMEs using these numbers throughout the rest of this
article.

The first transient, CME1, is noticeably south of the ecliptic in Figure 1,
and is therefore identifiable in WISPR-I as a front only seen in the bottom
third of the images from about UT 14:27–22:27. The other transients, starting
with CME2, are all very close to the ecliptic plane. In the WISPR-I images,
CME2 is seen as a very broad, very faint front passing rapidly through the
field of view at UT 21:57–23:57, angled diagonally from the upper left to the
lower right side of the images. As will be shown explicitly in the 3-D
reconstruction described below in Section 3, this is what a CME front is
expected to look like when directly hitting the spacecraft, and it is CME2
that we infer to be the first to impact PSP. The third transient, CME3, is the
brightest in COR2-A, as shown most clearly in the movie version of Figure 1.
We connect it with the bright, semicircular WISPR-I front seen clearly in
Figure 3, which is the first of a series of similar looking fronts seen by
WISPR-I, the most distinct of which we connect with the CME4-CME6 sequence of
transients, which are particularly small and jet-like in COR2-A.

We have not described the WISPR-O images much (see movie version of Figure 3),
as these data are even harder to interpret than the WISPR-I images. The
detailed reconstruction of the eruptions described in the next section
provides some clarity. From this it is clear that CME1 is not seen by WISPR-O
at all, being too far south; CME2 is responsible for a bright, very broad
front passing through the field of view at the end of April 24, oriented
diagnonally through the field of view (as in WISPR-I); CME3 is a large semi-
circular front filling the field of view and passing rapidly through it from
about UT 4:34–6:34 on April 25; and the final three CME4–CME6 transients are
only partially seen along the top of the WISPR-O field of view starting at
about UT 11:34, and are hard to separate.

##  3 Three Dimensional Reconstruction of the Six CMEs

Modeling a CME’s appearance in images requires both a morphological analysis,
where the 3-D shape of the CME is approximated; and a kinematic analysis,
where the Sun-center distance and speed of the CME’s leading edge as a
function of time are determined. We first discuss the kinematic component. For
each CME, this begins with measurements of the elongation angle of the leading
edge from Sun-center, $\epsilon$, which then must be converted to actual
distances from Sun-center, $r$, using some geometric approximation. A
kinematic model can then be inferred from these $r$ values.

In past analyses, we have usually focused on measurements from one spacecraft,
with the choice of spacecraft based on which provides the best vantage point
and/or the most extensive set of measurements covering the widest range of
distances (e.g., Wood et al., 2017). However, a single spacecraft approach
does not work in this case, because we only have data from two sources,
STEREO/COR2-A and PSP/WISPR, which are observing the CMEs at different
distances from the Sun, with no overlap. Thus, we combine $\epsilon$
measurements from both COR2-A and WISPR, and infer $r$ from both. In practice,
we only consider WISPR-I from PSP, as for most of the CMEs WISPR-O is only
seeing parts of the CME fronts, with the leading edge not clearly in the field
of view. (The only clear exception would be CME3.) Even for WISPR-I,
meaningful measurements are not possible for CME1 or CME2. For CME1, the
leading edge is below the WISPR-I field of view, and with CME2 directly
hitting the spacecraft, PSP is not able to view its leading edge well.

Figure 4: (a) Logarithmic elongation angles ($\epsilon$, in degrees) are
plotted versus time for the leading edges of each of the six CMEs observed on
2021 April 24–25. The $t=0$ reference time is UT 01:23:30 on April 24.
Diamonds are STEREO COR2-A measurements, and squares are from PSP/WISPR-I. (b)
Sun-center distances, $r$, computed from $\epsilon$ assuming the Harmonic Mean
approximation. Note the severe discontinuity between the COR2-A and WISPR-I
measurements for CME3–CME6. (c) Sun-center distances computed assuming the
Fixed-$\phi$ approximation, with better COR2-A/WISPR-I agreement. Simple
kinematic models are fitted to the data (see text).

The elongation angle measurements are shown in Figure 4(a) for the six CMEs.
Converting $\epsilon$ to $r$ requires a geometric approximation. The simplest
one that is commonly used is the “Fixed-$\phi$” approximation, which assumes
the CMEs are narrow enough for their leading edges to be considered point-
like, in which case

$r=\frac{d\sin\epsilon}{\sin(\epsilon+\phi)},$ (1)

where $d$ is the distance from the observer to the Sun and $\phi$ is the angle
between the CME trajectory and the observer’s line of sight to the Sun (Kahler
& Webb, 2007; Sheeley et al., 2008). Alternatively, if the CME front is
approximated as a sphere centered halfway between the Sun and the leading
edge, one finds

$r=\frac{2d\sin\epsilon}{1+\sin(\epsilon+\phi)},$ (2)

which has been called the “Harmonic Mean” approximation (Lugaz et al., 2009).
More complicated relations exist (e.g., Davies et al., 2013), but we will only
consider the Fixed-$\phi$ and Harmonic Mean assumptions here. It should be
noted that the trajectory directions quantified by $\phi$ and shown explicitly
in Figure 2 are determined by the morphological part of the analysis described
below. In practice, the generation of a full CME reconstruction involves
iteration between the kinematic and morphological components of the analysis.

Figures 4(b) and 4(c) show the $r$ measurements that result when the Harmonic
Mean and Fixed-$\phi$ approximations are used, respectively. It is immediately
clear that the results are poor for the Harmonic Mean option. The distances
inferred from the WISPR measurements are too low relative to those inferred
from COR2-A. The situation is improved when Fixed-$\phi$ is assumed. Although
a more modest discontinuity still clearly exists, the implied kinematic
profile seems more reasonable, with the higher inferred WISPR distances.

In a survey of 28 Earth-directed CMEs well observed by STEREO, both Harmonic
Mean and Fixed-$\phi$ were considered (Wood et al., 2017), and in all but four
cases it was the Harmonic Mean approximation that was deemed better at
providing a reasonable kinematic profile that yielded correct Earth arrival
times and synthetic images that collectively best matched the multi-viewpoint
data. However, those CMEs are bigger events than the narrow transients being
modeled here. Because the 2021 April 24 CMEs are small and narrow, it is
ultimately not surprising that Fixed-$\phi$ does better in this particular
instance, given that the Fixed-$\phi$ approximation assumes a very narrow
transient morphology.

Following past analyses (e.g. Wood et al., 2017), we fit a simple multi-phase
kinematic model to the $r$ measurements in Figure 4(c), with each phase
corresponding to either a time of constant velocity or constant acceleration.
Full kinematic models such as this allow synthetic images of CMEs to be
computed for any time whatsoever after CME initiation, including times without
$r$ measurements (e.g., WISPR-O image times, or WISPR-I image times for CME1
and CME2). The resulting fits are represented as solid lines in the figure.
For CME1 and CME2, where we only have COR2-A measurements, we simply use a
single phase model assuming a constant velocity. It is fortunate that this
ultimately leads to synthetic images that reproduce the actual images
reasonably well (see below), which it should be mentioned would not be the
case for CME3–CME6, which require the WISPR $r$ constraints.

Table 1: CME Parameters  Parameter Description CME1 CME2 CME3 CME4 CME5 CME6
VCOR2A (km/s) COR2-A speed 312 415 323 363 376 366 VWISPR (km/s) WISPR-I speed
… … 233 222 225 274 $\lambda_{s}$ (deg) Trajectory longitude 175 171 181 182
182 182 $\beta_{s}$ (deg) Trajectory latitude -12 0 0 1 1 1 FWHMs (deg)
Angular width 29 25 25 15 15 15 $\alpha_{s}$ Leading edge shape 2.5 2.0 2.0
2.0 2.0 2.0

For CME3–CME6, we assume a three-phase model, with a period of constant
acceleration followed by a period of constant deceleration, and finally a
period of constant velocity. We have used this kind of model for many past
CMEs (e.g., Wood & Howard, 2009), in general representing the acceleration
phase of a fast CME, a gradual deceleration as the CME is slowed by
interaction with the slower ambient wind, and finally a constant velocity as
the terminal CME speed is reached. In this case, the meaning of the inferred
acceleration and deceleration phases is more ambiguous, as the fit will be
using the acceleration to try to account for the aforementioned discontinuity
between the COR2-A and WISPR-I measurements. For this reason, we do not bother
to show the velocity profiles that correspond to the fits in Figure 4(c).

However, in the first two rows of Table 1, we list mean velocities observed
for each of the CMEs in the COR2-A and WISPR-I fields of view, computed from
simple linear fits to the COR2-A and WISPR-I height-vs-time data points. It is
noteworthy that for CME3–CME6 we infer significantly lower velocities in
WISPR-I than for COR2-A. Although systematic uncertainties in this analysis
are high, we believe these decelerations are probably real, at least
qualitatively, as they can be explained as a natural consequence of these
later CMEs plowing into the slower trailing parts of CME1 and CME2.

Turning our attention to the morphological part of the analysis, we rely on
techniques used many times in the analysis of stereoscopic CME imaging, dating
back to Wood & Howard (2009), but more thoroughly described by Wood et al.
(2017). This is a forward modeling approach, involving the use of a
parametrized functional form for the assumed CME shape, with the parameters
adjusted to yield synthetic images that best match the observations. Synthetic
image calculation includes full computation of the Thomson scattering from the
model structure, although we only place mass on the surface of the structure
to outline its shape. We are not here attempting to reproduce image
brightnesses or measure CME masses, so the densities assumed at the surface
are arbitrary. A popular paradigm for CME structure is that of a magnetic flux
rope (MFR), i.e. a tube shape permeated with a helical magnetic field with
both legs streching back toward the Sun, so we have usually assumed an MFR
shape in 3-D CME reconstructions (Wood & Howard, 2009). However, for the
small, faint CMEs studied here, we really only perceive the transients as
narrow fronts, with nothing to suggest the presence or orientation of an MFR.
Thus, we instead assume here a simple lobular front shape, a prescription we
have generally used in the past to model CME shocks, when visible ahead of the
ejecta (e.g., Wood & Howard, 2009; Wood et al., 2017), but which has also been
used before to model the ejecta itself in other cases where an MFR shape is
not clearly apparent (Wood et al., 2009).

Figure 5: Three dimensional visualizations of the six reconstructed CME
fronts, for four different times. The fronts are displayed in an HEE
coordinate system, with the size of the Sun shown to scale in the figures. The
time-dependent location of PSP is shown as a small red sphere. The first two
frames show the positions of the first 3 CMEs (CME1–CME3) before and after
CME2 hits PSP, respectively. The final frame shows PSP entering the flank of
CME3. The last 3 very narrow CMEs (CME4-CME6) are buried inside CME3, as
labeled in the final frame.

Figure 5 shows the resulting model CME fronts, and their positions relative to
each other and to PSP. The parameters used to define the fronts are listed in
the last four rows of Table 1. The trajectory longitude and latitude of the
CMEs are $\lambda_{s}$ and $\beta_{s}$, respectively, based on a heliocentric
Earth ecliptic (HEE) coordinate system, with the positive x-axis pointed
toward Earth, and the positive z-axis pointed toward ecliptic north. The CME
longitudes are shown explicitly in Figure 2, which are within $11^{\circ}$ of
each other. Only CME1 has a latitude significantly out of the ecliptic
($\beta_{s}=-12^{\circ}$). The $\phi$ trajectory angles relative to STEREO-A
and PSP used in equations (1) and (2) above can be computed from
$\lambda_{s}$, $\beta_{s}$, and the known spacecraft positions. The
$\alpha_{s}$ parameter in Table 1 defines the shape of the leading edge, with
higher values of $\alpha_{s}$ being flatter and lower values being more
rounded (see, e.g., Wood et al., 2017).

The widths of the CME fronts are quantified in Table 1 as a full-width-at-
half-max, $FWHM_{s}$. The CMEs are all very narrow, with $FWHM_{s}$ ranging
from $15^{\circ}$ for CME4–CME6 to $29^{\circ}$ for CME1. These values can be
crudely compared with the widths inferred for CMEs modeled with MFRs, with
$FWHM_{s}$ in those cases defining the width of the CME in the plane of the
MFR. In the survey of 28 Earth-directed CMEs of Wood et al. (2017), the
narrowest CME had $FWHM_{s}=41.5^{\circ}$. Even the tiny streamer blob CME
observed by PSP/WISPR in 2018 November 5 had $FWHM_{s}=43.6^{\circ}$ (Wood et
al., 2020). All this emphasizes that we are not perceiving much lateral extent
for the small 2021 April 24 CMEs compared to past CMEs that we have studied
stereoscopically.

The CME trajectories and morphological parameters were determined by a trial
and error process of experimenting with different parameters to see which
parameters lead to synthetic images that best match the real data. The self-
similar expansion of the model CME fronts is defined by the kinematic models
in Figure 4(c). The synthetic images resulting from our best-fit models are
compared with the real images in Figures 1 and 3, with the movie versions of
these figures available in the online article providing a more thorough
comparison. Artificial noise is added to the synthetic images for aesthetic
purposes. The added noise level is the same in absolute terms in all the
images. The apparent lower noise in Figure 1(a) compared to 1(b) is illusory
due to the upper range of the color scale being extended upwards to match the
range of CME3, which is the brightest of the six CMEs.

It is difficult to estimate uncertainties in the parameters listed in Table 1.
There have been many attempts to estimate errors in CME reconstructions based
on multi-spacecraft imaging. A recent example is Verbeke et al. (2023), which
also summarizes much of the past work on the subject. A lot of this work
focuses on kinematic uncertainties, due to the interest in predicting 1 au
arrival times at Earth for CMEs (e.g., Kilpua et al., 2012; Mays et al., 2015;
Wood et al., 2017; Wold et al., 2018; Paouris et al., 2021), but morphological
uncertainties are studied as well (e.g., Thernisien et al., 2009; Mierla et
al., 2010; Jang et al., 2016; Balmaceda et al., 2018). Uncertainties in an
analysis of this nature will depend greatly on the characteristics of the
particular event in question, and on the particular viewing geometry for that
event. The applicability of these studies to the very unique viewing geometry
of the 2021 April 24 events directed right at PSP is unclear.

One characteristic of our reconstruction that is worthy of note is that the
synthetic WISPR images are remarkably sensitive to the assumed trajectory
direction, with shifts of only a degree or two making noticeable changes to
the images. This is simply due to the extreme close proximity of the CMEs to
the spacecraft (see Figure 2). Perhaps the best examples of this are the
latitude parameters for CME3–CME6. For CME3, we infer $\beta_{s}=0^{\circ}$,
while for CME4–CME6 we infer $\beta_{s}=1^{\circ}$. Such a small difference
would usually not be significant in a CME reconstruction of this nature, but
in this case it is significant. This is best indicated by the WISPR-O data,
where the CME3 front is fully visible in the WISPR-O field of view, but the
CME4–CME6 fronts are at least halfway above the field of view. Only a single
degree of latitude shift in the reconstruction is enough to reproduce this.

With the image-based CME reconstruction now complete, we can finally address
the central question of which CMEs are responsible for the in situ signatures
observed by PSP. Our conclusions are illustrated by the four time frames shown
in Figure 5. The first CME, CME1, is directed south of the ecliptic
($\beta_{s}=-12^{\circ}$) and therefore does not hit PSP. It is CME2 that
makes a direct hit to PSP at the end of April 24. The final panel of Figure 5
then shows PSP encountering the flank of CME3, with PSP’s orbital motion being
partly responsible for carrying the spacecraft into CME3, in addition to the
expansion of the CME itself. The last three CMEs, CME4-6, are highly blended
by the time they reach the vicinity of PSP. It is likely that PSP encounters
trailing parts of these little transients near the beginning of April 26.

##  4 Comparison with PSP In-Situ Data

Figure 6: Plasma and field measurements from PSP/SWEAP and PSP/FIELDS during
the passage of the 2021 April 24 transients over the spacecraft. The top two
panels are proton number density and velocity. The third panel is the magnetic
field at PSP, both the total field and the $B_{x}$, $B_{y}$, and $B_{z}$
coordinates in an RTN coordinate system. The final panel is the proton
temperature. Vertical dashed lines indicate the arrival times of CME2 and CME3
predicted by the 3-D reconstruction in Figure 5, and the density and velocity
profiles predicted by this reconstruction are explicitly shown as well. The
velocity profiles match the data well. The 3-D reconstruction implies that the
broad density increase early in April 26 may be due to the last 3 modeled
transients, CME4-6.

We will now compare the predictions of our image-based event reconstruction
with measurements from PSP’s plasma and field instruments, in order to see if
the in situ data provide support for the reconstruction. The in situ data are
shown in Figure 6. The proton density, velocity, and temperature measurements
are from the Solar Wind Electrons Alphas & Protons (SWEAP) instrument on PSP
(Kasper et al., 2016), while the magnetic field measurements are from the
FIELDS instrument (Bale et al., 2016).

Vertical dashed lines in Figure 6 indicate when the 3-D reconstruction has PSP
being hit by CME2 and CME3, the first right at the end of April 24, and the
second midway through April 25. The PSP in situ data broadly support these
predictions. The predicted CME2 encounter time corresponds to a substantial
increase in velocity, magnetic field, and temperature. Proton densities
actually show a broad peak ahead of this arrival time, suggesting a pile-up
region ahead of the CME. Possible signatures of the predicted CME3 encounter
are more subtle, but there is a significant density increase at that time,
along with a velocity decrease, and a modest temperature peak.

We do not try to show predicted encounter times for CME4-6, because these last
three little transients are rather blended and are predicted to have an
extended grazing incidence encounter with PSP, making predictions more
uncertain. However, we do note in Figure 6 the existence of a broad density
peak during the first half of April 26, which would correspond to the time
when we would expect PSP to be affected by the passage of CME4-6.

The density and velocity profiles predicted at PSP for CME2 and CME3 by our
3-D reconstruction are explicitly shown in the top two panels of Figure 6. The
predicted velocity profiles indicate decreasing velocity behind the CME
arrival, characteristic of self-similar expansion. These predictions agree
well with the PSP/SWEAP observations, particularly for CME2. This is
encouraging, considering the difficulties in deriving accurate kinematic
models for these events (see Section 3).

A precise replication of the observed density profile is not expected,
considering the very rudimentary density distribution assumed in the
reconstruction, with mass placed only on the surfaces of the CME fronts.
Gaussian density profiles across the surface are assumed, so the predicted PSP
density profiles are roughly Gaussian as well. For CME2 the predicted profile
is very narrow, indicative of the direct hit of CME2 on PSP, carrying PSP
quickly through the CME surface. In contrast, the predicted density profile
for CME3 is broad, as PSP enters CME3 with more of a grazing incidence (see
last panel of Figure 5).

As for the magnetic field behavior, before the transients arrive PSP sees a
positive-polarity Parker spiral behavior, with positive $B_{x}$ and much
weaker negative $B_{y}$. This changes in the pile-up region ahead of CME2,
with $B_{x}$ and $B_{y}$ reversing sign, and the overall field strength
increasing dramatically as CME2 arrives at the end of April 24. However, after
CME2’s arrival, the field behavior remains relatively consistent, despite the
passages of the later transients. Throughout this period, $B_{x}$ remains
consistently negative and $B_{y}$ generally positive.

One obvious thing to do with the in situ data is to look for evidence of MFRs,
as it is generally assumed that MFRs lie at the hearts of all CMEs (Vourlidas
et al., 2013). We have already noted that for these small, relatively faint
CMEs, we perceive the transients in the white light images as simple lobular
fronts, with no clear evidence of an MFR shape, but this does not preclude the
existence of undetectable MFRs behind the visible leading fronts. However, the
in situ data do not show obvious “magnetic cloud” signatures, which are
considered the clearest indication of an MFR; with smoothly rotating, strong
fields, low density, and low temperature (Burlaga et al., 1981; Marubashi,
1986; Burlaga, 1988; Lepping et al., 1990, 2015; Bothmer & Schwenn, 1998). The
absence of such a signature is particularly notable for CME2, for which our
3-D reconstruction implies a direct hit on PSP, which we would expect to yield
an MFR signature in the in situ data, if an MFR is in fact present.

Figure 7: (a) The magnetic field panel from Figure 6 is reproduced, but with
shaded regions indicating six intervals identified as possible MFRs by an
automated routine (see Table 1). Vertical dashed lines indicate the arrival
times of CME2 and CME3, as in Figure 6. (b) Grad-Shafranov reconstructed cross
sections of the six MFRs. The color scale is the out-of-plane axial field in
nT. The contours are the in-plane transverse magnetic field lines. The white
and green arrows along $y=0$ are the measured transverse field and remaining
flow vectors, respectively, along the spacecraft path.

Nevertheless, there is field variability during the CME encounter, which could
be indicative of distinct magnetic features that could be interpreted as MFRs.
In order to explore this, we use an automated routine for finding MFRs within
the solar wind, developed by Hu et al. (2018). Most MFRs found by this
technique are quite small, with a median encounter duration of 20 minutes in
Wind data. The basis of the routine is the Grad-Shafranov (GS) technique for
inferring MFR orientations and deriving two-dimensional configurations for
clear magnetic cloud events within in situ solar wind data (Hu & Sonnerup,
2002; Hu, 2017). Unlike more classic techniques for modeling such structures
(e.g., Lepping et al., 1990, 2015; Nieves-Chinchilla et al., 2016), the GS
methodology does not assume that the MFR has a symmetric circular or
elliptical MFR cross section. The broader, more general theoretical
requirements of the GS approach allows the automated Hu et al. (2018) MFR
search routine to identify a greater number of potential MFR candidates within
the solar wind than more restrictive techniques.

Table 2: Flux Rope Properties  MFR# Start Time Duration Orientation $F_{\phi}$
$F_{\psi}$ $F_{\psi}/F_{\phi}$ (s) $\phi$ (deg) $\theta$ (deg) ($10^{10}$ Wb)
($10^{10}$ Wb) 1 2021-04-24 19:22:08 1904 160 80 0.019 1.87 96.4 2 2021-04-25
01:18:12 4844 140 100 2.21 36.7 16.6 3 2021-04-25 02:56:40 2520 120 100 2.09
39.9 19.1 4 2021-04-25 04:08:32 1372 120 100 0.49 11.0 22.2 5 2021-04-25
05:08:16 21644 140 80 61.0 194 3.18 6 2021-04-25 15:08:24 2968 140 130 1.19
37.9 31.8

Hu et al. (2018) identified 74,241 small MFRs in Wind data from 1996–2016.
Chen & Yu (2022) applied these techniques to PSP data from the first six PSP
orbits, finding nearly 6000 MFRs. We now similarly apply the algorithm to our
time interval of interest in 2021 April 24-26. Focusing only on longer
candidates ($\geq 20$ min), we identify six possible MFR structures for which
we can compute a full GS reconstruction. The time intervals of these MFRs are
listed in Table 2 and displayed in Figure 7(a). The first MFR (MFR1) is in the
pile-up region ahead of CME2’s arrival at PSP based on the 3-D reconstruction.
The next four MFRs are sequential within the time period that we associate
with the CME2 encounter. These four MFRs collectively account for most of that
time period. The last of these (MFR5) is by far the largest, with an encounter
duration of 6.0 hr. The final MFR (MFR6) is in the CME3 encounter period.

The full GS reconstructions of the six MFRs are shown in Figure 7(b), with a
color scale that maps the axial field component of the MFR in the cross-
sectional plane. In the GS modeling, the MFR is assumed to extend infinitely
out of the plane in both directions, and the figures indicate the inferred
spacecraft path through these structures along $y=0$, projected into the
cross-sectional plane. The orientation of the MFR is indicated by the $\phi$
and $\theta$ directions in Table 2, in the spacecraft-based RTN coordinate
system. The $\phi$ angle indicates the azimuthal direction, with
$\phi=0^{\circ}$ and $\phi=180^{\circ}$ pointed radially away and toward the
Sun, respectively; while $\theta$ is the polar angle, with $\theta=0^{\circ}$
and $\theta=180^{\circ}$ pointed toward ecliptic north and south,
respectively.

Figure 8: (a) The 3-D magnetic field configuration of MFR5 (see Figure 7),
viewed from a positive-R direction, looking toward the Sun, with the green dot
indicating PSP’s path through the structure. The field lines are traced from
the cross-sectional map on the left. (b) Here MFR5 is viewed from a positive-N
direction, looking down on the MFR with the Sun toward the bottom, and the
green arrow indicating PSP’s path through the MFR.

With $\theta=80^{\circ}-100^{\circ}$, the five CME2-associated MFRs are all
oriented close to the ecliptic plane. With $\phi=120^{\circ}-160^{\circ}$,
they are all encountered with the right side of the MFR slanted toward the Sun
relative to the spacecraft-Sun line. Figure 8 explicitly shows the orientation
of MFR5, with panel (b) showing how the MFR is tilted relative to the
spacecraft’s path in the ecliptic plane. With actual MFR structures assumed to
have legs that stretch back toward the Sun, the inferred [$\phi$,$\theta$]
angles could be consistent with PSP encountering the right sides of roughly
E-W oriented MFRs. However, it should be noted that the GS analysis assumes
that the MFR extends infinitely in directions out of the cross-sectional
plane. This introduces significant uncertainty in the interpretation of the
MFR orientation angles on large scales. In the example of the large CME2 MFR
(e.g., MFR5), significant curvature would be expected for it to be contained
within the confines of the rather narrow CME outline defined by the CME
lobular front inferred from the 3-D reconstruction.

Besides the MFR orientation information, Table 2 also provides estimates of
the toroidal and poloidal magnetic fluxes of the GS-constructed MFRs,
$F_{\phi}$ and $F_{\psi}$, respectively. The toroidal or axial flux is easiest
to understand, as it is simply the magnetic flux through the cross-sectional
plane illustrated in Figure 7 for each MFR. The poloidal flux is the
integration of the azimuthal field around the central axis times the area of
the surface through which this field is passing. This flux is much more
uncertain than $F_{\phi}$, since one must first assume that the azimuthal
field times the MFR radius is unchanged along the MFR length (i.e., under a
2-D geometry), and also a total length of the MFR must be known or assumed. In
principle, this could come from an image-based reconstruction of the MFR
shape, but in this case the MFR shapes are not discerned in the images and we
can only infer an outline of the CME leading edge. Still, we can use the
length of the curve outlining the lobular fronts as our estimate for MFR
length. At a distance from the Sun of 45 R⊙, this length is 94.6 R⊙ based on
the CME2 front shown in Figure 5, so the $F_{\psi}$ values in Table 2 are
computed assuming this length. Finally, the last column of Table 2 lists the
ratio of $F_{\psi}$ to $F_{\phi}$. These ratios are higher for the smaller
MFRs.

##  5 Summary and Discussion

The 2021 April 24-25 event is the first CME that directly hits PSP and is
simultaneously imaged by it. We have performed a full 3-D reconstruction of
the event based on the available imagery, which unfortunately is limited to
COR2-A and WISPR. The reconstruction is also complicated because the activity
actually consists of multiple little CMEs erupting sequentially from the solar
corona, but we have successfully connected six distinct, small CMEs seen by
COR2-A near the Sun with fronts seen by WISPR-I and WISPR-O as the transients
pass either over or very near PSP. The analysis clearly identifies CME2 as the
transient that first impacts PSP directly, with PSP later predicted to
encounter CME3 more obliquely.

We have found support for the 3-D reconstruction from PSP in situ data, with
density, temperature, and field signatures of the expected CME2 and CME3
impacts present in the in situ data. The declining velocity profiles observed
in situ are also in excellent agreement with those predicted by the
reconstruction, which simply assumes self-similar expansion. A broad density
enhancement early in April 26 may be due to the collective influence of
CME4-6.

Finally, we have used the GS-based automated MFR finding routine of Hu et al.
(2018) to identify MFR candidates within the PSP encounter time, ultimately
focusing on six of the longer candidates, all but the last associated with the
CME2 encounter period. The fifth MFR is by far the largest, with an encounter
duration of about 6 hours. The five CME2 MFRs have orientations roughly
consistent with PSP encountering the right sides of roughly E-W oriented MFRs,
which are sloping back towards the Sun.

The multiple eruptions seen by COR2-A are suggestive of intermittent
reconnection occuring near the heliospheric current sheet, an interpretation
consistent with the PSP in situ data showing a clear change in field polarity
at the beginning of the CME encounter. Small-scale magnetic islands are
commonly observed near current sheets (Khabarova et al., 2015), which could be
interpreted as small MFR structures, and which are potential sources of
particle acceleration (Zank et al., 2017; Van Eck et al., 2022). After initial
formation, small MFRs are expected to merge into larger ones (Drake et al.,
2021). It is possible that the sequence of CME2 MFRs in Figure 7 is indicative
of such a process, with a large, dominant MFR (MFR5) preceded by much smaller
MFRs that never merged with the large MFR before flowing away from the
reconnection site.

It is becoming more popular to envision microscale interchange reconnection as
being responsible for both the ubiquitous switchback structures observed by
PSP (Drake et al., 2021), and potentially for the solar wind as a whole (Wang,
2020; Raouafi et al., 2023). Activity such as the 2021 April 24-25 events may
provide a way to study this important process on a more macroscopic scale,
large enough for the intermittent outflows from such reconnection to be
apparent in imaging. The 2021 April 24-25 events are potentially particularly
useful for such purposes, with PSP providing both imaging and in situ
observations close to the Sun. Thus, these data may be worthy of more detailed
analysis in the future.

Financial support was provided by the Office of Naval Research. QH and YC
acknowledge support from NASA through grants 80NSSC21K1763 and 80NSSC21K0003,
and from NSF through award AGS-2229065. Parker Solar Probe was designed,
built, and is now operated by the Johns Hopkins Applied Physics Laboratory as
part of NASA’s Living with a Star (LWS) program (contract NNN06AA01C). We
particularly acknowledge the WISPR instrument team, funded by NASA through
grant NNG11EK11I, the SWEAP team led by J. Kasper, and the FIELDS experiment
led by S. Bale. The STEREO/SECCHI data are produced by a consortium of NRL
(US), LMSAL (US), NASA/GSFC (US), RAL (UK), UBHAM (UK), MPS (Germany), CSL
(Belgium), IOTA (France), and IAS (France). In addition to funding by NASA,
NRL also received support from the USAF Space Test Program and ONR.

## References

  * Bale et al. (2016) Bale, S. D., Goetz, K., Harvey, P. R., et al. 2016, Space Sci. Rev., 204, 49
  * Balmaceda et al. (2018) Balmaceda, L. A., Vourlidas, A., Stenborg, G., & Del Lago, A. 2018, ApJ, 863, 57
  * Bothmer & Schwenn (1998) Bothmer, V., & Schwenn, R. 1998, Ann. Geophys., 16, 1
  * Braga et al. (2022) Braga, C. R., Vourlidas, A., Liewer, P. C., et al. 2022, ApJ, 938, 13
  * Brueckner et al. (1995) Brueckner, G. E., Howard, R. A., Koomen, M. J., et al. 1995, Sol. Phys., 162, 357
  * Burlaga (1988) Burlaga, L. F. 1988, JGR, 93, 7217
  * Burlaga et al. (1981) Burlaga, L., Sittler, E., Mariani, F., & Schwenn, R. 1981, JGR, 86, 6673
  * Chen & Yu (2022) Chen, Y., & Hu, Q. 2022, ApJ, 924, 43
  * Davies et al. (2013) Davies, J. A., Perry, C. H., Trines, R. M. G. M., et al. 2013, ApJ, 777, 167
  * Drake et al. (2021) Drake, J. F., Agapitov, O., Swisdak, M., et al. 2021, A&A, 650, A2
  * Good et al. (2020) Good, S. W., Kilpua, E. K. J., Ala-Lahti, M., et al. 2020, ApJ, 900, L32
  * Hess et al. (2020) Hess, P., Rouillard, A. P., Kouloumvakos, A., et al. 2020, ApJS, 246, 25
  * Howard et al. (2008) Howard, R. A., Moses, J. D., Vourlidas, A., et al. 2008, Space Sci. Rev., 136, 67
  * Howard et al. (2022) Howard, R. A., Stenborg, G., Vourlidas, A., et al. 2022, ApJ, 936, 43
  * Hu (2017) Hu, Q. 2017, Sci. China Earth Sci., 60, 1466
  * Hu et al. (2022) Hu, Q, He, W., & Chen, Y. 2022, Frontiers in Physics, 10, 960315
  * Hu & Sonnerup (2002) Hu, Q., & Sonnerup, B. U. Ö. 2002, JGR, 104, 6899
  * Hu et al. (2018) Hu, Q., Zhang, J., Chen, Y., le Roux, J., & Zhao, L. 2018, ApJS, 239, 12
  * Jang et al. (2016) Jang, S., Moon, Y. J., Kim, R. S., Lee, H., & Cho, K. S. 2016, ApJ, 821, 95
  * Kahler & Webb (2007) Kahler, S. W., & Webb, D. F. 2007, JGR, 112, A09103
  * Kasper et al. (2016) Kasper, J. C., Abiad, R., Austin, G., et al. 2016, Space Sci. Rev., 204, 131
  * Khabarova et al. (2015) Khabarova, O., Zank, G. P., Li, G., et al. 2015, ApJ, 808, 181
  * Kilpua et al. (2012) Kilpua, E. K. J., Mierla, M., Rodriguez, L., et al. 2012, Sol. Phys., 279, 477
  * Korreck et al. (2020) Korreck, K. E., Szabo, A., Nieves-Chinchilla, T., et al. 2020, ApJS, 246, 69
  * Lepping et al. (1990) Lepping, R. P., Jones, J. A., & Burlaga, L. F. 1990, JGR, 95, 11957
  * Lepping et al. (2015) Lepping, R. P., Wu, C. -C., Berdichevsky, D. B., & Szabo, A. 2015, Sol. Phys., 290, 2265
  * Liewer et al. (2020) Liewer, P. C., Qiu, J., Penteado, P., et al. 2020, Sol. Phys., 295, 140
  * Liewer et al. (2021) Liewer, P. C., Qiu, J., Vourlidas, A., Hall, J. R., & Penteado, P. 2021, A&A, 650, A32
  * Lugaz et al. (2009) Lugaz, N., Vourlidas, A., & Roussev, I. I. 2009, Ann. Geophys., 27, 3479
  * Marubashi (1986) Marubashi, K. 1986, Adv. Space Res., 6, 335
  * Mays et al. (2015) Mays, M. L., Taktakishvili, A., Pulkinnen, A., et al. 2015, Sol. Phys., 290, 1775
  * Mierla et al. (2010) Mierla, M., Inhester, B., Antunes, A., et al. 2010, Ann. Geophys., 28, 203
  * Möstl et al. (2022) Möstl, C., Weiss, A. J., Reiss, M. A., et al. 2022, ApJ, 924, L6
  * Niembro et al. (2023) Niembro, T., Seaton, D. B., Hess, P., et al. 2023, FrASS, in press
  * Nieves-Chinchilla et al. (2016) Nieves-Chinchilla, T., Linton, M. G., Hidalgo, M. A., et al. 2016, ApJ, 823, 27
  * Nieves-Chinchilla et al. (2020) Nieves-Chinchilla, T. N., Szabo, A., Korreck, K. E., et al. 2020, ApJS, 246, 63
  * Paouris et al. (2021) Paouris, E., Vourlidas, A., Papaioannou, A., & Anastasiadis, A. 2021, Space Weather, 19, e02617
  * Raouafi et al. (2023) Raouafi, N. E., Stenborg, G., Seaton, D. B., et al. 2023, ApJ, 945, 28
  * Rouillard et al. (2020) Rouillard, A. P., Poirier, N., Lavarra, M., et al. 2020, ApJS, 246, 72
  * Sheeley et al. (2008) Sheeley, N. R., Jr., Herbst, A. D., Palatchi, C. A., et al. 2008, ApJ, 675, 853
  * Sheeley et al. (1997) Sheeley, N. R., Jr., Wang, Y. -M., Hawley, S. H., et al. 1997, ApJ, 484, 472
  * Thernisien et al. (2009) Thernisien, A., Vourlidas, A., & Howard, R. A. 2009, Sol. Phys., 256, 111
  * Van Eck et al. (2022) Van Eck, K., le Roux, J., Chen, Y., Zhao, L. L., & Thompson, N. 2022, ApJ, 933, 80
  * Verbeke et al. (2023) Verbeke, C., Leila Mays, M., Kay, C., et al. 2023, Adv. Space Res., in press (https://doi.org/10.1016/j.asr.2022.08.056)
  * Vourlidas et al. (2016) Vourlidas, A., Howard, R. A., Plunkett, S. P., et al. 2016, Space Sci. Rev., 204, 83
  * Vourlidas et al. (2013) Vourlidas, A., Lynch, B. J., Howard, R. A., & Li, Y. 2013, Sol. Phys., 284, 179
  * Wang (2020) Wang, Y. -M. 2020, ApJ, 904, 199
  * Wang et al. (1998) Wang, Y. -M., Sheeley, N. R., Jr., Walters, J. H., et al. 1998, ApJ, 498, L165
  * Winslow et al. (2021) Winslow, R. M., Lugaz, N., Scolini, C., & Galvin, A. B. 2021, ApJ, 916, 94
  * Wold et al. (2018) Wold, A. M., Mays, M. L., Taktakishvili, A., et al. 2018, JSWSC, 8, A17
  * Wood et al. (2020) Wood, B. E., Hess, P., Howard, R. A., et al. 2020, ApJS, 246, 28
  * Wood et al. (2021) Wood, B. E., Braga, C. R., & Vourlidas, A. 2021, ApJ, 922, 234
  * Wood & Howard (2009) Wood, B. E., & Howard, R. A. 2009, ApJ, 702, 901
  * Wood et al. (2009) Wood, B. E., Howard, R. A., Thernisien, A., Plunkett, S. P., & Socker, D. G. 2009, Sol. Phys., 259, 163
  * Wood et al. (2017) Wood, B. E., Wu, C. -C., Lepping, R. P., et al. 2017, ApJS, 229, 29
  * Zank et al. (2017) Zank, G. P., Adhikari, L., Hunana, P., et al. 2017, ApJ, 835, 147


# Impact of chemical disorder on magnetic exchange interactions in L10-FeNi
(tetrataenite)

Ankit Izardar  Claude Ederer  Materials Theory, ETH Zürich, Wolfgang-Pauli-
Strasse 27, 8093 Zürich, Switzerland

###### Abstract

We investigate the effect of chemical disorder on the magnetic exchange
couplings and the Curie temperature ($T_{\text{c}}$) in L10-ordered FeNi using
first-principles-based calculations. We use supercells to model chemical
disorder, to account for the specific symmetry-broken local chemical
environments around the individual atoms. We find a very strong variation of
the most dominant first-nearest neighbor Fe-Fe interaction for different
inequivalent Fe-Fe pairs, ranging from around 5 meV to 37 meV, compared to a
coupling strength of 27 meV in the ordered state. To estimate the influence of
such strong variations of the magnetic coupling constants on the Curie
temperature of the disordered or partially ordered state, we study a simple
Heisenberg model with random Gaussian-distributed nearest neighbor couplings
on an fcc lattice. Our Monte Carlo simulations for this model indicate that
strongly varying exchange couplings, such as those obtained for FeNi, can lead
to a reduction of $T_{\text{c}}$ of around 10 % relative to the one obtained
using only the average coupling.

##  I Introduction

The chemically ordered ferromagnet L10-FeNi (tetrataenite) has recently
generated considerable interest as a rare-earth-free, low-cost permanent
magnet, due to its high magneto-crystalline anisotropy energy and large
saturation magnetization Paulevé  _et al._ (1968); Kojima  _et al._ (2014);
Lewis  _et al._ (2014a, b); Edström  _et al._ (2014); Werwiński and Marciniak
(2017); Cui  _et al._ (2018); Tian  _et al._ (2019, 2020); Izardar and Ederer
(2020); Tian  _et al._ (2021). Since its discovery by Néel and coworkers in
the early 1960s Pauleve  _et al._ (1962); Néel  _et al._ (1964), several
attempts have been made to synthesize L10-FeNi with a high degree of chemical
order Goto  _et al._ (2017); Makino  _et al._ (2015); Saito  _et al._ (2019);
Kojima  _et al._ (2014). However, the synthesis of a fully ordered structure
remains challenging, due to the rather low order-disorder transition
temperature, which is around 590 K Pauleve  _et al._ (1962); Néel  _et al._
(1964); Reuter  _et al._ (1989). At this temperature, the diffusivity of atoms
is too low for the ordered structure to form on reasonable timescales.
Therefore, “naturally occurring” tetrataenite has only been found in iron
meteorites Albertsen  _et al._ (1978); Petersen  _et al._ (1977); Danon  _et
al._ (1979, 1980); Clarke and Scott (1980). Due to the difficulties in
obtaining fully ordered samples, it becomes essential to investigate and
understand how deviations from the perfect order affect the magnetic
properties, in particular the Curie temperature, $T_{\text{c}}$, and magnetic
anisotropy, of L10-FeNi.

Experimentally, only the Curie temperature of the disordered system
($T_{c}\approx 785$ K Onodera  _et al._ (1981); Wei  _et al._ (2014)) is
accessible, since the ordered system disorders on heating at temperatures
above $\sim$ 700 K. Thereby, the effective “disordering temperature” depends
strongly on the heating rate Santos  _et al._ (2015). On time-scales typical
for magnetization measurements, disordering occurs around 820 K, i.e., above
the Curie temperature of the disordered system, and results in an abrupt
vanishing of the magnetization Lewis  _et al._ (2014a). This indicates that
the nominal Curie temperature of the ordered system would be noticeably higher
than that of the disordered system.

Several previous studies have used first-principles calculations to obtain the
Curie temperature in L10-FeNi. For example, Edström et al. Edström  _et al._
(2014) and Tian et al. Tian  _et al._ (2020) obtained values for the Curie
temperature of the ordered phase of 916 K and 780 K, respectively, using
slightly different electronic structure methods. They also found that chemical
disorder leads to a reduction of the Curie temperature.

Both of these studies have used the coherent-potential approximation (CPA)
Soven (1967); Gyorffy (1972) to incorporate chemical disorder in the material.
The CPA is based on an effective medium description of the atomic
environments, and thus provides a very efficient method for the treatment of
disorder effects in random alloys using only a single unit cell. However, CPA
does not include effects related, e.g,. to the local symmetry-breaking of a
specific chemical environment around an individual atom. Such effects beyond
CPA can be particularly relevant, e.g., for the magnetic anisotropy, as we
showed in our previous work, where we have used supercells with different
distributions of Fe and Ni atoms to investigate the effect of chemical
disorder Izardar and Ederer (2020); Si  _et al._ (2021). These calculations
indicate that, for example, a moderate increase in Fe content, while reducing
the degree of chemical order in the system, leads to an increase of the
magneto-crystalline anisotropy energy, an effect not captured within the CPA.

In the present work, we use first principles calculations based on density
functional theory (DFT) to provide further insights into the effect of
variations in the local chemical environment on the magnetic exchange
interactions, and consequently the Curie temperature, in partially ordered
FeNi. To model the chemical disorder, we follow a similar approach as in our
previous work Izardar and Ederer (2020); Si  _et al._ (2021), i.e., we
incorporate the effect of a disordered local atomic environment on the
magnetic coupling by employing supercells. We find that treating the disorder
on a local level gives rise to remarkably strong variations of the first
nearest neighbour Fe-Fe coupling. Our subsequent analysis of the correlation
between the specific local chemical environment and the corresponding magnetic
exchange couplings indicates that the magnetic coupling is governed by long
range effects that clearly go beyond the closest neighbor environment. This
makes it extremely challenging to consider such configuration-dependent
couplings for the calculation of $T_{\text{c}}$ and other thermodynamic
properties. In order to obtain a rough estimate of how such strong variations
of the magnetic coupling constants will affect the Curie temperature, compared
to using only an average coupling, obtained, e.g. from an effective medium
treatment of chemical disorder, we perform Monte Carlo simulations for a
simple Heisenberg model with random Gaussian-distributed coupling constants.
We find that variations of the same order as obtained in our DFT calculations
for FeNi can lead to a reduction of $T_{c}$ of around 10 %.

In the following, we first describe the computational method we use to obtain
magnetic exchange couplings, and then present our results for both ordered and
partially disordered FeNi.

##  II Computational method

###  II.1 Magnetic exchange interactions

Within the (classical) Heisenberg model, the energy of a magnetic system is
expressed as a sum over pairwise (bilinear) interactions between localized
magnetic moments:

$E=-\frac{1}{2}\sum_{i\neq j}J_{ij}\textbf{S}_{i}\cdot\textbf{S}_{j}\quad.$
(1)

We use the convention that $\textbf{S}_{i}$ is a normalized vector describing
only the direction of the magnetic moment at site $i$.

It is well known that, for an itinerant magnetic material such as FeNi, the
Heisenberg model is not necessarily a good approximation Moriya (1985), and
Eq. (1) is typically only valid for not too large fluctuations around the
ferromagnetic ground state. Thus, to calculate magnetic exchange couplings,
$J_{ij}$, we use the following well-known equation based on the magnetic force
theorem, which is obtained by considering the energy variation with respect to
infinitesimal rotations of the magnetic moments Liechtenstein  _et al._
(1987):

$J_{ij}=\frac{1}{2\pi}\mathrm{Im}\int_{-\infty}^{\varepsilon_{\mathrm{F}}}\mathrm{d}\varepsilon\sum_{mm^{\prime}m^{\prime\prime}m^{\prime\prime\prime}}\Delta_{i}^{mm^{\prime}}G_{ij,\downarrow}^{m^{\prime}m^{\prime\prime}}(\varepsilon)\Delta_{j}^{m^{\prime\prime}m^{\prime\prime\prime}}G_{ji,\uparrow}^{m^{\prime\prime\prime}m}(\varepsilon)\quad.$
(2)

Here, $\Delta_{i}^{mm^{\prime}}$
($\Delta_{j}^{m^{\prime\prime}m^{\prime\prime\prime}}$) is the local exchange
splitting on site $i$ ($j$), and $G_{ji,\uparrow}^{m^{\prime\prime\prime}m}$
($G_{ij,\downarrow}^{m^{\prime}m^{\prime\prime}}$) is the spin-up (spin-down)
intersite Green’s function. Both quantities are expressed within a tight-
binding-like basis, where each basis-orbital (with index $m$) is localized on
a specific site.

To evaluate the quantities in Eq. (2), we first obtain the electronic
structure from plane-wave-based density-functional theory (DFT) calculations,
and then transform the corresponding Kohn-Sham Hamiltonian into a basis of
localized Wannier functions Rudenko  _et al._ (2013); Korotin  _et al._
(2015); Logemann  _et al._ (2017); Zhu  _et al._ (2020). As described in more
detail in Sec. III.1, we use Wannier functions defined by orbital projection
and subsequent orthonormalization (corresponding to the “initial projections”
in the wannier90 code Mostofi  _et al._ (2014)). This leads to a set of atom-
centered basis orbitals. As shown in Sec. III.1, an excellent representation
of all occupied bands in FeNi can be achieved by using a full set of $s$, $p$,
and $d$ projections for each atom.

In contrast, constructing maximally localized Wannier functions (MLWFs)
Marzari and Vanderbilt (1997) for FeNi, results in a set of Wannier functions
where the Wannier orbitals corresponding to the $s$ and $p$ projections become
localized in between the atoms, and thus cannot be used to evaluate Eq. (2).
This is similar to what has been described for the nearly free-electron-like
bands in fcc Cu (and other $3d$ transition metals), see e.g., Ref. Souza  _et
al._ , 2001. Further details are presented in Sec. III.1.

After a suitable set of Wannier functions has been constructed, we follow the
approach outlined in Ref. Zhu  _et al._ (2020) to obtain exchange couplings
for different pairs of atoms.

###  II.2 Computational details

In order to accommodate both the fully ordered L10 structure of FeNi as well
as some configurations with (partial) chemical disorder, we use an 8-atom
cell, corresponding to a $\sqrt{2}\times\sqrt{2}\times 1$ supercell of the
conventional cubic 4-atom fcc unit cell, or, equivalently, to a $2\times
2\times 1$ supercell of the 2-atom tetragonal primitive unit cell of the L10
structure. For test purposes, we also perform some calculations for the
perfectly ordered L10 structure using the 2-atom primitive unit cell. In all
cases, except where otherwise noted, we fix the lattice parameters and atomic
positions to that of a perfectly cubic fcc lattice with lattice constant
$a=3.56$ Å, and then distribute Fe and Ni atoms over the available sites
within the cell in different ways.

We perform DFT calculations using the Vienna ab initio Simulation package
(VASP) Kresse and Furthmüller (1996), the projector-augmented wave method
(PAW) Blöchl (1994); Kresse and Joubert (1999), and the generalized gradient
approximation according to Perdew, Burke, and Ernzerhof (PBE) Perdew  _et al._
(1996). Brillouin zone integrations are performed using the tetrahedron method
with Blöchl corrections and a $\Gamma$-centered $12\times 12\times 16$ k-point
mesh for the 8-atom cell. The plane wave energy cut-off is set to 550 eV, and
the total energy is converged to an accuracy of $10^{-8}$ eV. Our PAW
potentials include $3p$, $4s$, and $3d$ states in the valence for both Fe and
Ni. All calculations are performed for the ferromagnetically ordered state.

A Wannier representation of the Kohn-Sham Hamiltonian is then obtained using
the wannier90 code Mostofi  _et al._ (2014), using the same $k$-point mesh as
for the DFT calculations. To check the convergence of the calculated magnetic
exchange couplings with respect to the k-point sampling, we perform
calculations using up to $14\times 14\times 18$ k-points and find our results
to be sufficiently converged using a $12\times 12\times 16$ k-point mesh.

To obtain the Curie temperature for the ordered case, using the Heisenberg
model, Eq. (1), with the coupling constants obtained from our DFT
calculations, we perform Metropolis Monte Carlo simulations, as implemented in
the UppASD package Skubic  _et al._ (2008). We consider magnetic exchange
couplings for all pairs of atoms up to a distance of 10$a$ ($\sim$ 25.17 Å).
To accurately determine $T_{\text{c}}$, accounting for potential finite size
effects due to the limited size of our simulation cells, we use the Binder
cumulant method Landau and Binder (2014). Thus, $T_{\text{c}}$ is obtained as
the temperature where the fourth-order Binder cumulants, obtained for three
different cell sizes, cross. We consider cell-sizes of $20\times 20\times 20$,
$30\times 30\times 30$, and $36\times 36\times 36$, relative to the primitive
tetragonal 2-atom cell of the L10 structure. For the test calculations
presented in Fig. 3(b), we use a cell size of $20\times 20\times 20$.

For the model study presented in Sec. III.4, we use a $15\times 15\times 15$
supercell of the conventional fcc cubic cell. We then initialize the magnetic
couplings for all nearest neighbor pairs within this cell individually by
drawing random numbers from a Gaussian distribution with varying standard
deviation, $\sigma>0$, and a mean value of $\mu=1$. For $\sigma=0.0$, the
coupling constants of all nearest neighbor pairs are identical and equal to
the mean of the Gaussian distribution ($\mu=1.0$). The Curie temperature for
each $\sigma>0.0$ is obtained by taking an average over 100 instances of the
Gaussian-distributed magnetic exchange couplings with different random seeds.
Convergence with respect to the number of instances of the randomized system
has been verified by monitoring the cumulative average of $T_{\text{c}}$ with
increasing number of instances over all 100 samples.

##  III Results and Discussion

###  III.1  Construction of Wannier functions

As outlined in Sec. II.1, we start by constructing a set of Wannier functions
from orthonormalized projections on $s$, $p$, and $d$ orbitals for each Fe and
Ni atom in the unit cell. We use an outer energy window ranging from $-10$ eV
up to about 31 eV, which contains all occupied valence bands plus a certain
number of empty bands. Furthermore, in order to accurately reproduce all
occupied bands, we employ a frozen (inner) energy window from $-10$ eV up to
about 1 eV above the Fermi energy. The band dispersion (only for the majority
spin component) obtained from the resulting Wannier functions for the minimal
2-atom cell of the fully ordered L10 structure is shown in Fig. 1, together
with the underlying Kohn-Sham bandstructure.

Figure 1: Majority spin bandstructure of chemically ordered L10-FeNi. Blue and
green dots represent the Kohn-Sham and Wannier-interpolated bands,
respectively. Outer and inner energy windows, $E_{\text{o}}$ and
$E_{\text{i}}$, are indicated on the right side of the plot. The Fermi level
defines zero energy.

One can see that all bands below $\sim$ 10 eV are well described by the
Wannier-interpolated bands. Some weak oscillations can be seen in the lowest
lying, free-electron-like band around the $\Gamma$-point. These are due to the
fact that the $k$-mesh along the high-symmetry lines used to obtain the
interpolated bandstructure is much finer than the homogeneous $k$-mesh used to
construct the Wannier functions, and that, in order to obtain atom-centered
Wannier functions, we do not apply the usual “disentanglement procedure” to
obtain an optimally $k$-connected subspace. We note that the calculation of
the magnetic exchange couplings is based on the original homogeneous $k$-point
mesh, where the Wannier bands are identical to the DFT Kohn-Sham bands by
construction. For the minority-spin bands (not shown), we obtain a similar
good agreement between the Wannier-interpolated and the occupied Kohn-Sham
bands.

###  III.2  Magnetic interactions in L10-ordered FeNi

Figure 2: Magnetic exchange couplings for (a) Fe-Fe, (b) Fe-Ni, and (c) Ni-Ni
pairs as a function of atomic distance in the ordered L$1_{0}$ phase of FeNi
as calculated in this work (black dots) and by Edström et al. Edström  _et
al._ (2014) (green triangles).

Next, we calculate the magnetic exchange couplings for the ordered
L$1_{0}$-phase of FeNi, using the method described in Sec. II.1. In Fig. 2, we
compare the magnetic exchange couplings obtained in the present study (shown
by the black dots) with those calculated by Edström et al. Edström  _et al._
(2014) (shown by the green triangles). The different sub-panels show the
couplings corresponding to Fe-Fe, Fe-Ni, and Ni-Ni pairs as a function of
distance. Overall, there is very good agreement between the two data-sets,
except for a few cases discussed further below. One can see that the Fe-Fe
couplings are strongest and rather long-ranged while the Fe-Ni and Ni-Ni
interaction is weaker and decays rather quickly with distance.

Note that for certain distances (e.g., for the one corresponding to second
nearest-neighbour distance, $d_{ij}=a_{0}=3.56$ Å), two distinct values for
$J_{ij}$ are obtained, depending on the orientation (in-plane versus out-of-
plane) of the corresponding pair relative to the tetragonal axis, i.e., the
axis defined by the long range order. In particular for the second nearest
neighbor Fe-Fe interaction, this difference is rather large. This already
indicates a strong configuration dependence of the magnetic coupling, which
will be further analyzed in Sec. III.3 using supercells with partial chemical
disorder.

For the second nearest neighbor Ni-Ni pairs, the coupling is rather weak and
thus the (absolute) difference between in-plane and out-of-plane coupling is
small. Furthermore, one can see that Fe-Ni couplings for certain distances are
“missing” (e.g., corresponding to second nearest neighbors on the fcc lattice,
$d_{ij}=a_{0}=3.56$ Å). This is due to the arrangement of Fe and Ni atoms in
the underlying L$1_{0}$ structure.

The good agreement between our results and the calculations of Edström et al.
(Ref. Edström  _et al._ , 2014) is remarkable, since rather different
electronic structure methods, involving different approximations and basis
sets, have been employed. Furthermore, our calculations are based on a
metrically cubic fcc lattice, whereas Edström et al. have used tetragonal
lattice vectors with a slightly different $c/a$ ratio. The biggest difference
is observed for the second nearest neighbor (in-plane) Fe-Fe coupling and the
first nearest neighbor Fe-Ni coupling, for which Edström et al. obtain a
noticeably stronger coupling (by around 7 meV and 5 meV, respectively). Our
calculated coupling constants also appear to be in good agreement with more
recent calculations by Tian et al. Tian  _et al._ (2020) 111Note that Tian et
al. use a different normalization in the Heisenberg model and therefore the
absolute values in Ref. Tian  _et al._ , 2020 are different from ours.

Figure 3: (a) Magnetic coupling constants as a function of atomic distance
multiplied with the cube of the corresponding atomic distance, $J_{ij}\cdot
d_{ij}^{3}$, calculated for the ordered phase of FeNi. Coupling constants
corresponding to Fe-Fe, Fe-Ni, and Ni-Ni pairs are shown by black, blue, and
green markers, respectively. (b) Ferromagnetic Curie temperature
$T_{\text{c}}$, obtained from Monte Carlo simulations, as a function of the
cutoff distance used for the Fe-Fe interactions (see text).

From Fig. 2(a) one can see that there are noticeable Fe-Fe couplings even for
rather large atomic distances. To further analyze this distance dependence,
Fig. 3(a) shows the same coupling constants as in Fig. 2, but multiplied with
the cube of the corresponding inter-atomic distances, i.e. $J_{ij}\cdot
d_{ij}^{3}$. For the Fe-Ni and Ni-Ni couplings, the corresponding data-points
still converge quickly towards zero for large distances, which means that
these couplings decay with distance faster than $d_{ij}^{-3}$. On the other
hand, for the Fe-Fe pairs, one can see that the data-points oscillate and do
not seem to decay even for very long atomic distances. As pointed out in
previous works Edström  _et al._ (2014); Belozerov  _et al._ (2020), this
indicates an approximate $d_{ij}^{-3}$ dependence of the Fe-Fe interaction
which, together with the oscillatory behavior, is typical for metals with
RKKY-like exchange interactions.

Fig. 3(b) shows the ferromagnetic Curie temperature, $T_{\text{c}}$, obtained
from Monte Carlo simulations of the Heisenberg model, where the Fe-Fe
interactions are considered only up to a certain maximum inter-atomic
distance. For the Fe-Ni and Ni-Ni pairs, all calculated coupling constants
have been included, i.e., up to a very large distance of $10a_{0}$. Note that
due to the fast decay of these couplings, the following results should be
unaffected by the specific cutoff-distance used for the Fe-Ni and Ni-Ni
interactions. It can be seen that, due to the long range of the Fe-Fe
interaction, the calculated $T_{c}$ exhibits strong variations as function of
the Fe-Fe cutoff distance, but seems to converge to a value around 700 K once
all Fe-Fe interactions up to a distance of around 8 Å are taken into account.
This shows that, in order to obtain a reliable estimate of $T_{\text{c}}$ in
this system, it is essential to include Fe-Fe interactions up to rather large
distances.

Based on these test calculations, we now obtain an accurate estimate for
$T_{\text{c}}$ from Monte Carlo simulations of the Heisenberg model including
all calculated coupling constants up to a maximum distance of $10a_{0}$ and
then perform a Binder cumulant analysis, as described in Sec. II.2. We obtain
a value of $T_{c}=736$ K. Note that for the test calculations shown in Fig.
3(b), $T_{\text{c}}$ is obtained simply from the peak position of the
calculated temperature dependence of the specific heat, and thus differs
somewhat from the more accurate value obtained via the Binder cumulants. Our
calculated $T_{\text{c}}$ agrees well with the value of about 780 K obtained
by Tian et al. Tian  _et al._ (2020), whereas the $T_{\text{c}}$ of 916 K
obtained by Edström et al. Edström  _et al._ (2014) is noticeably higher. This
is due to the stronger second nearest neighbor Fe-Fe and first nearest
neighbor Fe-Ni coupling constants obtained in Ref. Edström  _et al._ , 2014
(see Fig. 2). Note that both our and the value of Tian et al. are lower than
the experimental $T_{c}$ of the disordered system and thus seem to
underestimate the “true” Curie temperature of the ordered state.

###  III.3 Magnetic interactions for (partially) disordered FeNi

Figure 4: Different ordered and (partially) disordered configurations
considered in this work, depicted in the $2\times 2\times 1$ supercell
relative to the primitive tetragonal cell of the L$1_{0}$ structure: (a)
ordered, (b) 1-pair-exchanged, and (c) 2-pairs-exchanged. Fe and Ni atoms are
represented by red and green spheres, respectively. Dotted lines in (a)
indicate the underlying conventional cubic cell of the fcc lattice.

We now investigate the effect of chemical disorder on the magnetic exchange
couplings by starting from the fully ordered case, and then successively
exchanging the positions of one or two pairs of Fe and Ni atoms within an
8-atom supercell, resulting in the two configurations shown in Fig. 4(b) and
(c). Note that all other configurations that can be created by exchanging one
Fe-Ni pair in this 8-atom supercell are equivalent to the one shown in Fig.
4(b), whereas several distinct configurations can be created by exchanging two
Fe-Ni pairs. For simplicity we limit our study to the configuration depicted
in Fig. 4(c). In the following, we refer to these two configurations as
“1-pair-exchanged” and “2-pairs-exchanged”, respectively.

Figure 5: (a) Depiction of the shared first nearest and first/second nearest
neighbor environment for a pair of nearest neighbor sites (indicated as red
colored spheres) in an fcc lattice. The conventional cubic cell is also shown
for clarity. (b) Schematic showing the occupation of the shared first and
first/second nearest neighbour environments for all nearest neighbor Fe-Fe
pairs that are present in the ordered [(i)], 1-pair-exchanged [(ii)-(v)], and
2-pairs-exchanged configurations [(vi)-(vii)] along with the corresponding
values of the magnetic coupling constant $J$.

In Fig. 5(b), we list the magnetic exchange couplings obtained for all
inequivalent first nearest neighbor Fe-Fe pairs in the three different
configurations. The corresponding local atomic environments are also indicated
and will be discussed further below. One can see that the calculated values
vary drastically, from 4.7 meV to 36.5 meV, while the corresponding value in
the fully ordered structure is 26.5 meV. We note that, if the nearest neighbor
Fe-Fe coupling would be completely configuration-independent, then all values
listed in Fig. 5(b) would be identical. The large variation of the Fe-Fe
nearest neighbor coupling in the different cases thus shows that the local
chemical environment has a significant influence on the magnitude of the
magnetic exchange interactions in FeNi. As already discussed in Sec. I, such
variations are not captured by effective medium methods such as the CPA, which
are often used to model chemical disorder in alloys.

Table 1: First nearest-neighbour Fe-Fe magnetic exchange interaction (in meV) for relaxed and unrelaxed 1-pair-exchanged and 2-pairs-exchanged configurations. | Unrelaxed | Relaxed
---|---|---
1-pair-exchanged | 7.4 | 10.1
| 36.5 | 38.8
| 19.1 | 19.9
| 4.7 | 6.8
2-pairs-exchanged | 23.1 | 21.9
| 14.7 | 15.5

The magnetic exchange couplings shown in Fig. 5(b) are obtained by decorating
the sites within a perfect fcc lattice in different ways with Fe and Ni atoms,
without allowing the atomic positions to relax within the resulting lower
symmetry. In order to assess the effect of such relaxations, we now
recalculate the magnetic coupling constants for the 1-pair-exchanged and
2-pairs-exchanged configurations after allowing all atomic positions to relax,
while still keeping the lattice vectors of the supercell fixed. The results
for the first nearest neighbor Fe-Fe couplings are shown in Table 1 and are
compared to the corresponding values for the unrelaxed case. It can be seen
that the relaxation leads to changes in the magnetic coupling constants of up
to about 2 meV, but the effect is clearly significantly weaker than the effect
due to the different chemical environments.

The strong configuration dependence of the first nearest Fe-Fe coupling raises
the question of whether it is possible to identify simple rules on how the
strength of this coupling depends on the local atomic environment. For this
purpose, we analyze the distribution of atoms on the sites that are closest
neighbors to _both_ Fe atoms forming the pair under consideration. As shown in
Fig. 5(a) there are four such sites, which form a square in the mid-plane
perpendicular to the line connecting the two coupled sites. These four sites
represent the minimal environment to be considered in any model describing the
configuration dependence of the first nearest neighbor coupling. The next
“shell” around the coupled Fe-Fe pair is formed by those sites that are first
nearest neighbors to one of the coupled sites and second nearest neighbors to
the other site. There are again four such sites, which form a rectangle in the
plane parallel to the Fe-Fe distance vector and perpendicular to the plane
formed by the common first nearest neighbors (see Fig. 5(a)).

In Fig. 5, we schematically depict the occupation of both the shared first
nearest and the shared first/second nearest neighbor sites for all the
inequivalent nearest neighbor Fe-Fe pairs included in the ordered, 1-pair-
exchanged, and 2-pairs-exchanged configurations, along with the value of the
corresponding magnetic exchange couplings. It is obvious that the shared first
nearest neighbor environment is not sufficient to classify the different
coupling constants, since, e.g., cases (iii), (v), and (vi) all have an
equivalent shared first nearest neighbor environment but exhibit vastly
different magnetic coupling constants (including both the highest and lowest
calculated values of 36.5 meV and 4.7 meV). The same holds for cases (i) and
(ii).

Considering both the shared first and first/second nearest neighbor
environment, all inequivalent Fe-Fe pairs contained in our three
configurations exhibit different local environments, which is in principle
compatible with a local model for the exchange coupling based on this
environment. However, to really establish or disprove such a model, one has to
consider much larger supercells, that allow to sample more configurations, and
also include cases with identical first/second nearest neighbor environment
but different further neighbor environment. This would require an excessive
computational effort. Considering that, in general, the applicability of a
short-range local model for an itinerant magnetic system such as FeNi is
rather questionable, we therefore refrain from sampling further couplings
using larger and larger supercells. Instead, we try to estimate the effect of
a strong configuration dependence of the magnetic coupling constants on the
Curie temperature of a disordered magnetic system using a simple Heisenberg
model with random couplings.

###  III.4 Model study with random couplings

The very high sensitivity of the magnetic coupling constants on the specific
chemical environment, and the high computational effort to fully resolve this
configuration dependence (if at all possible), represents a big obstacle for
the reliable estimation of magnetic ordering temperatures for disordered
itinerant magnets such as FeNi from first principles calculations. In the
following, we therefore employ a strongly simplified model to obtain a rough
estimate of how the strong configuration-dependent variations of the exchange
couplings can affect the Curie temperature of a disordered magnetic system, in
comparison to the Curie temperature obtained using only configuration-
independent “average” magnetic coupling constants.

Specifically, we consider a Heisenberg model with only nearest neighbor
interactions on an fcc lattice, and we approximate the configuration-dependent
variations of the magnetic coupling constants by a Gaussian-distributed random
variable, where the mean value of the Gaussian distribution represents the
average coupling constant, and its standard deviation, $\sigma$, quantifies
the configuration-dependent variations. We then perform temperature-dependent
Monte Carlo simulations as outlined in Sec. II.2, and analyze how the obtained
Curie temperature depends on $\sigma$, i.e., on the strength of the variation
in the magnetic coupling constants.

Figure 6: (a) Susceptibility, and (b) specific heat as a function of
temperature (in relative units) obtained from Monte Carlo simulations for a
simple Heisenberg model with only nearest neighbor couplings on an fcc lattice
for the case when all coupling constants are identical ($\sigma=0$). $T_{0}$
is the critical temperature for $\sigma=0$.

Fig. 6 shows the calculated susceptibility ($\chi$) and specific heat
($C_{v}$) as a function of temperature for this model for the case when all
coupling constants are identical to the average one ($\sigma=0$). One can see
that both $\chi$ and $C_{v}$ exhibit clear peaks at the same critical
temperature $T_{0}$. However, since the peak in the susceptibility appears
much sharper than the one in the specific heat, in the following we use the
peak value of $\chi$ to accurately determine the Curie temperature as function
of $\sigma$.

Figure 7: Curie temperature $T_{c}$ of the Heisenberg model with Gaussian-
distributed coupling constants as a function of the standard deviation,
$\sigma$ (defined relative to the mean value of the Gaussian distribution).
$T_{0}=T_{c}(\sigma=0)$ is the Curie temperature obtained for the average
coupling, and the error bars indicate the standard deviation of the mean
obtained by averaging over 100 instances with different random seeds for each
$\sigma$ (see text). The black dashed line is a quadratic fit to the data.

Fig. 7 shows the variation of the average Curie temperature, $T_{\text{c}}$,
relative to $T_{0}$, as a function of the standard deviation $\sigma$ of the
Gaussian distributed random magnetic coupling constants. Note that $\sigma$
corresponds to a Gaussian distribution with a mean equal to the average
magnetic coupling, and is therefore defined relative to this average coupling.
One can see that $T_{\text{c}}$ decreases with increasing standard deviation
$\sigma$, i.e. with increasing “randomness” of the magnetic exchange
couplings, and that this decrease can be fitted well with a quadratic
dependence. This indicates that using configuration-independent “average”
magnetic exchange couplings obtained by effective medium approaches such as
CPA are expected to overestimate the critical temperature of a disordered
system.

If we take the seven different values for the configuration-dependent nearest
neighbor interaction obtained from the ordered, the 1-pair-exchanged, and the
2-pairs exchanged configuration and evaluate the empirical standard deviation,
we obtain $\sigma\approx 0.6$ (relative to an average value of 18.9 meV).
Comparing this with the data shown in Fig. 7, this would correspond to a
reduction of $T_{c}$ of about 10 % compared to the value obtained using only
an average coupling strength. We note that this is clearly a very naive
estimation based on a very small number of samples, but it shows that a strong
configuration dependence of the magnetic coupling, as observed for FeNi in
Sec. III.3, can indeed lead to a noticeable reduction of $T_{c}$ compared to
that obtained from average effective medium couplings.

Note that the effect of random magnetic exchange couplings on the magnetic
transition temperature of a Heisenberg model with nearest and next-nearest
antiferromagnetic coupling on a two-dimensional square lattice was studied
previously by Li et al. Li  _et al._ (2015). Thereby, the nearest neighbor
coupling was obtained from a homogeneous random distribution within a given
range, and it was also found that the transition temperature decreases with
increasing variation of the random magnetic exchange couplings.

##  IV Summary and Conclusions

In summary, we have investigated the effect of chemical disorder on the
magnetic exchange couplings in L10-FeNi using first principles DFT
calculations. Thereby, we have used supercells with different atomic
distributions, to specifically include effects due to the symmetry-broken
local environments around the individual atoms that are not included in
effective medium approaches such as CPA. We find that such effects can lead to
rather large variations of the magnetic exchange couplings, exemplified by our
analysis of the nearest neighbor Fe-Fe coupling, which exhibits values ranging
from 4.7 meV up to 36.5 meV, with a value of 26.5 meV obtained for the fully
ordered structure.

Our analysis of the shared first and shared first/second nearest neighbor
chemical environments of the different inequivalent Fe-Fe nearest neighbor
pairs included in our supercells indicates that the strength of the couplings
is affected by long range effects that go beyond the closest neighbor shell,
as can be expected for an itinerant metallic material such as FeNi. The lack
of a simple relation that relates the variation of an individual $J_{ij}$
(relative to the corresponding average value) to its local chemical
environment makes it very challenging to incorporate effects beyond CPA in the
calculation of $T_{\text{c}}$. However, our simple model study using random
Gaussian-distributed nearest neighbor couplings on an fcc lattice suggests
that such local variations of the coupling constants can lead to a reduction
of $T_{\text{c}}$ of up to 10 %, compared to that obtained using only an
average coupling.

Thus, there is a hierarchy of effects that in general tend to reduce the
magnetic ordering temperature in random alloys such as partially ordered FeNi.
First, the “average” coupling strength is affected by the chemical disorder.
For example, Tian et al. obtain a reduction of around 30-35 % of the nearest
neighbor Fe-Fe coupling obtained within CPA for the disordered system compared
to the fully ordered case Tian  _et al._ (2020). This is in principle
consistent with our supercell calculations, where the simple average of this
coupling over all inequivalent Fe-Fe pairs in our two disordered
configurations gives about 17.6 meV, i.e., a reduction by about 35 % compared
to the ordered case. However, one should note that this average is based on
only very few samples. Second, the random connectivity between different
magnetic atoms (Fe and Ni in our case) will also affect the Curie temperature
relative to the ordered case, which corresponds to a very regular network of
Fe-Fe, Fe-Ni and Ni-Ni bonds. Finally, the variation of the coupling strength
according to the specific chemical environment around the atoms in the random
alloy can lead to a further decrease of $T_{\text{c}}$, as indicated by our
simple model. It is mainly this last effect that we have quantified within
this work.

Furthermore, it appears that our first-principles-based results are
underestimating the (hypothetical) Curie temperature of the fully ordered
system. It is unclear whether this underestimation is related to the general
applicability of the (classical) Heisenberg model to FeNi or whether it is
caused by deficiencies of the generalized gradient approximation in the
underlying DFT calculations (or other approximations in the method). However,
we note that the comparison with the results obtained by Edström et al.
Edström  _et al._ (2014) also demonstrates that moderate changes in specific
calculated coupling constants can lead to rather strong differences in the
predicted Curie temperatures.

###### Acknowledgements.

We are grateful to Alexander Edström for many helpful discussions and for
providing the data from Ref. Edström  _et al._ , 2014 included in Fig. 2. We
also thank Maximilian Merkel and Alberto Carta for help with using the
wannier90 code. This work was supported by ETH Zürich. Calculations were
performed on the cluster “Piz Daint”, hosted by the Swiss National
Supercomputing Centre, and the “Euler” cluster of ETH Zürich.

##  Appendix A Orbital decomposition of nearest neighbor Fe-Fe coupling

Eq. (2) in principle allows for a decomposition of the coupling constants
$J_{ij}$ into different orbital contributions. However, for the case of
L10-FeNi, the quantities appearing in Eq. (2), in particular the exchange
splitting $\Delta_{i}^{mm^{\prime}}$, contain off-diagonal elements mixing the
$d$ and $s$ type Wannier orbitals. Nevertheless, by restricting the summation
in Eq. (2) to only the diagonal elements $\Delta_{i}^{mm}$ and considering
only contributions from the $d$-type orbitals, one can obtain a “$d$-only”
contribution to the magnetic coupling.

In this way, we obtain a $d$-only contribution to the Fe-Fe coupling of 29.6
meV (second row in Table 2). This value is larger than the full value of 27.2
meV (first row in Table 2), obtained by considering all contributions in Eq.
(2), showing that here the combination of $sp$ and orbitally mixed terms leads
to a small negative contribution to $J_{ij}$ 222The small difference to the
value of 26.5 meV listed in Fig. 5 is due to the the different unit cells
(8-atom cell versus 2-atom cell) used in the calculations.. Furthermore, it
demonstrates that, as probably expected, the main contribution to the magnetic
coupling stems from the $d$ orbitals.

Table 2: Magnetic exchange couplings, $J_{ij}$, for the nearest-neighbor Fe-Fe interaction in L10-FeNi obtained using different sets and subsets of Wannier functions. The first column indicates the different sets of Wannier functions that have been constructed, while the second column indicates which terms are considered in Eq. 2 when evaluating the corresponding $J_{ij}$. All calculations are performed for the primitive 2-atom unit cell. Wannier set | Terms in Eq. 2 |  $J_{ij}$ (meV)
---|---|---
$d+s+p$ (projections) | all | 27.2
$d+s+p$ (projections) |  $d$-only | 29.6
$d$ (projections) | all | 35.7
$d$ (MLWFs) | all | 43.4

This raises the question of whether it would be sufficient to consider only
the $d$-bands in the first place, i.e., construct a smaller set of Wannier
functions describing only the $d$ bands, and still obtain a good estimate for
$J_{ij}$. To test this hypothesis, we construct two additional sets of Wannier
functions, where we include only five $d$-orbitals per atom. For the first
set, we obtain the Wannier functions from orthonormalized atomic projections
as before, while for the second set we perform a subsequent minimization of
the quadratic spread functional to obtain MLWFs. In both cases we use an
(outer) energy window ranging from $-10$ eV to about 5 eV above the Fermi
level and obtain a set of atom-centered $d$-like orbitals suitable to evaluate
Eq. (2). In both cases also the Wannier-interpolated bands resemble the DFT
band-structure in the energy range of the $d$-bands, i.e., between
approximately $-5$ meV and the Fermi level for the majority spin channel.

The magnetic coupling constants for the nearest-neighbour Fe-Fe coupling
obtained from these two additional sets of Wannier functions are listed in the
third and fourth row of Table. 2. It can be seen that the corresponding values
(in particular for the set of MLWFs) are significantly larger than the
$d$-only contribution obtained from the full description using also $s$ and
$p$ bands. This shows that, even though the $d$-orbitals make up the main
contribution to the magnetic coupling constants, it is nevertheless important
to include $s$ and $p$ states to accurately account for their effect on the
$d$ band dispersion. We note that, due to the entanglement of $d$ and $sp$
contributions in the bandstructure of FeNi, the $d$ subset of the full $d+s+p$
Wannier basis and the two different $d$-only Wannier sets (projected and
MLWFs) are all describing slightly different subspaces of the occupied Kohn-
Sham states.

Only the $d+s+p$ Wannier set results is a complete and accurate description of
all occupied bands in FeNi, and therefore only the corresponding value of
$J_{ij}$ should be considered as “correct” (or most accurate). Nevertheless,
our analysis raises the question of a potential basis set dependence of the
magnetic coupling constants, for example in cases where a complete description
of all occupied bands can be achieved using different sets of Wannier
functions, e.g., corresponding to different degrees of localization. In the
present case, a more systematic analysis is hindered by the strong
entanglement of bands and the fact that spread minimization on the full set of
Wannier functions leads to orbitals that are not atom-centered. We therefore
leave this question open for future research.

## References

  * Paulevé  _et al._ (1968) J. Paulevé, A. Chamberod, K. Krebs, and A. Bourret, Journal of Applied Physics 39, 989 (1968).
  * Kojima  _et al._ (2014) T. Kojima, M. Mizuguchi, T. Koganezawa, M. Ogiwara, M. Kotsugi, T. Ohtsuki, T.-Y. Tashiro, and K. Takanashi, J. Phys. D Appl. Phys. 47, 425001 (2014).
  * Lewis  _et al._ (2014a) L. H. Lewis, F. E. Pinkerton, N. Bordeaux, A. Mubarok, E. Poirier, J. I. Goldstein, R. Skomski, and K. Barmak, IEEE Magnetics Letters 5, 1 (2014a).
  * Lewis  _et al._ (2014b) L. H. Lewis, A. Mubarok, E. Poirier, N. Bordeaux, P. Manchanda, A. Kashyap, R. Skomski, J. Goldstein, F. E. Pinkerton, R. K. Mishra, R. C. K. Jr, and K. Barmak, Journal of Physics: Condensed Matter 26, 064213 (2014b).
  * Edström  _et al._ (2014) A. Edström, J. Chico, A. Jakobsson, A. Bergman, and J. Rusz, Physical Review B 90, 014402 (2014).
  * Werwiński and Marciniak (2017) M. Werwiński and W. Marciniak, Journal of Physics D: Applied Physics 50, 495008 (2017).
  * Cui  _et al._ (2018) J. Cui, M. Kramer, L. Zhou, F. Liu, A. Gabay, G. Hadjipanayis, B. Balasubramanian, and D. Sellmyer, Acta Materialia 158, 118 (2018).
  * Tian  _et al._ (2019) L.-Y. Tian, H. Levämäki, O. Eriksson, K. Kokko, Á. Nagy, E. K. Délczeg-Czirják, and L. Vitos, Scientific Reports 9, 8172 (2019).
  * Tian  _et al._ (2020) L.-Y. Tian, O. Eriksson, and L. Vitos, Scientific Reports 10, 14766 (2020).
  * Izardar and Ederer (2020) A. Izardar and C. Ederer, Phys. Rev. Materials 4, 054418 (2020).
  * Tian  _et al._ (2021) L.-Y. Tian, O. Gutfleisch, O. Eriksson, and L. Vitos, Scientific Reports 11, 5253 (2021).
  * Pauleve  _et al._ (1962) J. Pauleve, D. Dautreppe, J. Laugier, and L. Néel, J. Phys. Radium 23, 841 (1962).
  * Néel  _et al._ (1964) L. Néel, J. Pauleve, R. Pauthenet, J. Laugier, and D. Dautreppe, Journal of Applied Physics 35, 873 (1964).
  * Goto  _et al._ (2017) S. Goto, H. Kura, E. Watanabe, Y. Hayashi, H. Yanagihara, Y. Shimada, M. Mizuguchi, K. Takanashi, and E. Kita, Scientific Reports 7 (2017).
  * Makino  _et al._ (2015) A. Makino, P. Sharma, K. Sato, A. Takeuchi, Y. Zhang, and K. Takenaka”, Scientific Reports 5 (2015).
  * Saito  _et al._ (2019) M. Saito, H. Ito, Y. Suzuki, M. Mizuguchi, T. Koganezawa, T. Miyamachi, F. Komori, K. Takanashi, and M. Kotsugi, Applied Physics Letters 114, 072404 (2019).
  * Kojima  _et al._ (2014) T. Kojima, M. Ogiwara, M. Mizuguchi, M. Kotsugi, T. Koganezawa, T. Ohtsuki, T.-Y. Tashiro, and K. Takanashi, Journal of Physics: Condensed Matter 26, 064207 (2014).
  * Reuter  _et al._ (1989) K. B. Reuter, D. B. Williams, and J. I. Goldstein, Metallurgical Transactions A 20, 711 (1989).
  * Albertsen  _et al._ (1978) J. F. Albertsen, J. M. Knudsen, and G. B. Jensen, Nature 273, 453 (1978).
  * Petersen  _et al._ (1977) J. Petersen, M. Aydin, and J. Knudsen, Physics Letters A 62, 192 (1977).
  * Danon  _et al._ (1979) J. Danon, R. B. Scorzelli, I. S. Azevedo, and M. Christophe-Michel-Lévy, Nature 281, 469 (1979).
  * Danon  _et al._ (1980) J. Danon, R. B. Scorzelli, I. Souza-Azevedo, J. Laugier, and A. Chamberod, Nature 284, 537 (1980).
  * Clarke and Scott (1980) R. S. Clarke and E. R. D. Scott, American Mineralogist 65, 624 (1980).
  * Onodera  _et al._ (1981) S. Onodera, Y. Ishikawa, and K. Tajima, Journal of the Physical Society of Japan 50, 1513 (1981).
  * Wei  _et al._ (2014) Q. Wei, S. A. Gilder, and B. Maier, Physical Review B 90, 144425 (2014).
  * Santos  _et al._ (2015) E. D. Santos, J. Gattacceca, P. Rochette, G. Fillion, and R. Scorzelli, Journal of Magnetism and Magnetic Materials 375, 234 (2015).
  * Soven (1967) P. Soven, Phys. Rev. 156, 809 (1967).
  * Gyorffy (1972) B. L. Gyorffy, Phys. Rev. B 5, 2382 (1972).
  * Si  _et al._ (2021) M. Si, A. Izardar, and C. Ederer, “Effect of chemical disorder on the magnetic anisotropy in L$1_{0}$ FeNi from first principles calculations,” (2021), arXiv:2111.12492 [cond-mat.mtrl-sci] .
  * Moriya (1985) T. Moriya,  _Spin Fluctuations in Itinerant Electron Magnetism_ ((Springer-Verlag, Berlin, Heidelberg, 1985).
  * Liechtenstein  _et al._ (1987) A. Liechtenstein, M. Katsnelson, V. Antropov, and V. Gubanov, Journal of Magnetism and Magnetic Materials 67, 65 (1987).
  * Rudenko  _et al._ (2013) A. N. Rudenko, F. J. Keil, M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. B 88, 081405 (2013).
  * Korotin  _et al._ (2015) D. M. Korotin, V. V. Mazurenko, V. I. Anisimov, and S. V. Streltsov, Phys. Rev. B 91, 224405 (2015).
  * Logemann  _et al._ (2017) R. Logemann, A. N. Rudenko, M. I. Katsnelson, and A. Kirilyuk, Journal of Physics: Condensed Matter 29, 335801 (2017).
  * Zhu  _et al._ (2020) X. Zhu, A. Edström, and C. Ederer, Phys. Rev. B 101, 064401 (2020).
  * Mostofi  _et al._ (2014) A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer Physics Communications 185, 2309 (2014).
  * Marzari and Vanderbilt (1997) N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
  * Souza  _et al._ (2001) I. Souza, N. Marzari, and D. Vanderbilt, Physical Review B 65, 035109 (2001).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Physical Review B 54, 11169 (1996).
  * Blöchl (1994) P. E. Blöchl, Physical Review B 50, 17953 (1994).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Physical Review B 59, 1758 (1999).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
  * Skubic  _et al._ (2008) B. Skubic, J. Hellsvik, L. Nordström, and O. Eriksson, Journal of Physics: Condensed Matter 20, 315203 (2008).
  * Landau and Binder (2014) D. P. Landau and K. Binder,  _A Guide to Monte Carlo Simulations in Statistical Physics_, 4th ed. (Cambridge University Press, 2014).
  * Note (1) Note that Tian et al. use a different normalization in the Heisenberg model and therefore the absolute values in Ref. Tian2020 are different from ours.
  * Belozerov  _et al._ (2020) A. S. Belozerov, A. A. Katanin, and V. I. Anisimov, Journal of Physics: Condensed Matter 32, 385601 (2020).
  * Li  _et al._ (2015) W. C. Li, X. Song, J. J. Feng, X. T. Jia, M. Zeng, X. S. Gao, and M. H. Qin, Journal of Applied Physics 118, 013901 (2015).
  * Note (2) The small difference to the value of 26.5+.1667emmeV listed in Fig. 5 is due to the the different unit cells (8-atom cell versus 2-atom cell) used in the calculations.