#  Ab initio description of magnetic properties of spin-glass pyrochlore
NaSrMn2F7

Mohammad Amirabbasi<EMAIL_ADDRESS>Department of Physics, Isfahan
University of Technology (IUT), Isfahan 84156-83111, Iran  Faculty of Physics,
University Duisburg-Essen, 47057 Duisburg, Germany

###### Abstract

The magnetic and critical properties of manganese pyrochlore fluoride
NaSrMn2F7 have been determined which exhibits a glass transition due to charge
disorder. Spin-Hamiltonian calculations show that only the Heisenberg exchange
plays a significant role in determining the magnetic order of the ground
state. The Monte Carlo simulation detects that without considering weak bond
disorder ($\delta J$), no phase transition is observed at temperature
variation of specific heat. The $\delta J$=0.34 meV has been calculated by
using Monte Carlo simulations in such a way that the obtained spin freezing
temperature will be equal to the experimental one. The microscopic calculation
of $\delta J$ within DFT for different Na/Sr distributions proves that $\delta
J$ should be around 0.31 meV. The (h0l) planes density plot of elastic neutron
scattering structure function at low-temperature shows that the ground state
has a degenerate manifold due to the different numbers of short-range degrees
of freedom. The data collapsing of specific heat for different lattice sizes
by using finite-size scaling theory confirms that the obtained freezing
temperature by considering $\delta J$=0.34 meV is in good agreement with
experimental measurements.

###### pacs:

71.15.Mb, 75.40.Mg, 75.10.Hk, 75.30.Gw

##  I introduction

Spin-glass (SG) is a type of magnetic phase in which the spatial arrangement
of spins at low-temperatures is completely random Mydosh (1993); Binder and
Young (1986); Taniguchi  _et al._ (2009); Shintani and Tanaka (2006).

Figure 1: (Color online) Crystal structure of SG pyrochlore NaSrMn2F7 with
corner-shared tetrahedrons for (Left) A-site (Na-Sr atoms) and (right) B-site
(Mn atoms). Green, yellow and purple spheres are Na1+, Sr2+ and Mn2+,
respectively. $x$=0.3331 is experimental value which is reported from the
X-ray measurement Sanders  _et al._ (2016). The experimental lattice constant,
bond length and bond angle in this material are 80.9 and 99.1 (∘), 2.06(Å) and
10.8012(3)(Å), respectively Sanders  _et al._ (2016).

The presence of frustration (coexisting and competing ferromagnetic and
antiferromagnetic interactions) and randomness (such as chemical disorders or
weak bond disorders) are essential for glass formation Mydosh (1993); Binder
and Young (1986); Kawamura and Taniguchi (2015). These materials undergo a
second-order phase transition which makes the study of SG materials attractive
to scientists Kawamura and Taniguchi (2015). Experimental evidences to detect
such a phase are the bifurcation of the field cooled (FC) and zero field
cooled (ZFC) of DC susceptibility, the frequency dependence of AC
susceptibility, no magnetic saturation even at low-temperature and high
magnetic field and diffusive peaks inelastic neutron scattering pattern Mydosh
(1993); Taniguchi  _et al._ (2009); Vincent and Dupuis (2018). In recent
years, some geometrically frustrated pyrochlores have shown the SG phase
Mitsumoto  _et al._ (2020); Shinaoka  _et al._ (2011). These materials are
composed of a network of corner-sharing tetrahedra Subramanian  _et al._
(1983); Gardner  _et al._ (2010). If the exchange interaction between the
nearest neighbor is antiferromagnetic, the system fails to reach the ground
state with a global minimum. This condition is called geometric frustration
Subramanian  _et al._ (1983); Gardner  _et al._ (2010); Lacroix  _et al._
(2011); Reig-i Plessis and Hallas (2021). Despite many studies and research,
there are still many questions about what is the truly thermodynamic SG
transition Binder and Young (1986). Determination of the ground state of
classical Heisenberg antiferromagnetic SG pyrochlores in which only Heisenberg
interactions play a key role can be a good candidate to answer these
questions.

In 1997, with the discovery of the spin-ice phase in Ho2Ti2O7 Harris  _et al._
(1997), a study on the magnetic properties of oxide pyrochlores increased
Slobinsky  _et al._ (2021); Moessner and Sondhi (2003); Marlton  _et al._
(2021); Yahne  _et al._ (2021); Shi  _et al._ (2021); Fennell  _et al._
(2009); Bramwell and Gingras (2001). One of the main problems in the study of
these materials is the difficult conditions of their crystal growth. As a
consequence, the availability of large crystal pyrochlores in large sizes has
made it impossible. There is no such problem for pyrochlore fluoride with
general chemical formula A2B2F7 in which Na1+/Sr2+(Ca2+) atoms occupy A-site
position randomly and magnetic atoms such as Mn2+, Fe2+, Co2+ and Ni2+ occupy
B-site Sanders  _et al._ (2016); Krizan and Cava (2014, 2015a, 2015b). The
result of this random occupation is the presence of the $\delta J$ (difference
in the bond angle between the magnetic ions). Therefore, all of these
materials exhibit SG properties at low-temperatures Shinaoka  _et al._ (2011);
Andreanov  _et al._ (2010). In this paper, I intend to investigate the
magnetic and critical properties of NaSrMn2F7 which has the largest spin
moment among the pyrochlore fluoride by first-principle methods in the
framework of DFT and Monte Carlo (MC) simulation. In this material, the Mn2+
terminates in 3d5. Therefore, S=5/2 and the orbital moment has a small
contribution to the effective moment (6.25 $\mu_{B}$/Mn Sanders  _et al._
(2016)). The absence of orbital moment and its associated complexities make
this material a suitable candidate to study the effect of bond weak disorder
in the phase transition of SG ground state.

The lattice of the NaSrMn2F7 consists of two separate sublattices with corner-
shared tetrahedra (Fig. 1). In one sublattice, the Na1+/Sr2+ ions and in the
other the magnetic ion (Mn2+) sit on the vertices of each tetrahedra. The
space group of this pyrochlore is Fd-3m (No. 227). Accordingly, the Na1+/Sr2+,
Mn2+ and F1- occupy the 16d (0.5, 0.5, 0.5), the 16c (0, 0, 0), 48f ($x$, 1/8,
1/8) and the 8b (3/8, 3/8, 3/8), respectively. The only variable parameter in
this structure is the $x$-positional parameter. This parameter controls the
bond angle (Mn-F-Mn) and bond length (Mn-F) which plays a decisive role in the
magnetic properties. Experimental experiments show that this material is an
insulator and shows a SG transition at 2.5 K Sanders  _et al._ (2016). The
most important reason for SG transition can be attributed to the change in the
bond angle due to the presence of charge disorder (Na1+/Sr2+) in the lattice
Sanders  _et al._ (2016). The interaction between Mn2+ ions is the
superexchange mediated by F1-. This mechanism is directly related to the bond
length and bond angle. The change in bond angle causes variations in exchange
parameters. As a result, the ground state of the system takes SG order. The
Curie-Weiss temperature ($\theta_{\mathrm{CW}}$) equals $-$89.72 K which
indicates the presence of strong antiferromagnetic interactions between
magnetic ions Sanders  _et al._ (2016).

Microscopic examination of the exchange interactions in this material and its
effect on the formation of the SG phase as well as the study of the critical
properties seem interesting. This paper aims to determine magnetic and
critical properties of NaSrMn2F7 by constructing a spin Hamiltonian model. It
includes various interactions such as the Heisenberg exchange parameters up to
the third neighbor, the bi-quadratic, single-ion anisotropy (SIA) and
Dzyaloshinskii-Moriya interactions (DMI) as obtained by the first-principle
methods. By using this information, one can obtain the critical properties of
this material, including transition temperature and spin correlation and
elastic neutron scattering by the MC simulation. The current calculations show
that only the $J_{1}$ \+ $\delta{J}$ is important in determining the magnetic
properties of this material, and the effect of other interactions is
negligible. As a result, this material is considered as a simple model of SG
and it can be a good candidate for further studies such as understanding the
electron-electron correlation influence on magnetic ground state
Alexandradinata  _et al._ (2020) especially in the SG phase, the effect of
$\delta{J}$ on thermal fluctuations and removing huge degeneracy Andreanov
_et al._ (2010) and what is the border between SG and liquid phase according
to the strength of $\delta{J}$ Moessner and Chalker (1998); Yang  _et al._
(2015); Cepas and Canals (2012).

The paper is structured as follows. In section II the details of the DFT and
MC computational methods have been presented. Section III is devoted to
deriving the spin Hamiltonian and the different aspects of exchange
interactions in determination of magnetic ground state of NaSrMn2F7 and in the
following, the critical properties are discussed. Finally, In section IV a
summary is given.

Figure 2: (Color online) Different collinear spin configurations for
calculation Heisenberg couplings. The purple spheres denote Mn2+. The arrows
indicates the direction of the spins for each atom. From left to right, the
colors of some of the spins are selected differently to distinguish between
different configurations. The ferromagnetic configuration has been chosen as
the reference to calculate the total energy difference of the various
configurations. In this paper, 10 different spin configurations have been used
to make sure that the $J_{1}$ converges with respect to the number of
configurations.

##  II Computational methods

In this paper, the experimental lattice constant and $x$-positional parameter
have been examined. One introduces the classical spin Hamiltonian ($H$) to get
magnetic properties of SG pyrochlore NaSrMn2F7:

$\begin{split}H_{\rm{spin}}=-\frac{1}{2}\sum_{i\neq
j}J_{ij}(\vec{S_{i}}\cdot\vec{S_{j}})+\frac{1}{2}B\sum_{\rm
n.n}(\vec{S_{i}}\cdot\vec{S_{j}})^{2}\\\ +\frac{1}{2}D\sum_{\rm
n.n}\hat{D}_{ij}\cdot(\vec{S_{i}}\times\vec{S_{j}})+\frac{1}{2}\Delta\sum_{i}(\vec{S_{i}}\cdot\vec{d_{i}})^{2}\end{split}$
(1)

where $\vec{S_{i}}$ denotes magnetic spins, vector $\hat{D}_{ij}$ shows
direction of DM which is determined by Moriya’s rules Moriya (1960); Elhajal
_et al._ (2005) and vector $\vec{d_{i}}$ is the single-ion easy-axis direction
at each site $i$. According to Moriya’s rules, if the mirror planes join two
sites and the middle point of the opposite bond then the DM vectors should be
perpendicular to the mirror plane Moriya (1960); Elhajal  _et al._ (2005). It
should be noted that in this model $|\vec{S}|$=1. Different terms of spin
Hamiltonian express different types of exchange interactions without
(isotropic terms) and with (anisotropic terms) considering spin-orbit coupling
(SOC). Isotropic terms, $J_{ij}$ and $B$ are Heisenberg and bi-quadratic
coupling, respectively. Heisenberg interaction decreases rapidly with distance
because the NaSrMn2F7 is an insulator. Therefore, Heisenberg terms are
calculated up to the third nearest neighbor($J_{\mathrm{1}},J_{\mathrm{2}}$
and $J_{\mathrm{3a}}$). To this end, one has to regard the conventional unit-
cell (88 atoms) with $4\times 4\times 4$ optimized Monkhorst-Pack k-point
mesh. Then, $J_{\mathrm{1}},J_{\mathrm{2}}$ and $J_{\mathrm{3a}}$ values can
be obtained by calculating the total energy of different collinear spin
configurations and mapping them to the Heisenberg model (Fig. 2). There are
two possible third nearest neighbor exchange interaction in this structure,
$J_{\text{3a}}$ and $J_{\text{3b}}$Kennedy (1997); Wills  _et al._ (2006);
Sadeghi  _et al._ (2015). Generally, in the cubic pyrochlore structures
$J_{\text{3b}}\ll J_{\text{3a}}$. Therefore, $J_{\text{3b}}$ is approximated
0Kennedy (1997); Wills  _et al._ (2006); Sadeghi  _et al._ (2015).

The nearest neighbor coupling is important for $B$ and anisotropic terms such
as Dzyaloshinskii-Moriya coupling ($D$) and single-ion anisotropy strength
($\Delta$). Therefore, the primitive cell (22 atoms) with $6\times 6\times 6$
optimized Monkhorst-Pack k-point mesh must be regarded. The strength of $B$
can be determined via the total energy of noncollinear spin configurations in
the absence of spin-orbit coupling (GGA+$U$). In this case, if one starts from
the all-in/all-out spin configuration (where all spins are oriented towards
in/out center of tetrahedra) and if one rotates the spins so that
$\vec{S_{1}}+\vec{S_{2}}+\vec{S_{3}}+\vec{S_{4}}=0$, it can be shown that
Heisenberg term is degenerated Sadeghi  _et al._ (2015). The total energy
difference of different noncollinear spin configurations can be attributed to
$B$ Sadeghi  _et al._ (2015). With the uniform rotation of magnetic moments in
the absence of SOC, the isotropic terms of the spin Hamiltonian are
degenerated. By considering the spin-orbit coupling (GGA+$U$+SOC), the total
energy difference of the various noncollinear spin configurations is entirely
attributed to $D$ provided that the SIA is negligible Sadeghi  _et al._
(2015). Given that Mn+2 ends in 3d5, it is expected that SIA has a very small
contribution. To calculate this term, one has to examine two different
noncollinear spin configurations. In the first spin configuration, one
magnetic moment points outward from the center of the tetrahedra. The other
three moments point inward toward the center of the tetrahedra. For the second
spin configuration, the magnetic moments of the first configuration rotate
about 120∘ around the $Z$ axis. It can be demonstrated that for these two spin
configurations in the absence of SOC, isotropic terms do not change. In the
presence of SOC, the $D$ remains unchanged and the total energy difference
between these two spin configurations is entirely attributed to $\Delta$ Xiang
_et al._ (2011).

Table 1: Calculated Heisenberg couplings up to third neighbors and $D$ for different $U_{\text{eff}}$ parameters. Negative and positive value denotes anti-ferromagnetic and ferromagnetic exchange interaction, respectively. Also, One can see calculated $\theta_{\text{CW}}$ by MC simulations. To account the effect of random A-site occupations, one can regard $\delta J=0.13$ meV as a weak bond disorder. It should be noted that I consider 200-300 K as temperature range for calculation of $\theta_{\text{CW}}$ similar to experimental range Sanders  _et al._ (2016). $U_{\text{eff}}$ (eV) |  $J_{\text{1}}$ (meV) |  $J_{\text{2}}$ (meV) |  $J_{\text{3}}$ (meV) |  $D$(meV) |  $\theta_{\text{CW}}$(K)
---|---|---|---|---|---
0.00 |  $-$8.37 |  $-$0.06 |  $-$0.02 | 0.100 |  $-$247.88
2.58 |  $-$4.12 |  $-$0.02 |  $-$0.02 | 0.010 |  $-$107.81
3.30 |  $-$3.46 |  $-$0.01 |  $-$0.02 | 0.005 |  $-$89.19
3.58 |  $-$3.22 |  $-$0.01 |  $-$0.02 | 0.002 |  $-$79.90
4.58 |  $-$2.53 |  $-$0.01 |  $-$0.02 |  $-$0.002 |  $-$64.05
EXP |  |  |  |  |  $-$89.70 Sanders  _et al._ (2016)

Figure 3: (Color online)Obtained the band structure of NaSrMn2F7 in (Left)GGA
and (Right)GGA+SOC approaches, respectively. The Fermi energy is set to zero.
The all-in/all-out magnetic configuration (noncollinear spin texture) has been
considered for these calculations. The figures exhibit that the SOC has
ignorable effect on energy bands near the Fermi level.

The derivation of Eq. (1) has been done by the full-potential linearized
augmented plane wave (FPLAW) method as implemented in the FLEUR FLEURgroup
code. For the exchange-correlation energy, I employ the Perdew Burke Ernzerhof
parametrization of the generalized gradient approximation (GGA) Perdew  _et
al._ (1996). The GGA+$U$ approximation Cococcioni and de Gironcoli (2005) is
used to account for the on-site Coulomb interaction for $3d$ orbitals of Mn
atoms. The optimized cut-off of wave function expansion in the interstitial
region is set to $k_{\mathrm{max}}=4.2\,\ \mathrm{a.u.}^{-1}$. The muffin-tin
radius of Na, Sr, Mn and F atoms are set to 2.6, 2.8, 2.3 and 1.4 a.u.,
respectively. To increase the accuracy of the calculations, I have considered
the 3s and 3p orbitals of the Mn as semicore electrons (FLAPW+LO).

The effective $U$ parameter ($U_{\text{eff}}$=$U-J_{\text{H}}$), where $U$ and
$J_{\text{H}}$ are screened on-site Coulomb and Hund exchange interactions, is
calculated by the Density Functional Perturbation Theory (DFPT) Timrov  _et
al._ (2018) using the Quantum-Espresso (QE) code (employing the plane waves
basis set) Giannozzi  _et al._ (2009, 2017). The optimized 40 Ry and 400 Ry
cutoff have been considered for expanding wavefunction and charge density in
plane wave, respectively. The approximation electron-ion interactions have
been controlled through GBRV ultra-soft pseudo-potential Garrity  _et al._
(2014). In order to investigate the microscopic strength of $\delta J$ within
DFT, one considers a conventional unitcell with 88 atoms. For this unitcell,
one can define 97 different configurations with different A-site distributions
by using the supercell program . Each of these distributions has a specific
and unique symmetry. For calculations of the total energy and $J$ parameter
for each configuration, QE code has been used in GGA+$U$ approach.
To capture the critical properties and low-temperature spacial spin
arrangement of NaSrMn2F7, one performs MC simulation. The effect of charge
disorder has been regarded in this simulation. Since SG materials have
multiple local minimums in their energy landscape, reaching the global minimum
especially at low-temperature is not possible. To tackle this problem, the
replica exchange method Hukushima and Nemoto (1996) should be employed as
implemented in the ESpinS package Rezaei  _et al._ (2022). The Gaussian
distribution has been considered the particular choice of disorder. In the
replica exchange algorithm, the high-temperature configurations are swapped
for the low-temperature configurations so that the system reaches thermal
equilibrium. The condition for the system to reach equilibrium is that the MC
steps should be sufficiently selected so that, first, the temperature behavior
of the specific heat and energy of the system is smooth and, second, the
position of the specific heat peak changes slightly by changing the seed and
third, the acceptance ratio of MC steps at each temperature reaches at least
to 65 percent. The different linear size of the simulation cell (L$=$8, 10,
12) is used to ensure the presence of phase transition. In this way, there are
N$\times\mathrm{L}^{3}$ spins in the three-dimensional lattice where N is the
number of spins (N$=$4 for one tetrahedron). For MC calculations, when $\delta
J$ is considered, the calculations are performed for ten different starting
points (different seeds) and then the final result is expressed as the average
of them. 2$\times 10^{6}$ MC steps per spin at each temperature are performed
for the thermal equilibrium and data collection, respectively. To reduce the
correlation between the successive data, measurements are done after skipping
every 10 MC steps.

Figure 4: (Color online) PDOS calculated for AFM-ordered NaSrMn2F7 using the
GGA functional. It illustrates the local 3$d$ DOS of a single Mn atom for
majority spins. The obtained Lowdin Charge analysis indicates that the
majority spin levels are fully occupied and minority spin levels are
negligibly occupied. That’s why the 3$d$ DOS for majority spins has not been
exhibited. It should be noted that this calculation has been done by using QE
code Giannozzi  _et al._ (2009, 2017).

Figure 5: (Color online) $\bf(Left)$ The calculated magnetic specific heat for
different $\delta J$ values. The freezing temperature is similar to the
experimental value for $\delta J$=0.34 meV. A broad anomaly can be seen just
for $\delta J$=0.34 meV which suggests the absence of long-range magnetic
ordering. For $\delta J$=0.34 meV, the red, green and black color denotes
L=12, L=10 and L=8, respectively. Changing the peak of specific heat for
different sizes of the simulation cell (L=8, 10, 12) indicates finite-size
scaling which supports the value of T${}_{\text{f}}$. The temperature of the
peak has been considered as T${}_{\text{f}}$. $\bf(right)$ The data collapsing
of different lattice sizes of specific heat by using finite-size scaling
theory to obtain T${}_{\text{f}}$ and relevant critical exponents.

##  III Results and Discussion

The quantitative derivation spin Hamiltonian (Eq. (1)) is effective to explain
or predict the magnetic properties of this pyrochlore. This information can
answer the fundamental question about how a novel ground state emerges. The
results of calculations of different spin Hamiltonian terms for different
Hubbard parameters are given in Table. LABEL:tab1. Although the obtained band
structure (Fig. 3) suggests that the material is an insulator in the GGA
approach, one needs to describe accurately the electron-electron Coulomb
correlation due to the 3d orbital of the Mn atom. Therefore, one should apply
the Hubbard potential to 3d electrons as an on-site Coulomb interaction. Based
on the ab initio calculations (DFPT method), the effective Hubbard parameter
$U_{\text{eff}}$ obtains 4.58 eV for this system.

Figure 6: (Color online) $\bf(Left)$ Ab-initio calculations of total energy
for different A-site distributions in GGA+$U$ approach. The purple circle
which has minimum total energy is related to homogeneous distribution in which
Na and Sr occupy A-site tetrahedral with two by two ratios. The total energy
of homogeneous distribution has been considered as E${}_{\text{min}}$.
$\bf(right)$ Ab-initio calculations of $J_{1}$ for different A-site
distributions in GGA+$U$ approach. According to definition of standard
deviation, one can obtain $\delta J$=0.31 meV.

In this paper, I have considered $J_{\text{H}}$$\sim$1.0 eV as it is suggested
for 3d materials Coey  _et al._ (2013); Vaugier  _et al._ (2012). To ensure
the exact value of $U_{\text{eff}}$, $\theta_{\text{CW}}$ obtained from the MC
simulation can be compared with the experimental value. Since the
$\theta_{\text{CW}}$ is strongly dependent on $J_{1}$ and $J_{1}$ is strongly
dependent to $U_{\text{eff}}$, one can control $\theta_{\text{CW}}$ value by
increasing or decreasing the Hubbard parameter. As shown in Table. LABEL:tab1,
I have started from $U_{\text{eff}}$=4.58 eV and decreased it until the
obtained $\theta_{\text{CW}}$ is equal to the experimental value. In this way,
the optimum $U_{\text{eff}}$ is 3.30 eV.

According to the chemical valence of different atoms in NaSrMn2F7 compound,
one expects that Na, Sr, Mn and F are in the +1, +2, +2 and -1 oxidation
states, respectively. Therefore, the external shell of Mn+2 is 3d5 and
Hund$\textquotesingle$s rule predicts a spin S=5/2, and thereby a magnetic
moment 5 $\mu_{B}$ for it. The DFT calculations obtain 4.8 $\mu_{B}$, not far
from the native ionic picture. With this simple picture, one expects the
quantum angular momentum of Mn to be zero and the magnetization of this
compound originated only from the spin part. The DFT calculations give some
residual magnetization 0.024 $\mu_{B}$, 0.0008 $\mu_{B}$ and 0.0005 $\mu_{B}$
for F, Sr and Na atoms, respectively.

In order to ensure that the bond disorder does not affect the crystal field
and its consequence on higher order exchange interactions, one has to
calculate other important and effective exchange couplings. The DFT
calculations give those other exchange parameters such as $B$, $D$ and
$\Delta$ are 3 to 4 times smaller than the $J_{1}$. As a result, these
interactions are negligible and do not play a significant role in determining
the magnetic properties of the ground state. The antisymmetric exchange
interactions (DMI and SIA) are strongly dependent on the orbital moment of the
magnetic ion. Experimental measurements Sanders  _et al._ (2016) show that the
Mn orbital moment has a very small role in the effective orbital of this
material. The 0.0001 $\mu_{B}$ obtains as the orbital moment of each Mn ion
from ab initio calculations (GGA+$U$+SOC), which is a very small amount. It
can be concluded that antisymmetric exchange interactions in this material are
ignorable. In addition, Fig. 3 compares band structures for noncollinear (all-
in/all-out) magnetic configuration with and without SOC. As can be seen, the
effect of the SOC on the energy bands is very small, especially near the Fermi
level. Therefore, SOC effect in this material is very weak and as a result.
The calculations of spin-projected DOS (PDOS) and Lowdin charge analysis for
an Mn atom confirm that SOC effect is negligible (Fig. 4). According to the
Lowdin charge analysis, the total charge for majority and minority spin levels
are $d_{z^{2}}$=0.9973, $d_{xz}$=0.9905, $d_{yz}$=0.9909,
$d_{x^{2}-y^{2}}$=0.9973, $d_{xy}$=0.9927 and $d_{z^{2}}$=0.0348,
$d_{xz}$=0.0723, $d_{yz}$=0.0707, $d_{x^{2}-y^{2}}$=0.0358, $d_{xy}$=0.0723,
respectively. Therefore, the majority spin levels are fully occupied and
minority spin levels are negligibly occupied. It can be concluded that the
orbital moment is quenched in this material.

The exchange parameters of the second ($J_{2}$) and third ($J_{3}$) nearest
neighbors are very small. Thus, A-site random occupation only changes $J_{1}$.
One can estimate this weak bond disorder by using $\delta
J=\sqrt{\frac{3}{8}}k_{\text{B}}$T${}_{\text{f}}$ Saunders and Chalker (2007).
Since the experimental value of the freezing temperature is 2.5 K Sanders  _et
al._ (2016), $\delta J$ is equal to 0.13 meV.

Figure 7: (Color online) The (h0l) planes density plot of elastic neutron
scattering structure function obtained by MC simulation at T=0.5 K. The change
in the scattering pattern indicates change in magnitude and intensity of the
exchange interaction which is necessary to form the SG state.

In the following, I perform MC simulations to obtain ab initio freezing
temperature by using Eq. 1 for optimum $U_{\text{eff}}$. In order to account
for the effect of A-site random occupation on MC simulation, the calculations
are performed taking into account $\delta J$. Fig. 5 illustrates the
temperature behavior of specific heat of the sample. The peak of specific heat
has been regarded as T${}_{\text{f}}$. For $\delta J$=0.13 meV,
T${}_{\text{f}}$ obtains 1.17 K. One increases the $\delta J$ to reach the
experimental T${}_{\text{f}}$. The calculations show that the ab initio
freezing temperature can be reached at 2.5 K for $\delta J$=0.34 meV. By
decreasing the amount of $\delta J$, the freezing temperature decreases so
that for $\delta J$=0, no phase transition is observed up to 0.5 K. In order
to show how much the freezing temperature (T${}_{\text{f}}$) depends on the
strength of the $\delta J$, the MC simulations are repeated for the case where
the magnitude of the $\delta
J=\sqrt{\frac{3}{8}}k_{\text{B}}$T${}_{\text{f}}$=0.13 meV is kept constant
and the value of $J_{1}$=$-$3.46 meV is increased. The results of my
calculations exhibit that the value of T${}_{\text{f}}$ has very tiny changes
from 1.17 K even for $J_{1}$ =$-$4.58 meV and it reaches $\sim$ 2 K for
$J_{1}$=$-$50 meV. Therefore, one can say that the $\delta J$ acts as a small
perturbation and must account for the system SG phase transition at low-
temperature by eliminating the degeneracy of the ground state Saunders and
Chalker (2007). One of the remarkable features of SG materials is a broad
anomaly in the specific heat data around freezing temperature. Such behavior
is observed for $\delta J$=0.34 meV.

Another question comes from how one can define the T${}_{\text{f}}$. To answer
this question, one uses the finite-size scaling theoryFisher  _et al._ (1972);
Landau (1976). The singular part of the specific heat is given by

$\displaystyle c(t,\rm L)$ $\displaystyle\sim$ $\displaystyle\rm
L^{\alpha/\nu}{\cal C}(tL^{1/\nu})$ (2)

where $t=\frac{\text{T}_{\text{f}}-\text{T}}{\text{T}_{\text{f}}}$ is the
reduced temperature and $\rm L$ is linear lattice size. Finite-size scaling of
$\alpha$ is shown in Fig. 5. This figure clearly shows the data collapse of
the different lattice sizes (L$=8,10,12$) by choosing
$\alpha$=$-$1.21$\pm$0.06, $\nu$=0.92$\pm$0.03 and
T${}_{\text{f}}$=2.6$\pm$0.05 K. This result shows that the estimation of
T${}_{\text{f}}$ by changing $\delta J$ has good accuracy. The negative value
for $\alpha$ has been reported for other SG and SG-frustrated materials
Gunnarsson  _et al._ (1991); Yeshurun  _et al._ (1980); Gingras  _et al._
(1997); Andreanov  _et al._ (2010); Hawkins  _et al._ (1979).

To prove the obtained $\delta J$=0.34 meV, one calculates the microscopic
strength of $\delta J$ within DFT. Since Na and Sr occupy A-site fully
randomly with equal parts, one defines different 97 geometries with different
Na/Sr distributions. Fig. 6 shows how different 97 configurations are arranged
in terms of total energy. According to Fig. 6, one distribution has the lowest
total energy. It can be called homogeneous distribution because Na-Sr occupies
A-site in such a way that in each A-site tetrahedral, there are two Sr and two
Na. Therefore, one can only make a primitive cell with 22 atoms for only this
distribution.

To determine the magnetic order of the system below T${}_{\text{f}}$, one
calculates the elastic neutron scattering structure function (S($\vec{q}$))
for different seeds. S($\vec{q}$) is given by Lacroix  _et al._ (2011):

$\begin{split}S(\vec{q})=&\sum_{i;j}\langle(\vec{S}_{i}-\frac{\vec{S}_{i}\bm{\cdot}\vec{q}}{\vec{q}\bm{\cdot}\vec{q}}\vec{q})\bm{\cdot}(\vec{S}_{j}-\frac{\vec{S}_{j}\bm{\cdot}\vec{q}}{\vec{q}\bm{\cdot}\vec{q}}\vec{q})\rangle\\\
&\exp[i\vec{q}\bm{\cdot}(\vec{R}_{i}-\vec{R}_{j})].\end{split}$ (3)

in which $\vec{R}_{i}$ is the position of each of the four spins in the
tetrahedron unit cell i = 1 · · · N${}_{\text{cell}}$, $\vec{q}$ is the wave
vector and $\vec{S}_{i}$ denotes the direction of each spin at the
tetrahedron. Fig. 7 show the density plots of S($\vec{q}$) for different seeds
at T=0.5 K in (h0l) planes. According to Fig. 7, S($\vec{q}$) has different
patterns. Therefore, one concludes that the ground state below $T_{\text{f}}$
has an SG order. Because the different scattering patterns show that the
distribution of $J$ differs in terms of strength and magnitude. Random
exchange interaction is one of the requirements for the emergence of the glass
spin ground state.

Another feature of SG materials is that they undergo a second-order phase
transition. For MC simulations, the order of transition is obtained by the
Binder forth energy cumulant given by Binder (1981):

$U(\rm T)=1-\frac{1}{3}\frac{\langle E^{4}\rangle}{\langle E^{2}\rangle^{2}},$
(4)

in which E is the total energy of the system which can be calculated by MC.

Figure 8: (Color online) The scaling behavior of the minimum of Binder forth
energy cumulant $U_{\text{min}}$ versus L-3. The dashed lines are obtained by
fitting Eq. (5) to the data. We use the MC simulation data on the lattices
with L$=$ 4, 6, 8, 10, 12 and 14.

The $\rm L$ dependence of minimum of $U(\rm T)$ at critical temperature is
given by Binder (1981):

$U_{\rm min}(\rm L)=U^{*}+A\rm L^{-3}+B\rm L^{-6}+O(\rm L^{-9})$ (5)

Fig. 8 illustrates the scaling of the Binder cumulant minima versus $\rm
L^{-3}$. For $U^{*}<\frac{2}{3}$ and $U^{*}$=$\frac{2}{3}$, the order of phase
transition is first and second-order, respectivelyBinder (1981). The obtained
$U^{*}$=0.666663(6) shows that this system goes to the low-temperature spin-
glass ground state with second-order phase transition.

##  IV Conclusion

In summary, a spin Hamiltonian for spin-glass pyrochlore NaSrMn2F7 by using
DFT calculations obtained. Spin Hamilton included isotropic and anisotropic
terms such as Heisenberg coupling up to third neighbors, nearest neighbor bi-
quadratic, Dzyaloshinskii-Moriya and single-ion anisotropy. Since the Hubbard
parameter played a decisive role in describing the electron-electron
correlation and consequently the magnetic properties, this value was
determined via the DFPT method. Then, I reduced it to such an extent that the
$\theta_{\mathrm{CW}}$ obtained from the MC simulations equaled to the
experimental $\theta_{\mathrm{CW}}$. The DFT calculations showed that only
$J_{1}$ (as nearest neighbor Heisenberg exchange) is the most important factor
in determining the magnetic properties of this material. In the following, the
specific heat temperature behavior was calculated by using the obtained spin
Hamiltonian via MC simulations. Since Na-Sr random occupancy on A-site causes
weak bond disorder($\delta J$), $\delta J$ was increased to such an extent
that the calculated freezing temperature equals the experimental freezing
temperature. To prove this claim, the microscopic calculation of $\delta J$
within DFT was done. To this end, $J_{1}$ distribution was calculated for
different Na/Sr distributions. The obtained result showed that $\delta J$ from
DFT calculations was in good agreement with calculated $\delta J$ by using the
Monte Carlo simulations. I found that $\delta J$ behaves like a perturbation
so that it is responsible for spin-glass phase transition at low-temperatures.
To prove the ground state below critical temperature is spin-glass, I
calculated the elastic neutron scattering structure function (S($\vec{q}$)) in
(h,0,l) planes for different seeds. I conclude that the distribution of $J$
differs in terms of strength and magnitude. By using finite-size scaling
theory, I found that the $\alpha$ as a specific heat critical exponent was
negative which is the feature of spin-glass systems. As a suggestion,
experimental measurements such as inelastic neutron scattering can help us to
understand the nature of fluctuations at low-temperatures and their effect on
electron-electron correlation.

###### Acknowledgements.

Jülich Supercomputing Centre (JSC) is gratefully acknowledged by authors as
the provider of needed computing facilities. The guidelines provided by Peter
Kratzer, Marjana Ležaić and Mojaba Alaei are highly valued and appreciated.

## References

  * Mydosh (1993) J. A. Mydosh,  _Spin glasses: an experimental introduction_ (CRC Press, 1993).
  * Binder and Young (1986) K. Binder and A. P. Young, Rev. Mod. Phys. 58, 801 (1986).
  * Taniguchi  _et al._ (2009) T. Taniguchi, T. Munenaka, and H. Sato, in  _Journal of Physics: Conference Series_, Vol. 145 (IOP Publishing, 2009) p. 012017.
  * Shintani and Tanaka (2006) H. Shintani and H. Tanaka, Nature Physics 2, 200 (2006).
  * Sanders  _et al._ (2016) M. Sanders, J. Krizan, K. Plumb, T. McQueen, and R. Cava, Journal of Physics: Condensed Matter 29, 045801 (2016).
  * Kawamura and Taniguchi (2015) H. Kawamura and T. Taniguchi,  _Chapter 1 - Spin Glasses_, edited by K. Buschow, Handbook of Magnetic Materials, Vol. 24 (Elsevier, 2015) pp. 1 – 137\.
  * Vincent and Dupuis (2018) E. Vincent and V. Dupuis, in  _Frustrated Materials and Ferroic Glasses_ (Springer, 2018) pp. 31–56.
  * Mitsumoto  _et al._ (2020) K. Mitsumoto, C. Hotta, and H. Yoshino, Phys. Rev. Lett. 124, 087201 (2020).
  * Shinaoka  _et al._ (2011) H. Shinaoka, Y. Tomita, and Y. Motome, Phys. Rev. Lett. 107, 047204 (2011).
  * Subramanian  _et al._ (1983) M. Subramanian, G. Aravamudan, and G. Subba Rao, Progress in Solid State Chemistry 15, 55 (1983).
  * Gardner  _et al._ (2010) J. S. Gardner, M. J. P. Gingras, and J. E. Greedan, Rev. Mod. Phys. 82, 53 (2010).
  * Lacroix  _et al._ (2011) C. Lacroix, P. Mendels, and F. Mila,  _Introduction to frustrated magnetism: materials, experiments, theory_ , Vol. 164 (Springer Science & Business Media, 2011).
  * Reig-i Plessis and Hallas (2021) D. Reig-i Plessis and A. M. Hallas, Phys. Rev. Materials 5, 030301 (2021).
  * Harris  _et al._ (1997) M. J. Harris, S. T. Bramwell, D. F. McMorrow, T. Zeiske, and K. W. Godfrey, Phys. Rev. Lett. 79, 2554 (1997).
  * Slobinsky  _et al._ (2021) D. Slobinsky, L. Pili, G. Baglietto, S. Grigera, and R. Borzi, Communications Physics 4, 1 (2021).
  * Moessner and Sondhi (2003) R. Moessner and S. L. Sondhi, Physical Review B 68, 064411 (2003).
  * Marlton  _et al._ (2021) F. P. Marlton, Z. Zhang, Y. Zhang, T. E. Proffen, C. D. Ling, and B. J. Kennedy, Chemistry of Materials (2021), 10.1021/acs.chemmater.0c04515.
  * Yahne  _et al._ (2021) D. Yahne, D. Pereira, L. Jaubert, L. Sanjeewa, M. Powell, J. Kolis, G. Xu, M. Enjalran, M. Gingras, and K. Ross, arXiv preprint arXiv:2101.08361 (2021).
  * Shi  _et al._ (2021) Y. Shi, C. Nisoli, and G.-W. Chern, Applied Physics Letters 118, 122407 (2021).
  * Fennell  _et al._ (2009) T. Fennell, P. Deen, A. Wildes, K. Schmalzl, D. Prabhakaran, A. Boothroyd, R. Aldus, D. McMorrow, and S. Bramwell, Science 326, 415 (2009).
  * Bramwell and Gingras (2001) S. T. Bramwell and M. J. Gingras, Science 294, 1495 (2001).
  * Krizan and Cava (2014) J. W. Krizan and R. J. Cava, Phys. Rev. B 89, 214401 (2014).
  * Krizan and Cava (2015a) J. W. Krizan and R. J. Cava, Phys. Rev. B 92, 014406 (2015a).
  * Krizan and Cava (2015b) J. Krizan and R. J. Cava, Journal of Physics: Condensed Matter 27, 296002 (2015b).
  * Andreanov  _et al._ (2010) A. Andreanov, J. T. Chalker, T. E. Saunders, and D. Sherrington, Phys. Rev. B 81, 014406 (2010).
  * Alexandradinata  _et al._ (2020) A. Alexandradinata, N. Armitage, A. Baydin, W. Bi, Y. Cao, H. J. Changlani, E. Chertkov, E. H. Neto, L. Delacretaz, I. E. Baggari, _et al._ , arXiv preprint arXiv:2010.00584 (2020).
  * Moessner and Chalker (1998) R. Moessner and J. T. Chalker, Phys. Rev. Lett. 80, 2929 (1998).
  * Yang  _et al._ (2015) J. Yang, A. Samarakoon, S. Dissanayake, H. Ueda, I. Klich, K. Iida, D. Pajerowski, N. P. Butch, Q. Huang, J. R. Copley, _et al._ , Proceedings of the National Academy of Sciences 112, 11519 (2015).
  * Cepas and Canals (2012) O. Cepas and B. Canals, Physical Review B 86, 024434 (2012).
  * Moriya (1960) T. Moriya, Phys. Rev. 120, 91 (1960).
  * Elhajal  _et al._ (2005) M. Elhajal, B. Canals, R. Sunyer, and C. Lacroix, Phys. Rev. B 71, 094420 (2005).
  * Kennedy (1997) B. J. Kennedy, Physica B: Condensed Matter 241-243, 303 (1997), proceedings of the International Conference on Neutron Scattering.
  * Wills  _et al._ (2006) A. S. Wills, M. E. Zhitomirsky, B. Canals, J. P. Sanchez, P. Bonville, P. D. de Réotier, and A. Yaouanc, Journal of Physics: Condensed Matter 18, L37 (2006).
  * Sadeghi  _et al._ (2015) A. Sadeghi, M. Alaei, F. Shahbazi, and M. J. P. Gingras, Phys. Rev. B 91, 140407 (2015).
  * Xiang  _et al._ (2011) H. J. Xiang, E. J. Kan, M.-H. Whangbo, C. Lee, S.-H. Wei, and X. G. Gong, Phys. Rev. B 83, 174402 (2011).
  * (36) FLEURgroup, “http://www.flapw.de/,” .
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Cococcioni and de Gironcoli (2005) M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71, 035105 (2005).
  * Timrov  _et al._ (2018) I. Timrov, N. Marzari, and M. Cococcioni, Phys. Rev. B 98, 085127 (2018).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ (2017) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Garrity  _et al._ (2014) K. F. Garrity, J. W. Bennett, K. M. Rabe, and D. Vanderbilt, Computational Materials Science 81, 446 (2014).
  * Hukushima and Nemoto (1996) K. Hukushima and K. Nemoto, Journal of the Physical Society of Japan 65, 1604 (1996).
  * Rezaei  _et al._ (2022) N. Rezaei, M. Alaei, and H. Akbarzadeh, Computational Materials Science 202, 110947 (2022).
  * Coey  _et al._ (2013) J. M. D. Coey, M. Venkatesan, and H. Xu, “Introduction to magnetic oxides,” in  _Functional Metal Oxides_ (John Wiley and Sons, Ltd, 2013) Chap. 1, pp. 1–49.
  * Vaugier  _et al._ (2012) L. Vaugier, H. Jiang, and S. Biermann, Phys. Rev. B 86, 165105 (2012).
  * Saunders and Chalker (2007) T. Saunders and J. Chalker, Physical review letters 98, 157201 (2007).
  * Fisher  _et al._ (1972) M. E. Fisher and M. N. Barber, Physical review letters 28, 1516 (1972).
  * Landau (1976) D. P. Landau, Phys. Rev. B 13, 2997 (1976).
  * Gunnarsson  _et al._ (1991) K. Gunnarsson, P. Svedlindh, P. Nordblad, L. Lundgren, H. Aruga, and A. Ito, Phys. Rev. B 43, 8199 (1991).
  * Yeshurun  _et al._ (1980) Y. Yeshurun, M. B. Salamon, K. V. Rao, and H. S. Chen, Phys. Rev. Lett. 45, 1366 (1980).
  * Gingras  _et al._ (1997) M. J. P. Gingras, C. V. Stager, N. P. Raju, B. D. Gaulin, and J. E. Greedan, Phys. Rev. Lett. 78, 947 (1997).
  * Andreanov  _et al._ (2010) A. Andreanov, J. T. Chalker, T. E. Saunders, and D. Sherrington, Phys. Rev. B 81, 014406 (2010).
  * Hawkins  _et al._ (1979) G. F. Hawkins, R. L. Thomas, and A. M. de Graafa), Journal of Applied Physics 50, 1709 (1979).
  * Binder (1981) K. Binder, Phys. Rev. Lett. 47, 693 (1981).


# Simulated analogues II: a new methodology for non-parametric matching of
models to observations

Rami Al-Belmpeisi,1 Vito Tuhtan,1 Mikkel Bregning Christensen,1 Rajika
Kuruwita,1,2 and Troels Haugbølle1
1Niels Bohr Institute, University of Copenhagen, Øster Voldgade 5, DK-1350
Copenhagen, Denmark
2Heidelberg Institute for Theoretical Studies, Schloß-Wolfsbrunnenweg 35,
69118 Heidelberg, Germany
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Star formation is a multi-scale problem, and only global simulations that
account for the connection from the molecular cloud scale gas flow to the
accreting protostar can reflect the observed complexity of protostellar
systems. Star-forming regions are characterised by supersonic turbulence and
as a result, it is not possible to simultaneously design models that account
for the larger environment and in detail reproduce observed stellar systems.
Instead, the stellar inventories can be matched statistically, and best
matches found that approximate specific observations. Observationally, a
combination of single-dish telescopes and interferometers are now able to
resolve the nearest protostellar objects on all scales from the protostellar
core to the inner $10\,\mathrm{AU}$. We present a new non-parametric
methodology which uses high-resolution simulations and post-processing methods
to match simulations and observations using deep learning. Our goal is to
perform a down-selection from large data sets of synthetic images to a ranked
list of best-matching candidates with respect to the observation. This is
particularly useful for binary and multiple stellar systems that form in
turbulent environments. The objective is to accelerate the rate at which we
can do such comparisons, remove biases from hand-picking matches, and
contribute to identifying the underlying physical processes that drive the
creation and evolution of observed protostellar systems.

###### keywords:

stars: formation – stars: protostars – methods: numerical – binaries: general

††pubyear: 2023††pagerange: Simulated analogues II: a new methodology for non-
parametric matching of models to observations–References

##  1 Introduction

With the advent of long wavelength interferometers like ALMA and NOEMA,
observations with tens of AU resolution have become feasible for the nearest
protostars. This gives detailed information about how stars are formed, and in
recent years, observations with high sensitivity have revealed how streamers,
outflow cavities, and connections to the larger scale filaments play an
important role in regulating a non-steady mass flow to young stars from their
surrounding gas reservoir (Murillo et al., 2022; Garufi et al., 2022;
Valdivia-Mena et al., 2022; Alves et al., 2020; Ginski et al., 2021). In order
to reproduce the complex gas kinematic structures of observed filaments or
large-scale inflows and outflows, only models employing global simulations
that properly account for the connection from the molecular cloud (MC) scale
to the accreting protostar can be used to understand the systems further
(Kuffmeier et al., 2017; Frimann et al., 2016; Kuffmeier et al., 2023).

Figure 1: Schematic overview of the Simulated Analogues method for non-
parametric matching of simulated systems to an observation. _Zoom-in
technique_ : Target binary systems from molecular cloud simulations are
selected for zoom-in simulations. _Simulated analogues catalogue_ : Synthetic
images are produced from the zoom-in simulations, and system properties are
derived (e.g. semi-major axis, masses, inclination etc). _Machine learning_ :
The simulated analogues catalogue is partitioned into a test and training data
set for the implemented DCNNs, and the best matches to observation are found.

Recently, the problem of matching non-parametric simulations to observation
has been addressed with a manual, labour-intensive procedure, which is not
scalable to larger data sets (Jørgensen et al., 2022). The interpretation of
the observation was facilitated by finding simulation frames that reproduce
similar structures to the observed line-of-sight moment-1 velocity maps.
Synthetic moment-1 maps were constructed and outputs and observer viewpoints
were identified that reproduced the observations of the binary system IRAS-2A.
The analysis was instrumental in understanding the kinematics of IRAS-2A and
framing the observation in a temporal context. It also highlighted the complex
history of protobinary systems.

Computing advances have allowed the era of machine learning with a significant
impact on data-driven discovery in astrophysics. Deep Convolutional Neural
Networks (DCNNs) have found several applications in analysing rich
astrophysical data sets. Examples include inferring Warm Dark Matter particle
masses from image field data obtained from N-body cosmological simulations
(Rose et al., 2023), determining cosmological parameters from 3D light-cone
data from dark matter halos (Hwang et al., 2023), morphologically classifying
radio galaxies (Brand et al., 2023), and detecting pulsars from numerous
candidate sources in the X-ray (Wang et al., 2021) and radio wave regimes (Guo
et al., 2019). In the field of star formation, machine learning has been
utilised to identify and characterise protostellar outflows (Xu et al., 2022,
2020b) and signatures of stellar feedback bubbles in molecular line spectra
(Xu et al., 2020a) observed from CO emission in nearby star forming regions.

In this paper we introduce the Simulated Analogues toolbox, with which we aim
to automate the process of matching simulations with observations. We use
pipelined data processing including creation of synthetic observations, and
employ DCNNs for down-selection. We use the analysis and synthetic
observations tools and pipelines introduced in Tuhtan et al. (2023) (Hereafter
referred to as Paper I). In Section 2, we provide an overview of the Simulated
Analogues Toolbox. We briefly describe the MHD simulations and production of
synthetic images in Section 2.1 and Section 2.2 respectively. In Section 2.3
our data, the simulated analogues catalogue, are discussed, while in Section
2.4, we explain in detail how Deep Learning is applied to finding matches
between simulation as observation. In Section 3 we present the results
obtained from applying Deep Learning models to observations of IRAS-2A.
Specifically, in Section 3.1 and 3.2, we apply pre-trained models to image
data from molecular and dust-continuum emission respectively. Later in Section
3.3, we train an ensemble of neural networks and employ them once again for
the analysis of the dust-continuum emission. Finally, in Section 4 and 5, we
assess the effectiveness of the method, and its application on future data.

##  2 Methods and data

In our pipeline, zoom-in simulations (Kuffmeier et al., 2017) are used to
create synthetic observations of young stars. In Figure 1 a flowchart for the
method is shown, summarising the steps. The first step involves running zoom-
in simulations on target binary star systems. The second step concerns
producing synthetic observations and deriving physical properties to build the
‘Simulated Analogues catalogue’. This catalogue then serves for training and
testing machine learning algorithms to match the simulations with the
observations. The first two steps, the zoom-in model and building the
catalogue, are described in Paper I. Here we will provide a brief overview in
Section 2.1 and 2.2, while the catalogue itself is described in Section 2.3.
The remainder of this paper is dedicated to demonstrating the validity of
using Deep Learning to match simulation and observation non-parametrically.

###  2.1 Simulations - tracking star formation inside a molecular cloud

We used ‘zoom-in’ MHD simulations (Nordlund et al., 2013; Kuffmeier et al.,
2017) of three isolated binary stars, selected from a simulation of star
formation and early evolution inside a turbulent molecular cloud with a
realistic initial mass function (Haugbølle et al., 2018) and multiplicity
distribution (Kuruwita & Haugbølle, 2023). The simulations were run with the
adaptive mesh refinement, ideal MHD code RAMSES (Teyssier, 2002; Fromang et
al., 2006).

In the zoom-in models presented here, radiation and thermodynamics during the
collapse of the protostellar core are approximated using a polytropic equation
of state with varying gamma (Masunaga & Inutsuka, 2000). We allow for the
finest resolution of either $0.8\,\,\mathrm{AU}$ or $3.1\,\,\mathrm{AU}$ in a
region which is $20\,000\,\mathrm{AU}$ across. The minimum cell size outside
this region is rapidly increased to 201 $\,\mathrm{AU}$ and then gradually
increased to 403, 806, 1611, 3223, and 6446 $\,\mathrm{AU}$ outside distances
of 0.0625 pc (13.000 $\,\mathrm{AU}$), 0.125, 0.25, 0.5, and 1 pc respectively
to create an onion-like resolution structure that allows the gas and magnetic
fields to be adequately resolved when approaching the protostellar system. We
store outputs from the evolution for up to $260\,\mathrm{kyr}$ with a cadence
of $100\,\mathrm{yr}$. This way we can accurately track the early evolution of
protostellar systems from a realistic starting condition. These models are
also described in Jørgensen et al. (2022). For a more in-depth description of
the simulations additional details can be found in Paper I.

###  2.2 Post-processing and extracting observables

To compare with observations, we construct synthetic images, spectral energy
distributions and bolometric temperatures calculated for each simulation frame
using RADMC-3D (Dullemond et al., 2012) and convolved with a Gaussian to make
the resolution comparable to the observations. This could in principle be done
with any radiative transfer tool available in the community, and the effects
of instrumentation could also be treated in more detail with e.g. CASA (Bean
et al., 2022). For each snapshot of the simulations, we place the observer
around a cone chosen such that the projected separation matches the separation
of the observed protobinary, addressing the degeneracy that arises from
projecting two points in 3D onto a 2D plane (see Figure 2). To provide a
consistent frame of reference for interpreting the images, we set the north
vector of the image plane such that it places the primary star exactly above
the secondary star, just as in the observations. We sample the full 2$\pi$
angle by evenly distributing eight observers around the cone.

Figure 2: Schematic of how lines of sight are derived from a binary, such that
the projected separation equals $200\,\mathrm{AU}$. Eight lines of sight are
sampled evenly around the cone which gives the desired projected separation.
The synthetic images below illustrate how the system would appear to observers
located along the lines of sight 1 and 4 around the cone. The primary and
secondary stars are indicated in red and yellow colours, respectively. Image
Credit: Evangelia Skoteinioti Stafyla.

The simulated analogues catalogue consisting of synthetic images is expanded
with physical quantities that describe the system, such as the protostellar
masses, accretion rates, luminosities, and disk sizes. Details on extracting
observables from the simulation can be found in Paper I. This catalogue is the
data set used for both training and testing the inference using DCNNs.

###  2.3 Synthetic observations catalogue

The simulated analogues catalogue used in this paper is generated by analysing
all output from three simulated binary protostellar systems, namely B1, B2,
and B3. The parameters of the simulated protostars can be found in Table 1.
Given an output cadence of 100 years, and eight observers, the resulting
simulated analogues catalogue encompasses a total of $45\,000$ synthetic
images for each observable. These include the dust-continuum emission
discussed in Section 3.2 and nine kinematic molecular tracers discussed in
Section 3.1, giving in total $450,000$ images for the full catalogue. For each
output, we attach the following metadata: stellar ages ($t_{1}$, $t_{2}$),
masses ($M_{1}$, $M_{2}$), accretion rates ($\dot{M}_{1}$, $\dot{M}_{2}$),
luminosities ($L_{1}$, $L_{2}$), disk sizes ($d_{size}$, $d_{size_{2}}$), the
direction of the angular momentum of the circumstellar disks
$\mathbf{J}_{1},\mathbf{J}_{2}$ and disk masses ($M_{d1}$, $M_{d2}$) for the
primary and secondary stars, respectively. We also include radial profiles of
the gas computed as spherically binned mass weighted averages centred on each
star for the angular momentum ($\vec{L}_{1}(r),\vec{L}_{2}(r)$) and radial
velocity $\vec{v}_{1}(r),\vec{v}_{2}(r)$. Finally, we attach to each output
the gas mass within a 100 AU sphere from the primary star, $M_{gas,r\leq
100\,\mathrm{AU}}$.

Name | $\Delta x_{min}$ | Final Age | M${}_{\textrm{prim}}$ | M${}_{\textrm{sec}}$ | Periods
---|---|---|---|---|---
Units | AU | kyr | $M_{\odot}$ | $M_{\odot}$ |
B1 | 3.1 | 148 | 1.07 | 0.34 | 131
B2 | 3.1 | 260 | 2.93 | 1.47 | 103
B3 | 3.1 | 148 | 0.77 | 0.68 | 265
Table 1:  A short summary of the zoom-in simulations used in this work.
$\Delta x_{min}$ is the smallest cell size on the highest level of refinement.
The final age refers to that of the oldest star (usually the primary). The
masses are measured at the end of the simulation. Periods correspond to the
number of binary orbits a given system made measured from the first periastron
crossing until the end of the simulation.

###  2.4 Deep Learning

Deep Convolutional Neural Networks (DCNNs) are incredibly versatile
mathematical models that have revolutionised the field of computer vision by
significantly improving image recognition and analysis tasks. They are
designed to mimic the visual processing capabilities of the human brain,
allowing them to learn multiple levels of representation from data
automatically. The hierarchical structure of DCNNs, with multiple layers of
interconnected neurons, enables them to create progressively more abstract
representations of the image data.

When a query image is fed into a DCNN, the output of the previous layer serves
as the input for the subsequent layer, repeatedly passing it through a
combination of mathematical operations such as 2D convolution along the
image’s spatial extent. Every layer is followed by an activation function
which ensures non-linear mapping and generates the layer output. This
representation is commonly referred to as the image’s latent representation,
or feature vector, or feature map. In the initial layers of a DCNN, the layer
outputs typically consist of small-scale features which are characterised by
high spatial dimensions. As the layers go deeper into the network, the
features become progressively less spatially refined but can identify larger
and more complex patterns. The feature maps from the final layers of the
network are capable of abstracting high-level features and shapes, ultimately
such as spirals or high contrast jumps in density. While training DCNNs
requires large amounts of labelled data and computationally intensive
processes, their impressive performance justifies the investment.

Using the last layer of a neural network configuration, we generate feature
maps for all the images and compare them to the observation feature map using
element-wise distance measures to assess a quantifiable similarity measure.
This approach allows us to rank the entire list of synthetic images from least
to most similar with respect to the observation. We can then select the most
similar images to the observation by down-selecting to a list of best
candidates, which serve as a sample for exploring patterns in the simulated
observable and physical quantities and as a basis to put the observation into
a dynamical evolutionary context.

Apart from the similar image retrieval algorithm, we also utilise feature maps
as input to train a Multi-Layer Perceptron (MLP) on top of the pre-trained
DCNNs. MLP is a feed-forward type of artificial neural network where the data
flows from the input to the output layer in one direction through one or more
"hidden layers" which are fully connected (Rosenblatt, 1958). By combining a
variant of the triplet loss scheme (Schroff et al., 2015) and a scoring
function, we expand the capability of the neural network. The goal of training
the MLP neural network is to connect prominent visual features with the binary
systemic parameters to produce a model that retrieves the most similar images
in terms of systemic parameters by design.

In both approaches, the feature maps are compared to each other through
element-wise distance measures. The aim is to retrieve simulation frames that
share maximal visual similarity with the observation. We have explored the
following distance metrics:

  1. 1.

Cosine Distance

$d_{1}(\vec{u},\vec{v})=1-\frac{\vec{u}\cdot\vec{v}}{|\vec{u}||\vec{v}|}\,,$
(1)

  2. 2.

Euclidean Distance

$d_{2}(\vec{u},\vec{v})=\sum_{i}(|u_{i}-v_{i}|^{2})^{1/2}\,,$ (2)

  3. 3.

Braycurtis Distance

$d_{3}(\vec{u},\vec{v})=\sum_{i}\left|u_{i}-v_{i}\right|/\sum_{i}\left|u_{i}+v_{i}\right|\,,$
(3)

  4. 4.

Canberra Distance

$d_{4}(\vec{u},\vec{v})=\sum_{i}\frac{\left|u_{i}-v_{i}\right|}{\left|u_{i}\right|+\left|v_{i}\right|}\,,$
(4)

  5. 5.

L1 Distance 111 L1 distance is also referred to as city-block distance or
Manhattan distance.

$d_{5}(\vec{u},\vec{v})=\sum_{i}\left|u_{i}-v_{i}\right|\,,$ (5)

where $\vec{u}$ and $\vec{v}$ are the latent representations of the images
that are being compared.

In order to systematically match an observation to the members of the
simulated analogues catalogue, several DCNN were implemented using Python
supplied with the packages of Tensorflow (Abadi et al., 2015) and Keras
(Chollet et al., 2015).

####  2.4.1 Similar image retrieval using pre-trained convolutional layers

In the initially implemented approach, the used neural networks were
initialised by loading the weights from the training on the ImageNet dataset
(Deng et al., 2009). The used networks include DenseNet121 (Huang et al.,
2017), EfficientNetB7 (Tan & Le, 2019), InceptionV3 (Szegedy et al., 2016),
MobileNet (Howard et al., 2019), ResNetRS200 (He et al., 2016a), Xception
(Chollet, 2017), VGG16 and VGG19 (Simonyan & Zisserman, 2014). Despite using
pre-trained features on a different application domain, the DCNNs are expected
to be a powerful tool in retrieving a statistical sample of the most similar
images which can later prove useful in characterising several aspects of the
observed systems.

This method employs a DCNN to extract feature maps from both the simulated
images and the query image. These feature maps are then compared using the
distance metrics mentioned in Equations 1–5 to determine the visual similarity
between the two images. The similarity ranking is obtained by sorting the
feature maps according to their computed distance measures, smaller distances
for similar images and vice versa. This approach allows for the efficient
identification of the simulated images that are most similar to the query
image.

####  2.4.2 Transfer Learning

For the results of this section we used the following DCNNs: EfficientNetB7
(Tan & Le, 2019), MobileNet (Howard et al., 2019), ResNetRS200 (He et al.,
2016a), VGG16 and VGG19 (Simonyan & Zisserman, 2014). In this approach, we do
not use the pre-trained neural network models as-is, but add additional
trained layers on top them.

To effectively train a neural network it is essential to tune its
hyperparameters, which refer to model parameters that are not learned directly
from the data, but inserted manually by the user before the training. We use a
common technique to optimise our model’s hyperparameters, grid-search, which
results in a model with the Adam optimiser (Kingma & Ba, 2014), which is
configured with: $lr=10^{-4},b_{1}=10^{-3},b_{2}=0.999,\varepsilon=10^{-7}$.
In order to create neural networks specialised for astrophysical images, we
expand the neural networks by incorporating the MLP network. Aiming for
simplicity, the implemented MLP consists of the input layer, only one hidden
layer and the output layer with sizes of 256, 128 and 100 neurons
respectively. The MLP takes synthetic image data as input, derived from the
output of the used DCNN, along with metadata which contain the physical
parameters of the system. The objective of this implementation is to further
specialise the network from 2.4.1 in connection to the system’s physical
parameters in an astrophysical context. To prevent overfitting, while using
large neural networks, we employ transfer learning techniques. This involves
"freezing" the weights of the DCNN rendering them non-updatable during back
propagation, thus effectively training exclusively the MLP network. However,
the novelty of this method relies on the interplay of two components: the
_triplet loss_ and the _scoring function_.

The triplet loss scheme is an effective training method in deep learning,
where a triplet of synthetic images (1, 2, 3) with known metadata parameters
are randomly selected at each step of the training process. By default, the
first member of the triplet is chosen to be the anchor ($a=1$). The two
remaining analogues (2, 3) are compared with respect to the anchor based on
the triplet metadata parameters. To achieve this, a scoring function,
$g(q_{i},q_{j})$, is defined to measure the similarity between the metadata
parameters
$\vec{q}=({M_{1},M_{2},L_{1},L_{2},d_{size},d_{size_{2}},M_{gas,r\leq
100\,\mathrm{AU}}})$ of a pair of simulated analogues ($i,j$). The scoring
function is constructed to satisfy the specific constraints:

  1. 1.

if two images are similar in terms of metadata parameters, the scoring
function should be small, close to zero;

  2. 2.

if two images differ significantly in terms of metadata parameters, the
scoring function should be large;

  3. 3.

The scoring function should be non-negative.

The simulated analogue with the lowest scoring function is labelled as
positive (p), while the analogue with the highest scoring function is labelled
as negative (n). The suggested scoring function is calculated element-wise for
each component of metadata $q$ by dividing the absolute difference of the
metadata parameters by the sum of the metadata parameters:

$g(q_{a},q_{j})=\frac{|q_{a}-q_{j}|}{q_{a}+q_{j}}\,,$ (6)

where the index $j$ takes values from the two remaining members of the triplet
set of images (2,3).

A score which takes into consideration the scoring functions is then used to
determine which image is considered positive or negative in the training
process. The score is defined as the difference in scoring functions of the
two competing members of the triplet (2,3) with respect to the anchor (a):

$Score=\sum_{i=1}^{k}(g(q_{a}^{i},q_{2}^{i})-g(q_{a}^{i},q_{3}^{i}))\,,$ (7)

where the summation over $i$ is performed over all the metadata parameters to
get the final score. As earlier mentioned, since our metadata consists of 6
systemic parameters for the description of a binary system, $k=6$. Given a
triplet of synthetic images along with their metadata, the score calculation
proceeds, which sets the positive (p) and negative (n) analogue according to
the following relation:

$(p,n)=\left\\{\begin{array}[]{ll}(2,3)&Score>0\\\ (3,2)&Score<0\\\
\text{Undefined, resampling triplet}&Score=0\end{array}\right.$

The mathematical expression for the triplet loss algorithm is usually
expressed in terms of the following quantity:

$\left[\left\|\vec{f^{a}}-\vec{f^{p}}\right\|-\left\|\vec{f^{a}}-\vec{f^{n}}\right\|\right]_{+}=\left[\left\|d_{ap}-d_{an}\right\|\right]_{+}\,,$
(8)

where $\vec{f^{a}}$,$\vec{f^{p}}$ and $\vec{f^{n}}$ are the 100-dimensional
feature maps from the output of the MLP for the anchor, positive, and negative
image respectively and $d_{ap}$, and $d_{an}$ are the $L_{1}$ distances
between the anchor image and positive and negative scoring images
respectively. If the $L_{1}$ distance between the anchor and positive
$\left\|\vec{f^{a}}-\vec{f^{p}}\right\|$ is smaller than the $L_{1}$ distance
between the anchor and negative $\left\|\vec{f^{a}}-\vec{f^{n}}\right\|$, then
the difference of those terms will be negative and the total result will be
zero; no learning will happen at this training step, since for this triplet
the neural network is already predicting 100-D feature maps which are closer
for the pair (a-p) compared to the pair (a-n). However, if the opposite was
true, and the negative feature map $L_{1}$ distance is closer to the anchor
than the positive is, learning is applied in this training step according to
the following algorithm: the two distances are transformed into a probability
distribution with the Softmax function, which is defined as

$P(p)=\frac{e^{d_{ap}}}{e^{d_{ap}}+e^{d_{an}}}\,,$ (9)
$P(n)=\frac{e^{d_{an}}}{e^{d_{ap}}+e^{d_{an}}}\,.$ (10)

After the probabilities are computed, the probability $P(n)$ of having a
larger distance with respect to $a$ than $p$ does is maximised (i.e.
$d_{an}\gg d_{ap}$), through the minimisation of the loss function
$Loss=P(p)+(1-P(n))$.

To accurately evaluate the model, we split the aforementioned code into eight
projections and use the ones with $\varphi=(0,\pi/2,\pi,3\pi/2)$ for training
and the projections with $\varphi=(\pi/4,3\pi/4,5\pi/4,7\pi/4)$ for testing.
We employ regularisation strategies, which are used in machine learning to
reduce overfitting and avoid our model performing poorly on unseen data. L2
regularisation includes a penalty term based on the L2 norm of the weights in
the MLP layers. Large weights are discouraged by this penalty, resulting in
more generalisable representations within the network. L2 regularisation can
be mathematically defined as:

$\text{ Loss }^{{}^{\prime}}=\text{ Loss }+\lambda\sum_{i=1}^{n}w_{i}^{2}\,,$
(11)

where $w_{i}$ are the weights of the network and $\lambda=10^{-5}$ is the
tuning parameter that sets the amount that we want to penalise large weights
within the MLP network. We use the initialisation of He et al. (2015) to set
the initial values for the weights in the MLP layers, which allows deep
networks to efficiently converge. As a last approach to address overfitting,
we use dropout (Srivastava et al., 2014), a technique that involves assigning
a probability, denoted as P, for the presence of a layer of neurons/neuron
during a training step. We implement 2D spatial dropout when passing the
feature maps from the DCNNs to the MLP. Additionally, dropout is applied to
the first two layers of the MLP. For all approaches, we use $P=0.5$.

The objective of this transfer-learning design is to develop a machine
learning model capable of assessing the similarity between observed
protostellar systems based on their images and physical parameters. To
accomplish this, we employed a variant of the triplet loss training algorithm
that maps the DCNN-produced feature maps through the MLP to generate a new
trained 100-D feature map. These 100-D latent representations are constructed
to contain a similar encoding for query images of protostellar systems with
similar physical parameters, through the non-linear mappings of a MLP network
and the incorporation of metadata. The resulting model allows for the
comparison of protostellar systems based on image and parameter similarity,
facilitating further understanding through their statistical properties.

Figure 3: Moment-1 maps close to the primary star of IRAS-2A. In each panel,
the transition lines are integrated across velocities between 2 and 12 $km$
$s^{-1}$ ($\pm 5km$ $s^{-1}$ relative to the systemic velocity) at positions
where the integrated emission exceeds 10 $\sigma$. The figure shows a
selection of transition lines for complex organic molecules and displays
consistently a quadrupole velocity structure. These data are taken from
Jørgensen et al. (2022).

##  3 Results

In this section, we use both generic pre-trained DCNNs and DCNNs specifically
trained for our purpose to identify good matches between synthetic images in
the catalogue and the high-resolution observation of the protostellar binary
IRAS-2A. We use both line and continuum tracers at different scales to match
both the small-scale kinematics in the immediate vicinity of the primary and
the larger-scale environment mapped in the continuum. Ultimately, these
matches can then be used to extract probable metadata parameters that inform
us about the age, evolutionary stage, dynamics and physical parameters of the
observed system.

Figure 4: Parameter evolution plots for system B1 (Paper I) showing
separation, mass, accretion rate, and luminosity. The vertical lines indicate
the respective time that each best match was observed using VGG16 and VGG19 on
the molecular line emission (Section 3.1), while the colour reflects the
molecule moment map that produced the best match as can be seen in the legend.
Figure 5: Distance distributions of each combination of pre-trained DCNNs and
distance metrics in Section 3.2. The distances are scaled to unit variance
($\sigma$) and zero mean. The simulated analogues in the red regions represent
matches with the smallest distances to the observations and have a distance
that is at least $2\sigma$ from the mean distance.

###  3.1 Pre-Trained DCNNs for Image Retrieval: moment-1 maps

Moment-1 maps derived from molecular tracers are a powerful tool for
understanding the kinematic structure of young protostars. Different tracers
have different critical densities, optical depths, and reaction pathways
allowing for mapping at different scales of the system. Complex kinematics and
corresponding moment-1 maps resulting from a combination of (filamentary)
inflow, outflows, and a nascent disk are observed and expected based on
modelling. This is particularly true for binary and multiple systems.
Jørgensen et al. (2022) observed the protostellar binary IRAS-2A and derived
nine moment-one maps at the 100 au scale, based on transitions in different
molecular tracers (see Figure 3).

We use the pre-trained convolutional networks on this data set to explore how
well the feature maps can be used to automatically extract similar images from
incomplete data sets covering relatively few resolution elements in a small
field of view. In general, using DCNNs on incomplete data is challenging and
has only recently been addressed (Likowski et al., 2022). However, as one can
see from the observed moment-1 maps shown in Figure 3, the emission is coming
from specific regions. To overcome this issue, we first rotate, then zoom, and
finally crop a square region from each moment map image separately such that
only a few areas in the processed image have empty values. This could also be
accomplished by including the missing data directly from the observations, if
we have had access to the raw data. Afterwards, we perform 2D nearest-
neighbour interpolation to approximate the missing pixels. As a last step, we
use cubic interpolation to scale the image to a resolution of $224^{2}$
pixels, which matches the original input size of the neural networks. Note
that the neural network can be scaled to analyse higher resolution images,
with the cost of significant up-scaling of feature map sizes, and training
time.

The cropping and rotation-specific parameters that are used for each moment-1
map are also applied to the simulation outputs, to produce synthetic moment-1
maps that capture the same scales (cropping), and orientation (rotating) as
the moment maps with respect to the binary system in the observations.

In this initial application, for simplicity, we experiment with the VGG16 and
VGG19 architectures, due to their relatively simple feed-forward architecture,
in contrast to more recent ones, such as ResNet which incorporate residual
connections between non-neighbouring layers. Residual connections are
empirically known to effectively ease the optimisation procedure during back
propagation, resulting in models that have additional shortcut identity
connections without increasing model complexity or including extra parameters
(He et al., 2016b).

After experimenting with the different distance measures we have chosen to
define the distance $D_{ij}$ between two feature vectors of the $j^{th}$
observed moment-1 map ($\vec{f}^{j}$) and the $i^{th}$ synthetic observation
($\vec{f}^{i}_{SO}$) is computed as the squared sum of normalised Cosine and
Euclidean distances

$D_{ij}=d_{1}^{2}(\vec{f}^{j},\vec{f}^{i}_{SO})/N_{ij1}+d_{2}^{2}(\vec{f}^{j},\vec{f}^{i}_{SO})/N_{ij2}\,,$
(12)

where $N_{ijk}=\max(d_{k}(\vec{f}^{j},\vec{f}^{i}_{SO}))$.

In Figure 4 we show the 50 synthetic images with the lowest distance to each
of the nine moment-1 maps as vertical lines. We see clear clustering of the
different best matches, and cross consistency, with different molecular line
maps resulting in similar best matches. The best matches are similar though
not identical to the best matches identified by hand by Jørgensen et al.
(2022), which serves as an independent validation of the approach.

###  3.2 Pre-Trained DCNNs for Image Retrieval: dust-continuum emission

To evaluate the performance of DCNNs on large-scale data sets with many
resolution elements we have applied the methodology to find simulated matches
to the dust-continuum observation of the protostellar young stellar IRAS-2A on
the thousand AU scale. Additionally, we upscaled the networks to analyse
images of $800^{2}$ resolution which results in a pixel size of
$2.5\,\mathrm{AU}$.

As for the moment-1 maps, each network ranks the synthetic images in terms of
similarity by measuring the distance between the corresponding feature vector
and the feature vector of the observed continuum flux. Figure 5 shows the
histogram distribution of distances shifted to the average value and
normalised by the standard deviation. A priori this abstract similarity
measure for the feature vectors does not necessarily correspond to a
similarity of the synthetic images to the observed image. To test this, we
evaluated the similarity of the 8 best-ranked images for each combination of
deep neural networks and distance metrics subjectively. From this, we conclude
that the Euclidean distance measure applied to the feature vectors produced by
DenseNet121 is the best combination for selecting images that are similar to
the observation, as shown in Figure 6.

Figure 6: Ranked selection of synthetic observations and the continuum image
of IRAS-2A (Jørgensen et al., 2022). IRAS-2A is shown in the middle panel of
the third row. The ranking of the images is annotated above each panel. The
ranking is determined using DenseNet121 and an Euclidean distance metric.
Figure 7: Stacked histograms and probability density functions for the
predicted system parameters using pre-trained DCNNs as described in 2.4.1. The
parameters are derived from the synthetic observations that are in the red
region in Figure 5. We use the DCNNs listed in the legend together with an
Euclidean distance metric. The overall variations of the histograms (black
line) are captured with a kernel density estimator. The predicted parameters,
comprising the age of the primary $t_{1}$, mass of the primary $M_{1}$, mass
of the secondary $M_{2}$, total mass $M$, binary separation $d$, and disk size
of the primary $d_{size}$, are arranged from top-left to bottom-right. Figure
8: Stacked histograms and probability density functions for the predicted
system parameters using the trained DCNNs which are augmented with the MLP
layers described in Section 3.3. The parameters are derived from the synthetic
observations that are at least $-2\sigma$ away from the mean distance (similar
to the red region in Figure 5). We use the DCNNs listed in the legend together
with an Euclidean distance metric. The overall variations of the histograms
(black line) are captured with a kernel density estimator. The predicted
parameters, comprising the age of the primary $t_{1}$, mass of the primary
$M_{1}$, mass of the secondary $M_{2}$, total mass $M$, binary separation $d$,
and disk size of the primary $d_{size}$, are arranged from top-left to bottom-
right.

Based on the ranking, we create a system parameter estimator tool that uses
the metadata of the best-ranked simulation outputs to construct a distribution
of plausible derived parameters for the observation. We have made a cut by
choosing all simulation outputs that have a distance in feature vectors from
the observation that is less than -2$\sigma$ (see Figure 5). The distribution
of these parameters are shown as histograms and Kernel Density Estimators
(KDE) for the binary total mass ($M$), mass of the primary ($M_{1}$), mass of
the secondary ($M_{2}$), age of the primary ($t_{1}$), separation ($d$), and
disk size of the primary ($d_{size}$) respectively in Figure 7.

Apart from an outlier generated by ResNetRS200, we see that the results
cluster around a well-defined peak corresponding to an age interval of 25 to
60 kyr. Typically, the true distance is similar to the observed distance
though for some simulation outputs, there is a high angle between the image
plane and the orbital plane. The disk size of the primary is predicted to be
somewhere between $140\,\mathrm{AU}$ and $200\,\mathrm{AU}$ and the total mass
is between 1.0 and $1.5\,\mathrm{M}_{\odot}$.

###  3.3 Trained DCNNs for Image retrieval: scoring and triplet loss

In the previous sections we have used DCNNs that are designed to distinguish
between different images but optimised for astrophysical images. Here we
introduce the technique of similar image retrieval, which includes the use of
metadata scoring and triplet loss, to train a neural network. During the
training process, each network completes 10 epochs, each of which involves
loading 18000 triplets of images in batches of 64 triplets at a time for each
weight update. The training of each network takes approximately 5 hours on a
Nvidia A100 GPU. The networks learn by minimising the loss function.

Based on the neural network ensemble predictions, we implement the system
parameter estimator tool that was previously mentioned in Section 3.2 which
produces probability density functions for the predicted system parameters, as
shown in Figure 8. These results consist of analogues within the $2\sigma$
best matches.

##  4 Discussion

Systematic matching of observations to simulations is necessary to properly
take into account degeneracies produced by projection effects and time
evolution. The method aims to determine which members of the simulated
analogues catalogue better describe a given observation. Essentially, the
problem we are solving is finding the point in time and viewing angle from
which our simulations best match the observation. “Similarity to observations”
is difficult both to define and quantify. Nevertheless, using our neural
network algorithm, we can reduce the number of possible images from tens of
thousands to a handful of simulated analogues based on either their visual
similarity alone or combining their visual similarity and systemic parameters.
Analogues such as those described in Section 3.1 and Section 3.2, are ranked
in terms of visual similarity, i.e., gas morphology. This can serve as a
useful tool for interpreting recent observations, which reveal an increasingly
complex variety of gas structures around YSO binaries (Zhang et al., 2019;
Kumar et al., 2020; Alves et al., 2020; Ginski et al., 2021). Our neural
network predictions for the age of the primary could also be used as a
constraint on the evolutionary class of each observed system, according to the
class-age relation that is described in Dunham et al. (2015). Currently, a
major bias and uncertainty is the underlying simulations catalogue. Using
zoom-in models is important for creating a catalogue with a realistic core
structure, infall history, and anchoring of magnetic fields; this is
particularly the case for binary or multiple systems. However, currently, the
available data set is limited by the number of such zoom-in models. In this
paper, we have chosen to use only three models to compare with IRAS-2A. This
could be seen as a major bias, but fortunately, these systems were already
chosen to match the overall system parameters of the observation (Jørgensen et
al., 2022) and therefore well-suited for such a comparison.

To ensure the accuracy and reliability of deep learning models within the
Simulated Analogues toolbox, it is crucial to have access to large and diverse
data sets that capture the evolution of numerous simulated systems originating
from different regions of the molecular cloud. This diversity is necessary to
capture the full range of conditions under which single, binary, and multiple
young stellar objects (YSOs) form and to enable the models to derive more
precise conclusions on the system parameters. The focus in this paper is on
binary YSOs, where the classical after-the-fact modelling of the collapse of
an isolated core is challenging at best, and the methodology is particularly
well-suited. But the method is general. We anticipate that Deep Learning
models within Simulated Analogues can be successfully extended and applied to
assess the similarity of simulations and astrophysical observations of all
kinds, as long as appropriate databases are available that contain both image
data and systemic parameters.

The IRAS-2A system observation contains some striking visual features, those
of the spiral trails which are located at the southwest and northeast. Only
some of the deep learning configurations retrieve images which also appear
similar to our eyes, while most of the combinations of distance measures and
DCNNs produce matches that do not appear to be as useful for the task of image
retrieval. Several neural networks like the pre-trained Densenet121 can down-
select to a handful of synthetic images which are similar to the observation
both at the centre and in the extended structure, like the bulge-like feature
west of the binary system, the connecting bridge between the young stellar
objects or the spiral trailing gas structures. Additionally what is
interesting to observe is that DenseNet121 can capture structures that share
similar spiral trailing features but counter-rotating. This validates that it
is a network which displays a high degree of translational, rotational, and
mirroring invariance in its features.

Focusing on small-scale features and using Deep Learning matching with VGG16
and VGG19 for the molecular line emission (see Section 3.1), the models
predict matches which cluster around distinct short periods throughout the
evolution, mainly around 38, 51 and 56 kyr, as shown by the vertical dashed
lines in Figure 4. However, most of those best matches are found in clustered
regions within the $30-60\,\mathrm{kyr}$ interval.

The results of Section 3.2 demonstrated in Figure 7 show more variety. The
distribution of prediction system parameters, such as the age of the primary,
are similar for most DCNNs, but there are some outliers. As seen in Figure 7
ResNetRS200 picks up images corresponding to a very early evolutionary state
when the two stars are still approaching each other, while EfficientNetB7
picks up best matches at very late times. The majority of the networks cluster
around the ages of 25 to 60 kyr, regardless of the exact DCNN used. The shape
of the PDF with respect to the age of the primary is reflected on the PDFs of
$M_{1}$, $M_{2}$ and $M$ since mass is increasing monotonically with time. On
the other hand, the separation between the protostars ($d$) and the disk size
($d_{size}$) follow different temporal dependencies which suggest that the
most probable 3D separation between the two stars is around $\approx
2300\,\mathrm{AU}$ and the most probable circumstellar disk size for the
primary is estimated to be $100-200\,\mathrm{AU}$.

Lastly, using a trained network, where the metadata similarity is taken into
account, Figure 8 shows that the system IRAS-2A is most likely less than
$75\,\mathrm{kyr}$ old. The individual DCNNs produce similar results, and the
best matches have a broader distribution in all system parameters. One notable
difference compare to the untrained results is that the network preferentially
find matches where the stars are closer to periastron or times where the
system has already dispensed their orbital angular momentum to the gas and
have separations smaller than $1000\,\mathrm{AU}$. While we still observe a
local maximum on the probability density distribution of the circumstellar
disk size of the primary at $150-200\,\mathrm{AU}$, there is also a strong
peak at epochs where no disk is present ($d_{size}=0$).

In conclusion, according to our current small database, most methods
qualitatively agree and suggest that the age of the primary of IRAS-2A $\sim
30-60\,\mathrm{kyr}$. However, there is variation in the predicted system
masses and separations, which depends on the neural network configuration and
training methods used. The obtained results are in good agreement with a more
observational approach and taking into account other observational indicators,
such as the class 0 classification of the system, the absence of disk
detection and the strong and misaligned outflows.

##  5 Conclusions

In this paper we have introduced a novel method to efficiently process massive
simulation data sets with the end goal of identifying good matches with
observed YSOs.

As a first pilot application, we have used machine learning and neural
networks to characterise the age and system properties of the YSO binary
IRAS-2A using simulations of binary star formation. Based on our analysis,
this YSO binary can be characterised as a young protostellar system, which is
consistent with Jørgensen et al. (2022). The predictions from our Deep
learning models, qualitatively agree that the system IRAS-2A is most likely
less than $60\,\mathrm{kyr}$.

This is a first pilot study, and the predictions for the stellar parameters
should not be over-interpreted, given that the database of simulated analogues
on which the models were trained is limited, and we are still exploring the
best practises for how to use machine learning to identify best matches. The
evolution of YSO system parameters can be significantly influenced by the
kinematic structures and gas flow present in the interstellar environment
where the simulated system is located. This suggests that by including
additional simulated binary systems, the metadata parameter space can be
expanded, as more evolutionary curves of systems are added. Expanding to tens
or hundreds of simulated systems is crucial in building a robust and versatile
database that enables machine learning models to achieve generalisation rather
than over-fitting the limited available data, and for automating the process
in conjunction with e.g. larger observational surveys such as eDisk (Ohashi et
al., 2023). While constructing models using a data set that consists of only
three systems and their corresponding evolution, we expect that these models
will perform well only in analysing observations that correspond to similar
physical conditions to our simulations.

In this paper, we have shown the viability and potential of the method, but
much still has to be done in terms of characterising the different DCNNs and
their applicability to astronomical images.

We expect that in the coming years, several generic DCNNs will be developed,
which will enable machine learning models to extract even more meaningful
latent representations from the input images, allowing a more robust image-
analysis pipeline. These can readily be incorporated into the work-flow, but
have to be quality controlled through comparison with other networks and
direct inspection of both the selected synthetic images and the discarded
images.

The ALMA archive contains a vast collection of high-resolution images and data
of young protostars. Combined with an increasing set of simulated protostellar
systems, finding simulated analogues is a very valuable tool. It is our hope
that the simulated analogues toolbox can ultimately accelerate the scientific
learning rate from observations, by speeding up the analysis, allowing us to
learn more about the dynamical state of the systems, and providing robust
interpretations of the data.

## Acknowledgements

The research leading to these results has received funding from the
Independent Research Fund Denmark through grant No. DFF 8021-00350B (TH,RLK).
This project has received funding from the European Union’s Horizon 2020
research and innovation program under the Marie Sklodowska-Curie grant
agreement No. 847523 ‘INTERACTIONS’. The astrophysics HPC facility at the
University of Copenhagen, supported by research grants from the Carlsberg,
Novo, and Villum foundations, was used for carrying out the simulations, the
analysis, and the long-term storage of the results. RLK also acknowledge
funding from the Klaus Tschira Foundation. yt (Turk et al., 2011), MATPLOTLIB
(Hunter, 2007) and NumPy (Harris et al., 2020) were used to visualise and
analyse these simulations. Variants of VGG16, VGG19 (Simonyan & Zisserman,
2014), InceptionV3 (Szegedy et al., 2016), Resnet200 (He et al., 2016a),
Xception (Chollet, 2017), DenseNet121 (Huang et al., 2017), MobileNet (Howard
et al., 2019) were used to implement our algorithms using the packages of
TensorFlow (Abadi et al., 2015) and Keras (Chollet et al., 2015). Special
thanks to Evangelia Skoteinioti Stafyla for the illustration in Figure 2.

## Data Availability

All data and tools are available from the authors at reasonable request.

## References

  * Abadi et al. (2015) Abadi M., et al., 2015, TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, https://www.tensorflow.org/
  * Alves et al. (2020) Alves F. O., Cleeves L. I., Girart J. M., Zhu Z., Franco G. A. P., Zurlo A., Caselli P., 2020, ApJL, 904, L6
  * Bean et al. (2022) Bean B., et al., 2022, Publications of the Astronomical Society of the Pacific, 134, 114501
  * Brand et al. (2023) Brand K., Grobler T. L., Kleynhans W., Vaccari M., Prescott M., Becker B., 2023, MNRAS, 522, 292
  * Chollet (2017) Chollet F., 2017, in Proceedings of the IEEE conference on computer vision and pattern recognition. pp 1251–1258
  * Chollet et al. (2015) Chollet F., et al., 2015, Keras, https://github.com/fchollet/keras
  * Deng et al. (2009) Deng J., Dong W., Socher R., Li L.-J., Li K., Fei-Fei L., 2009, in 2009 IEEE conference on computer vision and pattern recognition. pp 248–255
  * Dullemond et al. (2012) Dullemond C. P., Juhasz A., Pohl A., Sereshti F., Shetty R., Peters T., Commercon B., Flock M., 2012, Astrophysics Source Code Library, p. ascl:1202.015
  * Dunham et al. (2015) Dunham M. M., et al., 2015, ApJSS, 220, 11
  * Frimann et al. (2016) Frimann S., Jørgensen J. K., Haugbølle T., 2016, A&A, 587, A59
  * Fromang et al. (2006) Fromang S., Hennebelle P., Teyssier R., 2006, A&A, 457, 371
  * Garufi et al. (2022) Garufi A., et al., 2022, A&A, 658, A104
  * Ginski et al. (2021) Ginski C., et al., 2021, ApJL, 908, L25
  * Guo et al. (2019) Guo P., et al., 2019, MNRAS, 490, 5424
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Haugbølle et al. (2018) Haugbølle T., Padoan P., Nordlund Å., 2018, ApJ, 854, 35
  * He et al. (2015) He K., Zhang X., Ren S., Sun J., 2015, in Proceedings of the IEEE international conference on computer vision. pp 1026–1034
  * He et al. (2016a) He K., Zhang X., Ren S., Sun J., 2016a, in European conference on computer vision. pp 630–645
  * He et al. (2016b) He K., Zhang X., Ren S., Sun J., 2016b, in Proceedings of the IEEE conference on computer vision and pattern recognition. pp 770–778
  * Howard et al. (2019) Howard A., et al., 2019, in Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV).
  * Huang et al. (2017) Huang G., Liu Z., Van Der Maaten L., Weinberger K. Q., 2017, in Proceedings of the IEEE conference on computer vision and pattern recognition. pp 4700–4708
  * Hunter (2007) Hunter J. D., 2007, Computing In Science & Engineering, 9, 90
  * Hwang et al. (2023) Hwang S. Y., Sabiu C. G., Park I., Hong S. E., 2023, arXiv e-prints, p. arXiv:2304.08192
  * Jørgensen et al. (2022) Jørgensen J. K., Kuruwita R. L., Harsono D., Haugbølle T., Kristensen L. E., Bergin E. A., 2022, Nature, 606, 272
  * Kingma & Ba (2014) Kingma D. P., Ba J., 2014, arXiv e-prints, p. arXiv:1412.6980
  * Kuffmeier et al. (2017) Kuffmeier M., Haugbølle T., Nordlund Å., 2017, ApJ, 846, 7
  * Kuffmeier et al. (2023) Kuffmeier M., Jensen S. S., Haugbølle T., 2023, European Physical Journal Plus, 138, 272
  * Kumar et al. (2020) Kumar M. S. N., Palmeirim P., Arzoumanian D., Inutsuka S. I., 2020, A&A, 642, A87
  * Kuruwita & Haugbølle (2023) Kuruwita R. L., Haugbølle T., 2023, A&A, 674, A196
  * Likowski et al. (2022) Likowski M., Smieja M., Struski L., Tabor J., 2022, in 2022 IEEE/CVF Winter Conference on Applications of Computer Vision (WACV). IEEE Computer Society, Los Alamitos, CA, USA, pp 2917–2926, doi:10.1109/WACV51458.2022.00297, https://doi.ieeecomputersociety.org/10.1109/WACV51458.2022.00297
  * Masunaga & Inutsuka (2000) Masunaga H., Inutsuka S.-i., 2000, ApJ, 531, 350
  * Murillo et al. (2022) Murillo N. M., van Dishoeck E. F., Hacar A., Harsono D., Jørgensen J. K., 2022, A&A, 658, A53
  * Nordlund et al. (2013) Nordlund Å., Haugbølle T., Küffmeier M., Padoan P., Vasileiades A., 2013, Proceedings of the International Astronomical Union, 8, 131
  * Ohashi et al. (2023) Ohashi N., et al., 2023, ApJ, 951, 8
  * Rose et al. (2023) Rose J. C., et al., 2023, arXiv e-prints, p. arXiv:2304.14432
  * Rosenblatt (1958) Rosenblatt F., 1958, Psychological Review, 65, 386
  * Schroff et al. (2015) Schroff F., Kalenichenko D., Philbin J., 2015, in Proceedings of the IEEE conference on computer vision and pattern recognition. pp 815–823
  * Simonyan & Zisserman (2014) Simonyan K., Zisserman A., 2014, arXiv preprint arXiv:1409.1556
  * Srivastava et al. (2014) Srivastava N., Hinton G., Krizhevsky A., Sutskever I., Salakhutdinov R., 2014, The journal of machine learning research, 15, 1929
  * Szegedy et al. (2016) Szegedy C., Vanhoucke V., Ioffe S., Shlens J., Wojna Z., 2016, in Proceedings of the IEEE conference on computer vision and pattern recognition. pp 2818–2826
  * Tan & Le (2019) Tan M., Le Q., 2019, in International conference on machine learning. pp 6105–6114
  * Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
  * Tuhtan et al. (2023) Tuhtan V., Al-Belmpeisi R., Christensen M., Kuruwita R., Haugbølle T., 2023, MNRAS
  * Turk et al. (2011) Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2011, ApJS, 192, 9
  * Valdivia-Mena et al. (2022) Valdivia-Mena M., et al., 2022, A&A, 667, A12
  * Wang et al. (2021) Wang L.-Q., Jin J., Liu L., Shen Y., 2021, Research in A&A, 21, 171
  * Xu et al. (2020a) Xu D., Offner S. S., Gutermuth R., Van Oort C., 2020a, ApJ, 890, 64
  * Xu et al. (2020b) Xu D., Offner S. S., Gutermuth R., Van Oort C., 2020b, ApJ, 905, 172
  * Xu et al. (2022) Xu D., Offner S. S., Gutermuth R., Tan J. C., 2022, ApJ, 941, 81
  * Zhang et al. (2019) Zhang Y., et al., 2019, Nature Astronomy, 3, 517