# Physics-Guided and Physics-Explainable Recurrent Neural Network for Time
Dynamics in Optical Resonances

Yingheng Tang  1,2 Jichao Fan  1 Xinwei Li  3 Jianzhu Ma  4 Minghao Qi  2
Cunxi Yu  1,∗ and Weilu Gao1,∗

###### Abstract

Understanding the time evolution of physical systems is crucial to revealing
fundamental characteristics that are hidden in frequency domain. In optical
science, high-quality resonance cavities and enhanced interactions with
matters are at the heart of modern quantum technologies. However, capturing
their time dynamics in real-world scenarios suffers from long data acquisition
and low analysis accuracy due to slow convergence and limited time window.
Here, we report a physics-guided and physics-explainable recurrent neural
network to precisely forecast the time-domain response of resonance features
with the shortest acquired input sequence being 7% of full length, and to
infer corresponding resonance frequencies. The model is trained in a two-step
multi-fidelity framework for high-accuracy forecast, where the first step is
based on a large amount of low-fidelity physical-model-generated synthetic
data and second step involves a small set of high-fidelity application-
oriented observational data. Through both simulations and experiments, we
demonstrate that the model is universally applicable to a wide range of
resonances, including dielectric metasurfaces, graphene plasmonics, and
ultrastrongly coupled Landau polaritons, where our model accurately captures
small signal features and learns essential physical quantities. The
demonstrated machine learning algorithm offers a new way to accelerate the
exploration of physical phenomena and the design of devices under resonance-
enhanced light-matter interaction.

Department of Electrical and Computer Engineering, University of Utah, Salt
Lake City, UT 84112, USA

Department of Electrical and Computer Engineering, Purdue University, West
Lafayette, IN 47907, USA

Division of Physics, Mathematics and Astronomy, California Institute of
Technology, Pasadena, CA 91125, USA

Institute of Artificial Intelligence, Peking University, Beijing, China ∗To
whom correspondence should be addressed; E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
Forecasting time dynamics has been a key area in many contexts of scientific
research and commercial decision making, such as climate modeling[1], medical
science[2], and finance[3]. Classical parametric models informed by domain
expertise, such as autoregressive integrated moving average[4], enjoy merits
of model simplicity and easy solutions while suffer from large prediction
errors in strongly nonlinear and stochastic data. Machine leaning approaches
provide new perspective of learning time dynamics and underlying complex
representation in a data-driven manner without fixed parameters or
structures[5]. Especially, given the impressive success of recurrent neural
networks (RNNs) in natural language processing and logical interpretation of
time series as sequences, various RNN architectures, such as long short-term
memory and gated recurrent unit (GRU), have been employed for time series
forecast applications[6, 7, 8]. However, the overfitting and local minimum
problems often lead to low model efficacy[9]. Moreover, the predictive power
of these models is limited to forecasting short sequences, while the high
accuracy long sequence forecast is more desirable for practical scenarios[10].
In addition, conventional machine learning approaches generally require
significant amount of training data[11, 12] that is infeasible in practical
context and they are also unable to extract explainable information and
knowledge from the dataset. The fusion of machine learning models with
empirical models is a promising effective and efficient learning philosophy to
address these limitations[13, 14], especially physics-guided machine learning
methods in physical science[15, 16].

High-quality optical cavities for enhancing light-matter interaction form the
backbone of optical quantum technologies[17, 18, 19, 20, 21, 22, 23]. Among
them, one example showing the significance of capturing and understanding time
dynamics is the terahertz time-domain spectroscopy (THz-TDS) of THz resonance
structures, which continue to stimulate much interest in diverse disciplines
and applications such as health[24], sensing and imaging[25, 26],
security[27], computing[28, 29], and communication[30], in addition to quantum
applications[31, 32]. The uniqueness of THz-TDS includes the simultaneous
acquisition of amplitude and phase information, broad spectral coverage, and
time-resolved capability[33, 34, 35]. However, the energy decaying in
resonance features is notoriously slow, leading to long data acquisition time.
Moreover, theoretically infinitely long THz time signal is inevitable to be
practically truncated for reasonable experimental measurement and numerical
simulation time or preventing undesired phenomena. Thus, scientific and
technical trade-offs between data accuracy and acquisition time have to be
carefully crafted.

Here, we describe a physics-guided and physics-explainable cascaded GRU
networks model to forecast long-sequence time-domain signal using short-
sequence input obtained from THz-TDS numerical simulations and experiments.
Instead of directly training the model using a significant amount of high-
fidelity data from either computationally expensive electromagnetic solvers or
time consuming experiments, we employ a two-step multi-fidelity training
approach. A large number of low-fidelity physical-model-generated synthetic
data is employed to first pretrain the model, which reduces the hypothesis and
search space of the network for fast and efficient learning, alleviates local
minimum problems for high-accuracy forecast, and generalizes model
applicability. Through transfer learning using only a small set of high-
fidelity application-oriented data, the pretrained model is tailored to a
broad range of resonance features, including resonant dielectric
metasurfaces[36], electro-optic graphene plasmonics[37], and the ultrastrong
coupling between electron cyclotron resonance in a high-mobility two-
dimensional electron gas (2DEG) and photons in a high-quality-factor
cavity[32]. The cascaded networks enable precise long sequence forecast with
the shortest input time sequence being 7% of full sequence length, suggesting
the best $15\times$ data acquisition speedup. Furthermore, this model
accurately captures signal features that only occupy 0.01 % of total signal
energy in the experimentally collected data of Laudau polaritons, and
simultaneously learns resonance energies of spectral features. The polariton
dispersions obtained from experimentally measured and forecast time series, as
well as model-learned quantities, show excellent agreement, and the derived
cooperativity from forecast spectra matched experimentally reported value
well[32].

Figure 1a illustrates the model of cascaded GRU networks taking short input
sequence (length $k$) and forecasting long output sequence (length $L$). The
initial input sequence is used in the first GRU network and then combined with
forecast output sequence for the input of the next-stage network.
Theoretically, the minimum required number of GRU networks is
$O$(log2($L/k$)). In a middle stage of GRU networks, a physics-explainable GRU
decoder is branched out for simultaneously learning corresponding physical
quantities, resonance energies, associated with time signals (golden dashed
rectangle in Fig. 1a). Hidden states from the GRU encoder and the time series
data from previous stage GRU network make the input for both physics-
explainable GRU decoder and the forecast GRU generator to produce the time
series output of this GRU stage. The physics-explainable GRU decoder suggests
that hidden states obtained from time series forecast GRU networks bear actual
physical meaning, and is also practically useful in experiments for
accelerating not only the acquisition but also the explanation and
understanding of time dynamics data.

Figure 1b displays the training process of GRU networks. Conventional training
approach (orange path in Fig. 1b) requires a large training dataset, which
generally takes long time to be generated through simulations or experiments
and frequently leads to suboptimal trained networks because of local minimum
issues. In contrast, we utilize a physics-guided multi-fidelity two-step
training approach (blue path in Fig. 1b), where random networks are first
pretrained with a large number of low-fidelity synthetic data instantaneously
generated from analytical physical models based on domain expertise. The
pretrained network is then fine trained through transfer learning with a small
high-fidelity specific dataset from various applications, which is obtained
from either numerical simulations or experiments. The fully trained network is
optimal and has superior performance over the one trained through conventional
approach.

We first demonstrate our model and training methodology for time series
forecast in an example of resonant dielectric metasurfaces, which consists of
a periodic array of unit cells with four dielectric cylindrical pillars of
varying heights ($h$) and diameters ($d$)[36]; see Fig. 2a. We generated
training data by numerically simulating electrical field time response of
different structures with randomly generated combinations of $h$ and $d$;
detailed data generation procedure is described in _Methods_. In addition to
bright photonic modes, the coupling between neighboring pillars can generate
sharp Fano resonances. Figure 2b displays the training loss as a function of
training epochs using conventional and two-step approaches. Despite extensive
hyper-parameter tuning in conventional training approaches (see more details
on hyper-parameters in _Methods_), we clearly observe local minimum issue that
is prevailing in non-convex optimization problems and leads to poor forecast
(blue dashed line with cross markers).

In order to solve this problem and improve forecast accuracy, we employ a
multi-fidelity training framework[15]. Most physical resonance features and
structures render themselves as a sum of damped oscillations in time signals,
which follow a general mathematical form
$\Sigma_{i}A_{i}e^{-\alpha_{i}t}\text{sin}\\\ (\omega_{i}t+\phi_{i})$, where
$A_{i}$ is the amplitude factor, $e^{-\alpha_{i}t}$ describes harmonics
envelope decay, $\text{sin}(\omega_{i}t+\phi_{i})$ is an oscillating carrier
with a resonance frequency $\omega_{i}$ and initial phase $\phi_{i}$, and $i$
enumerates all resonance features. We employed this analytical model to
generate low-fidelity synthetic data, by randomly choosing $\omega_{i}$ within
the frequency of interest, and arbitrarily selecting $\alpha_{i}$, $\phi$ and
the number of resonances; detailed generation procedures, data pre-processing,
and exemplary synthetic data are in _Supplementary Information Sections 1, 2,
and Fig. S1_. All GRU networks were pretrained with these synthetic data and
then fine trained via transfer learning with a small set of high-fidelity
time-domain data of the dielectric metasurface obtained through full-wave
finite domain time difference (FDTD) simulations (see _Methods_ for training
data generation). The ratio of the length of final sequence over that of input
sequence, $L/k$, is $\sim 15$. As clearly indicated in Figs. 2b and c, Mean-
Squared-Error (MSE) and the accuracy of forecast time series utilizing the
multi-fidelity training approach outperforms those obtained using conventional
approach. Furthermore, in the frequency domain, the spectra that is calculated
through the Fourier transformation of forecast time signals obtained using the
model trained by the multi-fidelity approach also show excellent agreement
with the target spectra calculated by running full-time simulations.

The transfer learning technique greatly generalizes the applicability of
pretrained networks. The similar damped oscillation signature of the time-
domain signals in most resonance features physically guarantees the
feasibility of such generalization. We demonstrated the approach
generalization to two physically distinct resonance features, active graphene
plasmonics[37] and ultrastrongly coupled Laudau polaritons[31, 32], in
addition to dielectric metasurfaces. Specifically, periodically patterned
monolayer graphene ribbons, which are expected to support localized plasmonic
THz resonance from bounded carriers by ribbon boundaries, are simulated to
obtain time-domain signals; see Fig. 3a. The ribbon width and graphene Fermi
level are arbitrarily selected; details can be found in _Methods_. With
similar pretrained networks and a small set of short input sequences ($L/k\sim
4$), both time and frequency domain forecast show excellent agreement with the
target time signal (Fig. 3b) and spectra (Fig. 3c) obtained from FDTD
simulations.

Moreover, an ultrahigh-mobility 2DEG inside a high-quality-factor one-
dimensional THz photonic cavity under magnetic fields displayed ultra-narrow
Landau polaritons[32]; see Fig. 3d. Their spectra were experimentally measured
with THz-TDS under high magnetic fields. There are in total 71 spectra under
various magnetic fields from 0 T to 4.5 T, and we used 24 spectra as training
data to fine-train synthetic-data-trained model and the rest 47 spectra as
test data. These 24 spectra were chosen as ones taken under magnetic fields
either uniformly or randomly distributed between 0 T and 4.5 T. We then used
the model trained by selected 24 measurements to forecast the rest 47 time
signals under other magnetic fields with corresponding short input sequence to
reveal full dispersion of Laudau polaritons. Compared with the other two
simulation examples, experimental data are much more noisy and the number of
dataset is very limited. Moreover, our Laudau polariton features locate inside
the defect mode of the stopband of photonic crystals, and the signal energy of
interest only occupies 0.01 % of total signal energy; see Fig. S2 in
_Supplementary Information_. All these factors make the forecast using
conventional direct training approach nearly impossible in practice. The
multi-fidelity two-step approach instead significantly improves the predictive
power and accuracy, and thus small but important resonance features are
captured. The loss function is also the MSE loss calculated from predicted and
target time signals; details on the loss function selection are in
_Supplementary Information Section 3 and Fig. S3_. Figures 3e and f display
one representative time signal and corresponding spectra from Fourier
transformation both aligning well with target experimental results.

For Laudau polariton experiments, in addition to cascaded GRU networks for
time-domain signal forecast, we added a physics-explainable GRU network as
shown in Fig. 1a in the middle (2nd) GRU network to infer resonance energies
associated with time signals. The cascaded GRU networks are first trained for
time series forecast, and then the physics-explainable GRU network is trained
with the hidden states and intermediate sequences from time series forecast
GRU networks as input; see _Methods_ for detailed training process. After
training, the short time signals from test data are input into the full model,
and both forecast time signals and corresponding resonance energies are
simultaneously generated.

Figures 4a and b display all 47 experimental and predicted spectra obtained
from the Fourier transformation of corresponding time signals, where $L/k$ is
$\sim 6$ for predicted time sequences. All essential physical features
expected in linearly polarized transmission spectra are well reproduced,
including cyclotron-resonance-active lower polaritons (CRA-LPs), CRA upper
polaritons (CRA-UPs), and CR-inactive (CRI) modes[32]; see _Supplementary
Information Section 4 and Figs. S4 and S5_ for more data in both time and
frequency domains. In a stark contrast, the spectra obtained directly from
short input sequences with appropriate zero padding, as shown in Fig. 4c,
display completely random patterns with all features missing. The distinct
difference between Figs. 4a, 4b, and 4c highlights the necessity of long data
acquisition to capture essential THz features if no prediction is employed, as
well as the high predictive power and accuracy of our model of cascaded GRU
networks and training approach. Note that the results do not depend on how 24
training spectra are selected and the average validation loss displays small
variance when training spectra are randomly shuffled for different magnetic
fields; see _Fig. S6 and Supplementary Information Section 4_.

Furthermore, yellow stars in Fig. 4d, orange circles in Fig. 4e, and gold
triangles in Fig. 4f show extracted peak positions from experimental measured
spectra that are the Fourier transformation of time signals, model-learned
resonance frequencies, and extracted peaks from predicted spectra,
respectively. We utilized transfer matrix formalism to calculate transmission
spectra for samples under various magnetic fields, through which the coupling
ratio $g$ can be extracted; details can be found in _Supplementary Information
Section 5_ and Ref. [32]. Figs 4d, e, and f also display calculated
transmittance colour contour map to match extracted and learned peaks. Despite
noisy CRA-UP branch in predicted spectra, the fitting can be done uniquely.
The extracted coupling ratios, $g_{\text{exp}}$, $g_{\text{explain}}$, and
$g_{\text{pred}}$ from experimental spectra, physics-explainable GRU decoder,
and predicted spectra, are 150.1 GHz, 145.6 GHz and 147.9 GHz. This excellent
agreement confirms not only the precise forecast of cascaded GRU networks but
also generated hidden states bear actual and explainable physical meanings,
which are resonance energies.

With obtained predicted dispersion, the magnetic field corresponding to the
conditions at zero detuning and far away from zero detuning can be determined.
As shown in Figs. 4g, h, and i, we predicted spectra of CRI mode at 3 T and
the CRA-UP and CRA-LP peaks at zero detuning using short input sequence with
$L/k\sim 2$ using 24 training spectra and determined their spectral line
widths to be 4.6 GHz, 5.2 GHz and 5.2 GHz, respectively. Thus, the
cooperativity from predicted sequence is 3663, which is very close to the
experimental value 3513[32]. In addition to the acceleration of experimental
data acquisition for mapping out Laudau polariton dispersion and determining
corresponding physical quantities, our model is useful to assist and guide
experiments by dynamically forecasting time signals through active learning.
We utilized measured spectra from 4.5 T to a magnetic field close to zero
detuning (e.g., $B_{\text{m}}=1.4$ T) as training spectra to forecast the
spectrum under next lower magnetic field ($B_{\text{m-1}}$). We then augmented
the training spectra by adding the measured spectra under $B_{\text{m-1}}$ for
forecasting the spectra under $B_{\text{m-2}}$, and this process keeps moving
to the lowest magnetic field. The anomalous increase of training loss during
this active learning process indicates new time signal and spectral features,
such as additional transmission peaks in the stopband around zero detuning,
and helps experimentalists to make their decisions; see _Fig. S7 and
Supplementary Information Section 4_ for more details.

Our physics-guided and physics-explainable cascaded GRU networks combine
complementary advantages of both physics-based models and machine learning
approaches for fast, efficient, and accurate forecast of long time signals in
optical resonances with a minimum of 7 % input sequence, which corresponds to
a best $15\times$ speedup of data acquisition. The multi-fidelity two-step
training framework enables the model to be generalized to varieties of
contexts of resonance features and structures because of the similarity of
their physical appearance. This model is especially efficient in practical
experimental scenarios, where long data acquisition is inevitable and the
collection of a large number of data is almost infeasible. The incorporation
of synthetic data in training process together with transfer learning
technique significantly reduces the required number of training data to 24
spectra for forecasting experimental spectra of Laudau polaritons, and
improves the predictive power and accuracy for capturing signal features that
only occupy 0.01 % of total signal energy. This approach is promising to
accelerate the discovery of new phenomena in complex systems through the
analogy with more accessible optical system[38] and the exploration of device
functionalities under resonant light-matter interactions.

##  0.1 Training data generation for dielectric metasurface and graphene
plasmonics.

Finite domain time difference (FDTD) time-domain simulations implemented in
commercial Ansys Lumerical software are used for the training data generation
in the examples of the dielectric metasurface and graphene plasmonics. For the
dielectric metasurface example, we used four silicon cylindrical rods as the
base structure with periodic boundary condition. Different data were generated
by randomly selecting the radii of four cylindrical rods. The radius is chosen
in the range from $39.5\,\mu$m to $44.5\,\mu$m with a step resolution of
$0.25\,\mu$m. We generated a total number of $5500$ samples and using $5000$
of them as the training set and the rest $500$ as the validation set. For the
graphene plasmonics example, graphene monolayer is modeled as a 2D rectangle
conducting sheet in Lumerical material library, including both interband and
intraband contributions. Fermi level and scattering rate are two parameters
used to calculate dielectric constants used for the software. The dataset is
generated by randomly sweeping graphene ribbon width and chemical potential.
The width range is between $3.8\,\mu$m to $13.8\,\mu$m with a step of
$0.2\,\mu$m. The chemical potential ranges from $0.18$ eV to $0.41$ eV with a
step of $0.01$ eV. We generated a total number of $1000$ samples and using
$800$ of them as the training set and the rest $200$ as the validation set.

##  0.2 Training of recurrent neural networks with GRU.

The multi-fidelity two-step training contains two steps: pretrain and fine-
train (transfer learning). During the pretrain stage, the synthetic data are
divided into pieces with different sequence length and fed into the model.
During the fine-train process, the actual data (simulations and experimental
measurements) are used to further update all parameters of the model. The GRU
network we used here is a four-layer structure with a size of $200$ hidden
states. In the pretrain stage, batch size is set as $128$ and the total epoch
number is $80$. The Adam optimizer is used with the initial learning rate set
at $5\times 10^{-4}$. The learning rate decays every 30 epochs with the decay
rate 0.01. The detailed hyper-parameters are summarized in the _Supplementary
Information Section 6 Table 1_.

In the example of Landau polaritons, once the model for time dynamics forecast
is fully trained, a physics-explainable GRU decoder is connected to the
trained encoder of the 2nd GRU to receive the encoded hidden states and
processed time signal outputs. We trained this physics-explainable GRU network
with resonance frequency labels associated with time-domain signals, while we
kept weights of the rest cascaded GRU networks unchanged and only updated this
branch GRU decoder. The decoder is a four-hidden layer GRU network similar to
the decoders in time signal forecasting networks. In limited 71 experimental
spectral data, we utilized 24 data for training and 47 data for testing. This
physics-explainable network was also trained using the Adam optimizer with a
learning rate of $5\times 10^{-5}$. The total epoch for the training is 300
with the batch size of 1. The learning rate decays every 100 epochs with the
decay rate 0.1.

Upon publication, source data will be provided with this paper. All other data
that support the plots within this paper and other findings of this study will
be available from the corresponding authors upon reasonable request.

Upon publication, the codes relevant to neural network models and associated
simulation data will be publicly available on _GitHub_. The experimental data
will be available from the corresponding authors upon reasonable request. All
other codes that support the plots within this paper and other findings of
this study are available from the corresponding authors upon reasonable
request.

J. F. and W. G. thank the support from the University of Utah start-up fund.
C. Y. thanks the support from grants NSF-2019336 and NSF-2008144.

C. Y. and W. G. conceived the idea and designed the project. Y. T. performed
the modeling and calculations with the help of J. F., J. M., M. Q., C. Y. and
W. G. X. L. helped with the analysis of Landau polariton data. Y. T. and W. G.
wrote the manuscript. All authors discussed the manuscripts and provided the
feedback.

The authors declare that they have no competing financial interests.

Correspondence and requests for materials should be addressed to Cunxi Yu
(email<EMAIL_ADDRESS>and Weilu Gao (email: weilu.gao@utah.edu).

Figure 1:  Cascaded GRU networks and two-step training approach. (a) Forecast
of a $L-$long time signal with a $k-$long time signal using $O$(log2($L/k$))
cascaded gated recurrent units (GRUs). The output sequence of the GRU at each
stage is combined with the input sequence to serve as the input for the next-
stage GRU. In experiments, a physics-explainable GRU network is connected to a
GRU encoder in the middle stage of cascaded GRU networks that are for time
series forecast. (b) Slow and suboptimal conventional training approach
(orange path) and our fast, broadly applicable, and optimal two-step training
approach (blue path). Conventional training approach needs a large number of
training datasets. The generation of such large datasets requires long time of
numerical simulations and experiments and sometimes is not even feasible in
experimental settings. In contrast, in our multi-fidelity training approach,
the model is first trained with fast generated low-fidelity synthetic data
from physical models and then tailored to specific applications using only a
small number of high-fidelity datasets from either numerical simulations or
experiments.

Figure 2:  Dielectric metasurfaces. (a) A generic dielectric metasurface with
the unit cell consisting of four cylindrical pillars of different dimensions.
(b) The escape from local minimum in the two-step training approach. (c)
Forecast time signal using conventional training approach and two-step
approach. The clear better forecast performance is observed in the two-step
approach. (d) THz spectra obtained through the Fourier transformation of the
full time signal and the forecast time signal using short input time sequence.
Excellent agreement between forecast spectra and target spectra confirms the
superior performance of our two-step approach.

Figure 3:  Model transfer to graphene plasmonics and Laudau-polaritons. (a) –
(c) Graphene plasmonics. (a) Schematics of graphene ribbons supporting
localized plasmonic resonance. (b) Time domain and (c) corresponding frequency
domain response for predicted time signal and full target time signal,
respectively. (d) – (f) Laudau-polaritons in strongly coupled photons in an
one-dimensional photonic crystal cavity with the electron cyclotron resonance
in high-mobility two-dimensional electron gases (2DEG). (e) Time domain and
(f) corresponding frequency domain response under specific magnetic field for
predicted time signal and full target time signal, respectively.

Figure 4:  Experimental verification of models in Laudau-polaritons. Linearly
polarized THz transmittance spectra under various magnetic fields obtained
through the Fourier transformation of (a) experimentally acquired time
signals, (b) predicted time signals with a much shorter input experimental
time sequence, and (c) short input experimental time signals with zero
padding. CRA-LP and CRI branches clearly observed in experimental spectra are
reproduced in predicted spectra, where the CRA-UP branch is more noisy. The
spectra obtained from zero-padded input time signal is completely random.
Simulated transmittance colour contour map using transfer matrix method to fit
(d) experimental spectra shown in (a), (e) model-learned resonance
frequencies, and (e) predicted spectra shown in (b). Yellow stars in (d),
orange circles in (e), and gold triangles in (f) mark the peak positions
extracted from experimental spectra, learned from physics-explainable GRU
network, and extracted from predicted spectra. The extracted coupling rates
$g$ are 150.1 GHz, 145.6 GHz, and 147.9 GHz, respectively. Lorentzian fits
(red dashed line) of predicted spectra for (g) the UP peak at zero detuning
(1T) (h) the LP peak at zero detuning, and (i) the CRI mode at 3T. The
obtained full-width at half-maximum values are indicated by green arrows.

## References

  * [1] Mudelsee, M.  Trend analysis of climate time series: A review of methods.  _Earth-Sci. Rev._ 190, 310–322 (2019).
  * [2] Topol, E. J.  High-performance medicine: the convergence of human and artificial intelligence.  _Nat. Med_ 25, 44–56 (2019).
  * [3] Lim, B. & Zohren, S.  Time series forecasting with deep learning: A survey.  _arXiv preprint arXiv:2004.13408_ (2020).
  * [4] Box, G. E., Jenkins, G. M., Reinsel, G. C. & Ljung, G. M.  _Time series analysis: forecasting and control_ (John Wiley & Sons, 2015).
  * [5] Goodfellow, I., Bengio, Y., Courville, A. & Bengio, Y.  _Deep learning_ , vol. 1 (MIT press Cambridge, 2016).
  * [6] Fu, R., Zhang, Z. & Li, L.  Using lstm and gru neural network methods for traffic flow prediction.  In _2016 31st Youth Academic Annual Conference of Chinese Association of Automation (YAC)_ , 324–328 (IEEE, 2016).
  * [7] Kong, W. _et al._ Short-term residential load forecasting based on lstm recurrent neural network.  _IEEE Transactions on Smart Grid_ 10, 841–851 (2017).
  * [8] Nelson, D. M., Pereira, A. C. & de Oliveira, R. A.  Stock market’s price movement prediction with lstm neural networks.  In _2017 International joint conference on neural networks (IJCNN)_ , 1419–1426 (IEEE, 2017).
  * [9] Hyndman, R. J.  A brief history of forecasting competitions.  _Int. J. Forecast._ 36, 7–14 (2020).
  * [10] Zhou, H. _et al._ Informer: Beyond efficient transformer for long sequence time-series forecasting.  _arXiv preprint arXiv:2012.07436_ (2020).
  * [11] Jiang, J., Chen, M. & Fan, J. A.  Deep neural networks for the evaluation and design of photonic devices.  _Nat. Rev. Mater._ 1–22 (2020).
  * [12] Ma, W. _et al._ Deep learning for the design of photonic structures.  _Nat. Photonics_ 1–14 (2020).
  * [13] Rangapuram, S. S. _et al._ Deep state space models for time series forecasting.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems (NIPS 2018)_ , 7796–7805 (2018).
  * [14] Salinas, D., Flunkert, V., Gasthaus, J. & Januschowski, T.  Deepar: Probabilistic forecasting with autoregressive recurrent networks.  _Int. J. Forecast._ 36, 1181–1191 (2020).
  * [15] Willard, J., Jia, X., Xu, S., Steinbach, M. & Kumar, V.  Integrating physics-based modeling with machine learning: A survey.  _arXiv preprint arXiv:2003.04919_ (2020).
  * [16] Karniadakis, G. E. _et al._ Physics-informed machine learning.  _Nat. Rev. Phys._ 3, 422–440 (2021).
  * [17] Pellizzari, T., Gardiner, S. A., Cirac, J. I. & Zoller, P.  Decoherence, continuous observation, and quantum computing: A cavity qed model.  _Phys. Rev. Lett._ 75, 3788 (1995).
  * [18] Cirac, J. I., Zoller, P., Kimble, H. J. & Mabuchi, H.  Quantum state transfer and entanglement distribution among distant nodes in a quantum network.  _Phys. Rev. Lett._ 78, 3221 (1997).
  * [19] Duan, L.-M., Lukin, M. D., Cirac, J. I. & Zoller, P.  Long-distance quantum communication with atomic ensembles and linear optics.  _Nature_ 414, 413–418 (2001).
  * [20] Mabuchi, H. & Doherty, A.  Cavity quantum electrodynamics: coherence in context.  _Science_ 298, 1372–1377 (2002).
  * [21] Englund, D. _et al._ Controlling cavity reflectivity with a single quantum dot.  _Nature_ 450, 857–861 (2007).
  * [22] O’brien, J. L., Furusawa, A. & Vučković, J.  Photonic quantum technologies.  _Nat. Photonics_ 3, 687–695 (2009).
  * [23] Forn-Díaz, P., Lamata, L., Rico, E., Kono, J. & Solano, E.  Ultrastrong coupling regimes of light-matter interaction.  _Rev. Mod. Phys._ 91, 025005 (2019).
  * [24] Siegel, P. H.  Terahertz technology in biology and medicine.  _IEEE Trans. Microw. Theory Techn_ 52, 2438–2447 (2004).
  * [25] Tonouchi, M.  Cutting-edge terahertz technology.  _Nat. Photonics_ 1, 97–105 (2007).
  * [26] Mittleman, D. M.  Twenty years of terahertz imaging.  _Opt. Express_ 26, 9417–9431 (2018).
  * [27] Liu, H.-B., Zhong, H., Karpowicz, N., Chen, Y. & Zhang, X.-C.  Terahertz spectroscopy and imaging for defense and security applications.  _Proc. IEEE_ 95, 1514–1527 (2007).
  * [28] Lin, X. _et al._ All-optical machine learning using diffractive deep neural networks.  _Science_ 361, 1004–1008 (2018).
  * [29] Li, Y., Chen, R., Sensale-Rodriguez, B., Gao, W. & Yu, C.  Real-time multi-task diffractive deep neural networks via hardware-software co-design.  _Sci. Rep._ 11 (2021).
  * [30] Kleine-Ostmann, T. & Nagatsuma, T.  A review on terahertz communications research.  _J. Infrared Millim. Terahertz Waves_ 32, 143–171 (2011).
  * [31] Zhang, Q. _et al._ Collective non-perturbative coupling of 2D electrons with high-quality-factor terahertz cavity photons.  _Nat. Phys._ 12, 1005–1011 (2016).
  * [32] Li, X. _et al._ Vacuum bloch–siegert shift in landau polaritons with ultra-high cooperativity.  _Nat. Photonics_ 12, 324–329 (2018).
  * [33] Ulbricht, R., Hendry, E., Shan, J., Heinz, T. F. & Bonn, M.  Carrier dynamics in semiconductors studied with time-resolved terahertz spectroscopy.  _Rev. Mod. Phys._ 83, 543 (2011).
  * [34] Jepsen, P. U., Cooke, D. G. & Koch, M.  Terahertz spectroscopy and imaging–modern techniques and applications.  _Laser & Photonics Reviews_ 5, 124–166 (2011).
  * [35] Neu, J. & Schmuttenmaer, C. A.  Tutorial: An introduction to terahertz time domain spectroscopy (thz-tds).  _J. Appl. Phys._ 124, 231101 (2018).
  * [36] Nadell, C. C., Huang, B., Malof, J. M. & Padilla, W. J.  Deep learning for accelerated all-dielectric metasurface design.  _Opt. Express_ 27, 27523–27535 (2019).
  * [37] Ju, L. _et al._ Graphene plasmonics for tunable terahertz metamaterials.  _Nat. Nanotechnol._ 6, 630–634 (2011).
  * [38] Li, X. _et al._ Observation of dicke cooperativity in magnetic interactions.  _Science_ 361, 794–797 (2018).


# Design of soft magnetic materials

Ananya Renuka Balakrishna  Aerospace and Mechanical Engineering, University of
Southern California, Los Angeles, CA 90089 Corresponding author, Office: (213)
740-8762, Email<EMAIL_ADDRESS>Richard D. James  Aerospace and Engineering
Mechanics, University of Minnesota, Minneapolis, MN 55455

###### Abstract

We present a strategy for the design of ferromagnetic materials with
exceptionally low magnetic hysteresis, quantified by coercivity. In this
strategy, we use a micromagnetic algorithm that we have developed in previous
research and which has been validated by its success in solving the “Permalloy
Problem”—the well-known difficulty of predicting the composition 78.5% Ni of
lowest coercivity in the Fe-Ni system— and by the insight it provides into the
“Coercivity Paradox” of W. F. Brown. Unexpectedly, the design strategy
predicts that cubic materials with large saturation magnetization $m_{s}$ and
large magnetocrystalline anisotropy constant $\kappa_{1}$ will have low
coercivity on the order of that of Permalloy, as long as the magnetostriction
constants $\lambda_{100},\lambda_{111}$ are tuned to special values. The
explicit prediction for a cubic material with low coercivity is the
dimensionless number $(c_{11}-c_{12})\lambda_{100}^{2}/\kappa_{1}=81$ for
$\langle 100\rangle$ easy axes. The results would seem to have broad potential
application, especially to magnetic materials of interest in energy research.

## Introduction

A longstanding puzzle in materials science is understanding the origins of
magnetic hysteresis in ferromagnetic materials. Hysteresis in this domain
refers to the differing behaviors obtained when a demagnetized specimen is
subject to an increasing magnetic field to saturation vs. that obtained when
decreasing the field from saturation to zero. The effect is typically
characterized by the final (absolute) value of the magnetic field after such a
test, termed the coercivity. Informally, soft magnets have low coercivity.
This paper is concerned with the prediction of coercivity from micromagnetic
theory. An unexpected prediction of our study is that coercivity can be made
very small even in materials with large magnetocrystalline anisotropy
constant, as long as the magnetostrictive constants are tuned appropriately.

Aside from basic scientific interest on the origins of hysteresis and the
traditional application to transformers, a strategy for the discovery of new
soft magnetic materials is desirable for rapid power conversion electronics,
all-electric vehicles, and wind turbines, especially in cases where induction
motors are favored. Magnetic hysteresis has also become critical to the
adoption of proposed spintronic and storage devices, as requirements for
limiting energy consumption have moved to the forefront [1, 2]. These
requirements impact a wide range of applications, from hand held electronic
devices to storage systems and servers at data centers. While our analysis is
aimed primarily at bulk applications, the key idea that micromagnetic theory
with magnetostriction, together with a well-chosen potent defect, can be used
to predict hysteresis suggests a strategy for the lowering of coercivity also
in these small-scale applications. In fact, in certain film based devices, the
likely potent defects, such as threading dislocations in epitaxial films, are
often better characterized than in bulk material.

Currently, a widely accepted strategy to lower the hysteresis in cubic
ferromagnetic alloys is based on changing composition so as to reduce the
magnitude of the anisotropy constant $|\kappa_{1}|$. This has the effect of
flattening the graph of magnetocrystalline anisotropy energy vs. magnetization
and reducing the penalty associated with magnetization rotation. Intuitively,
this makes sense, as it apparently makes available additional low energy
pathways of an alloy in a metastable state on the shoulder of the hysteresis
loop, as an applied field is being lowered. A related idea is the known
strategy of tuning the composition so as to be precisely at the point where
two different symmetries coincide, again leading to a flattening of the of the
magnetocrystalline anisotropy energy vs. magnetization and the lowering of
hysteresis. The latter is the strategy used by Clarke and collaborators [3, 4]
that led to the particular composition of Terfenol: TbxDy1-xFe2, $x=0.3$. We
add that modern research on these RFe2 cubic Laves phase materials has focused
on the benefits of exploiting a nearby morphotropic phase boundary in these
systems to enhance magnetostrictive response under small fields [5, 6, 7].

However, these strategies cannot be the whole story behind coercivity. For
example, in the iron-nickel system a sharp drop in coercivity occurs at the
Permalloy composition, at which $\kappa_{1}=-161\mathrm{Jm^{-3}}$ [8]. Tuning
the magnetocrystalline anisotropy constant to zero in iron-nickel alloys in
fact leads to an alloy composition with noticeably higher hysteresis.
Similarly, in Sendust (Fe0.85Si0.096Al0.054 alloy) the hysteresis is minimum
when the magnetocrystalline anisotropy and magnetostriction constants are
close to zero, and not precisely at $\kappa_{1}=0$ [9]. Finally, there are
quite a few isolated examples of alloys which have very large
magnetocrystalline anisotropy constants but low coercivity: an example is the
uniaxial phase of Ni51.3Mn24.0Ga24.7 having a uniaxial magnetocrystalline
anisotropy constant of $2.45\times 10^{5}$ Jm-3 and coercivity less than $1$
kAm-1 in single crystals [10]. Another example is the Galfenol alloys
($\mathrm{Fe_{1-x}Ga_{x}}$ for $0.13\leq x\leq 0.24$) that have small magnetic
hysteresis despite their very large magnetocrystalline anisotropy constants
[11, 12]. These examples suggest that the magnetocrystalline anisotropy
constant is the not the only factor that governs hysteresis in magnetic
alloys.

The precise role of magnetostriction constants on magnetic hysteresis is not
well understood for two reasons: (a) Prior research has typically ignored the
contribution of magnetoelastic interactions on hysteresis, justified by the
small values of the magnetostriction constants, and (b) mathematical methods,
such as the linear stability analysis, overestimate the coercivity of bulk
alloys—by over three orders of magnitude in some cases—as compared to measured
experimental values [13, 14]. The latter is referred to as the “Coercivity
Paradox” [13, 15]. These factors limit our understanding of how the balance
between fundamental material constants—such as magnetocrystalline anisotropy
and magnetostriction constants—affects magnetic coercivity.

Figure 1: We model a 3D computational domain with a spike-domain
microstructure. (a) The 3D computational domain is $\Omega$ embedded inside an
oblate ellipsoid ${\cal E}$, which is several times larger in size than
$\Omega$. A spike-domain microstructure forms around a non-magnetic defect
that grows under an applied field. (b) Spike-domain microstructures are
commonly found in magnetic materials around defects such as cavities and pores
in Fe-film and Fe-Si crystal, respectively. [16, 17]. (Reprinted subfigure (b)
with permission from Ref. [16]. Copyright (2020) by the American Physical
Society)

In our recent work, we developed a computational tool based on micromagnetics,
including magnetoelastic terms, that is adapted to the prediction of
coercivity in bulk magnetic alloys [18, 19]. Micromagnetics, since its first
postulation in 1963 by W.F.Brown Jr. [13], has been applied to a wide range of
problems in ferromagnets and forms the basis to several computational
frameworks [20, 21, 22]. Our coercivity tool is based on the micromagnetics
theory [18], however differs from earlier works in the following ways: A key
feature of this tool is that it uses a large but highly localized
disturbance—in the form of a Néel-type spike domain—to predict magnetic
coercivity. Néel spike-like domains are frequently observed to form around
defects in various magnetic alloys (See [17, 23, 24, 25, 26] for examples) and
grow under an external field. In the absence of this spike-domain, the second
variation of the micromagnetic energy misses the energy barrier for
magnetization reversal [27, 13]. Instead, we model a large localized
disturbance (i.e., compute a non-linear stability analysis) to predict
coercivity on the shoulder of the hysteresis curve, see Fig. 1. Other features
of our coercivity tool are that we account for magnetoelastic interactions,
however small, in estimating coercivity of bulk magnetic alloys; and we use
the ellipsoid and reciprocal theorems to accelerate our coercivity
calculations, see [18]. Our computations are fully three-dimensional with the
magnetization evolving from (010) to (100) in all directions as one leaves the
spike domain.

To arrive at the Néel spike as a reasonable description of a potent defect, we
include in [18] a study of alternative nuclei, including rotated spikes and
multiple spikes. These spike domains serve as nuclei that grow during
magnetization reversal. In this previous study [18], we investigated the role
of defect geometry, defect orientation, and defect number on magnetic
coercivity, and we found that these features did not significantly affect
coercivity when compared to fundamental material constants. Mathematically,
they can be interpreted as localized disturbances that, under the right
combination of material constants and suitable applied field, are able to
surmount an energy barrier [28, 29, 30, 31] Small, smooth disturbances of a
homogeneous state are not able to surmount this barrier at realistic applied
fields due to the dominance of domain wall energy at small scales, and, in our
view, this is the essence of the Coercivity Paradox [18, 19]. Mathematically,
the situation is that the study of the second variation of the micromagnetic
energy misses the important energy barrier, while it is captured as a large
amplitude, but localized, disturbance.

Our viewpoint is consistent with the thesis of Pilet [28] who examines the
microscopic state on the shoulder of the hysteresis loop (far away from where
coercivity is measured) using magnetic force microscopy, and finds good
correlation with the presence of large localized disturbances. Another feature
of our tool is that we account for the magnetoelastic interactions—however
small—in estimating coercivity. Our results are in a form that is amenable to
alloy development, as in related searches for low hysteresis phase
transformations [32].

An alternative highly effective route to low hysteresis in magnetic materials
is the synthesis of nanocrystalline material [33]. For example, the powder
synthesis and rapid solidification techniques have resulted in nanocrystalline
and amorphous magnetic alloys that offer low hysteresis and enhanced
permeability [34, 2, 35]. Here we do not specifically analyze this situation,
but we think this should be possible in the nanocrystalline case, though
challenging due to the many grains that would likely have to be considered.

The aim of the present work is to identify combinations of material constants,
including saturation magnetization ($m_{s}$), magnetocrystalline anisotropy
constant ($\kappa_{1}$), elastic moduli ($c_{11},c_{12},c_{44}$), and
magnetostriction constants ($\lambda_{100},\lambda_{111}$) that give lowest
coercivity in a cubic material. We do not vary the exchange constant, which in
practice does not vary much: we also give arguments why, when combined with
the demagnetization energy, it should matter little. We find the striking
result that magnetostriction plays a critical role, and more importantly,
coercivities on the order of those found in very soft materials such as
Permalloy are predicted to be possible in materials with large
magnetocrystalline anisotropy constant, as long as the magnetostriction
constants are tuned to special values. We find that magnetocrystalline
anisotropy and magnetostriction constants play a particularly important role,
but neither has to be small.

The simulations are carried out using our newly developed coercivity tool.
Details on our micromagnetic algorithm are described in Methods. Specifically,
we apply this tool in two studies: In Study 1, we test the hypothesis that the
tuning of magnetostriction constants, in addition to the magnetocrystalline
anisotropy constant, is necessary to reduce coercivity in magnetic alloys.
Here we study how combinations of $\kappa_{1}$ and $\lambda_{100}$ affect
coercivity, and ignore the contribution from $\lambda_{111}=0$. In total, we
compute coercivity values from $N=2,163$ independent simulations. In Study 2,
we test our hypothesis that there exists a specific combination of material
constants—$\kappa_{1},\lambda_{100},\lambda_{111}$—at which magnetic
coercivity is the lowest. Here, we compute coercivity by systematically
varying the magnetocrystalline anisotropy and the magnetostriction constants,
in the ranges $0\leq\kappa_{1}\leq 2000\mathrm{Jm^{-3}}$,
$-2000\mu\epsilon\leq\lambda_{100}\leq 2000\mu\epsilon$, and
$0\leq\lambda_{111}\leq 600\mu\epsilon$ respectively ($N=605$). Overall, our
computations from over 2500 independent simulations show that the lowest
coercivity is attained when the dimensionless number
$\frac{(c_{11}-c_{12})\lambda_{100}^{2}}{2\kappa_{1}}\approx 81$. Here,
$c_{11},c_{12}$ are the elastic stiffness constants of the soft magnet
assuming a linear cubic relation. To our knowledge, this discovery of the
balance between material constants at which magnetic hysteresis is small has
not been proposed before. A theoretical analysis supporting the importance of
this dimensionless number is given below. This analysis further supports the
idea that, in other situations, a second dimensionless number
$\frac{c_{12}\lambda_{100}^{2}}{2\kappa_{1}}$ may become important, and tuning
the stiffness constants $c_{11}$ and $c_{12}$ so that both of these
dimensionless constants have particular values may be desirable [36].

## Results

### Computation of coercivity

Figure 2: The coercivity map as a function of the magnetocrystalline
anisotropy $\kappa_{1}$ and the magnetostriction constant $\lambda_{100}$. (a)
We carried out a total of $N=2163$ computations with the magnetostriction
constant $\lambda_{111}$ set to zero. We find that minimum coercivity is
achieved when $\frac{(c_{11}-c_{12})\lambda_{100}^{2}}{2\kappa_{1}}\approx
81$. For comparison, the solid red-dot corresponds to the coercivity of the
permalloy composition. Note, the parabolic relationship is a best fit
polynomial to the spread of coercivity data from our calculations. Alloys,
such as the permalloy, lie close to the vertex of this parabola and we
interpret this as the shortest distance to the parabola with a dimensionless
constant of 81. The axes (in red) indicates the crystallographic directions
along which coercivity values were measured for alloys with $\kappa_{1}>0$ and
$\kappa_{1}<0$, respectively. The 3D surface contours of the coercivity values
in (b) ‘inset A’ and (c) ‘inset B’ are shown. These plots show coercivity
wells (regions of minimum coercivity) as $\kappa_{1}\to 0$ and
$\kappa_{1}>>0$, respectively. Note, in sub-figure (b), the 3D surface contour
is discontinuous along $\kappa_{1}=0$—we attribute this discontinuity to the
transformation of easy axes from [100] to [111] crystallographic direction.

In Study 1, we test our hypothesis that the magnetostriction constants, in
addition to the magnetocrystalline anisotropy constant, is necessary to reduce
coercivity in magnetic alloys. Fig. 2(a) shows a heat map of coercivity values
as a function of the magnetocrystalline anisotropy constant $\kappa_{1}$ and
the magnetostriction constant $\lambda_{100}$. A key feature of this plot is
that the coercivity is minimum along a curve described by
$\frac{(c_{11}-c_{12})\lambda_{100}^{2}}{2\kappa_{1}}\approx 81$. As expected
the coercivity is small when $\kappa_{1}\to 0$ and the magnetostriction
constant is small $\lambda_{100}\to 0$, see Inset A. However, surprisingly, we
find that the coercivity value is also small for non-zero magnetocrystalline
anisotropy constants $\kappa_{1}>>0$ with suitable combinations of the
magnetostriction constant, see Inset B. Note the rather large
magnetocrystalline anisotropy constants being considered, well outside the
range associated to normal soft magnetism. Furthermore, the minimum coercivity
valleys are symmetric about the $\lambda_{100}=0$ axis. This symmetric
response arises from the even terms $\lambda_{100}^{2}$ in the free energy
expression. We observe a similar lowering of coercivity in magnetic alloys
with $\kappa_{1}<0$ at suitable values of $\lambda_{111}$ magnetostriction
constant (see Supplementary Figure 2).

Fig. 2(b-c) are 3D surface plots of the inset regions ‘A’ and ‘B’
respectively. The “wells” in these plots correspond to combinations of
material constants at which coercivity is minimum. Fig. 2(b) is a 3D plot of
the inset region ‘A’—here, we note a discontinuity or a jump in coercivity
values at $\kappa_{1}=0$. This discontinuity is because of the change in easy
axes for magnetic alloys with $\kappa_{1}>0$ and $\kappa_{1}<0$. For example,
we compute coercivities along $[100]$ and $[111]$ crystallographic directions
for magnetic alloys with $\kappa_{1}>0$ and $\kappa_{1}<0$, respectively.
These alloys have different values of magnetostriction constants along their
easy axes. We believe that the computed rapid change of coercivity at
$\kappa_{1}=0$ is real, and arises from the anisotropy of these
magnetostriction constants. Fig. 2(c) shows the 3D surface plot of the inset
region B. Here, coercivities comparable to that of the permalloy composition
are achieved at large magnetocrystalline anisotropy values $\kappa_{1}\approx
1700\mathrm{Jm^{-3}}$ with suitable combinations of the magnetostriction
constant $\lambda_{100}\approx 1200\mu\epsilon$.

Overall, the findings from Study 1 contradict the general understanding of
hysteresis in magnetism, i.e., the magnetocrystalline anisotropy constant
needs to be near zero for small hysteresis. Our findings show that magnetic
hysteresis (or coercivity value) is small not only when $\kappa_{1}\to 0$ but
also when $\kappa_{1}>>0$, given suitable values of the magnetostriction
constants. Fig. 2 shows that the magnetostriction constant, $\lambda_{100}$,
in addition to the magnetocrystalline anisotropy constant $\kappa_{1}$ plays
an important role in reducing coercivity in magnetic alloys.

Figure 3: The coercivity map as a function $\kappa_{1}$, $\lambda_{100}$ and
$\lambda_{111}$. (a) The minimum coercivity is achieved along a parabolic
relation, $\frac{(c_{11}-c_{12})\lambda_{100}^{2}}{2\kappa_{1}}\approx p$ .
Here, the value of coefficient $p$ depends on the value of the second
magnetostriction constant $\lambda_{111}$. For example, as $\lambda_{111}$
increases and the coercivity well loses symmetry across the $\lambda_{100}=0$
axes. (b) The 3D surface plot identifies a combination of material constants,
namely $\kappa_{1}$, $\lambda_{100}$ and $\lambda_{111}$, at which the
coercivity is small.

Thus far, we investigated coercivity values by setting one of the
magnetostriction constants $\lambda_{111}$ to be zero. In Study 2 we test the
hypothesis that there exists a specific combination of material constants,
namely $\kappa_{1},\lambda_{100},\lambda_{111}$, at which magnetic coercivity
is the lowest. We compute coercivities at every combination of
$\kappa_{1},\lambda_{100}$ and $\lambda_{111}$ in the parameter space
$0\leq\kappa_{1}\leq 2000\mathrm{Jm^{-3}}$ and
$-2000\mu\epsilon\leq\lambda_{100}\leq 2000\mu\epsilon$ and
$0\leq\lambda_{111}\leq 600\mu\epsilon$.

Fig. 3(a) shows the coercivity map as a function of $\kappa_{1}$ and
$\lambda_{100}$ for increasing values of $\lambda_{111}$. A key feature here
is that the minimum coercivity relationship
$\kappa_{1}\propto\lambda_{100}^{2}$ is unique at each value of the
$\lambda_{111}$ magnetostriction constant. For example, the minimum coercivity
valleys gradually widen and become asymmetric about $\lambda_{100}=0$, with
increasing value of the $\lambda_{111}$ magnetostriction constant. We identify
combinations of material constants at which minimum coercivity is achieved,
and then plot a 3D surface through these points in the
$\kappa_{1}-\lambda_{100}-\lambda_{111}$ parameter space, see Fig. 3(b). This
3D plot represents a surface through the material parameter space on which
coercivity is small.

### Theoretical analysis

The results of our studies above can be understood in the following way. We
begin from the free energy that we have used in the simulations of
micromagnetics in dimensional form [18],

$\displaystyle\int_{\Omega}\\{\nabla\mathbf{m}\cdot\mathrm{A\nabla\mathbf{m}+\kappa_{1}(\mathrm{m_{1}^{2}m_{2}^{2}}+\mathrm{m_{2}^{2}m_{3}^{2}}+\mathrm{m_{3}^{2}m_{1}^{2}})}$
(1) $\displaystyle+\
\frac{1}{2}[\mathbf{E}-\mathbf{E_{0}\mathrm{(\mathbf{m})}}]\cdot\mathbb{C}[\mathbf{E}-\mathbf{E_{0}\mathrm{(\mathbf{m})}}]-\mu_{0}m_{s}\mathbf{H_{\mathrm{ext}}\cdot
m}\\}\mathrm{d\mathbf{x}}+\frac{1}{2}\int_{\mathbb{R}^{3}}\mu_{0}\left|\mathbf{H_{d}}\right|^{2}\mathrm{d\mathbf{x}}.$

where ${\bf E}=\frac{1}{2}(\nabla{\bf u}+(\nabla{\bf u})^{T})$ and the energy
is to be minimized, or minimized locally, over the pair of functions ${\bf
u},{\bf m}\in H^{1}(\Omega)$. Here, ${\bf m}$ (with components on the cubic
axes $\mathrm{m_{1},m_{2},m_{3}}$) has been previously non-dimensionalized so
${\bf m}\cdot{\bf m}=1$, and the demagnetization field satisfies the
magnetostatic equations ${\bf H}_{{\bf d}}=-\nabla\zeta_{{\bf
m}},\nabla\cdot(-\nabla\zeta_{{\bf m}}+m_{s}{\bf m})=0$ on all of space, so
${\bf H}_{{\bf d}}$ has the dimensions of $m_{s}$, as does ${\bf H}_{\rm
ext}$. Please note that ${\bf m}$ is extended to $\mathbb{R}^{3}$ by making it
vanish outside $\Omega$. Typical accepted values from a large compositional
space appropriate to Fig. 3(a) including most of the Fe-Ni system are

$\displaystyle{\rm A}\sim 10^{-11}\mathrm{Jm^{-1}},\ \ \kappa_{1}\sim
0-6\times 10^{3}\mathrm{Nm^{-2}},\ \ \lambda_{100}\sim 0-2000\times 10^{-6},\
\ \lambda_{111}\sim 0-600\times 10^{-6},$ $\displaystyle m_{s}\sim
10^{6}\mathrm{Am^{-1}},\mu_{0}\sim 1.3\times 10^{-6}\mathrm{NA^{-2}},\ \
c_{11},c_{12},c_{44}\sim 10-24\times 10^{10}\mathrm{Nm^{-2}}$ (2)

A nondimensional form is obtained by dividing the micromagnetic energy by
$\mu_{0}m_{s}^{2}\sim 10^{6}\mathrm{Nm^{-2}}$, changing variables ${\bf
x}^{\prime}=\frac{1}{\ell}\,{\bf x}\in\Omega^{\prime}$, where $\ell$ is a
typical length scale and ${\bf x}^{\prime}$ is dimensionless. This gives a
typical nondimensional values of the micromagnetic coefficients

$\displaystyle\frac{A}{\ell^{2}\mu_{0}m_{s}^{2}}\sim 10^{-17},\ \
\frac{\kappa_{1}}{\mu_{0}m_{s}^{2}}\sim 10^{-3},\ \ \frac{(c_{11}-c_{12}\ {\rm
or}\ c_{44})(\lambda_{100}^{2}\ {\rm or}\
\lambda_{111}^{2})}{\mu_{0}m_{s}^{2}}\sim 10^{-5}-10^{-1},$
$\displaystyle\frac{1}{m_{s}}H_{\rm ext}\sim 10^{-3},\ \
\frac{1}{m_{s}^{2}}H_{\rm d}^{2}\sim 1$ (3)

The range of magnetostrictive coefficients is consistent with Fig. 3(a), and
for material constants with a size range, we choose a typical intermediate
value. A caveat with these numbers is the observation made by Brown [27]
(discussed also in [17]) that the formal linearization of geometrically
nonlinear micromagnetics that gives Eq. 1 implies a possible ${\bf m}$
dependence of the elastic moduli $c_{11},c_{12},c_{44}$. This dependence is
usually neglected, as is done here. Note from Eq. 3 the wide ranging values of
magnetostrictive energy and its relative importance generally.

It is seen from the nondimensionalized coefficients Eq. 3 that
magnetostriction and demagnetization energy are dominant, but it is important
to observe that each can be made negligible by suitable magnetization
distributions. The magnetostrictive energy can be made to vanish by choosing a
magnetization that satisfies curl$({\rm curl}\,{\bf E}_{0}({\bf m}))^{T}=0$,
that is, ${\bf E}_{0}({\bf m}({\bf x}))=\frac{1}{2}(\nabla{\bf u}+(\nabla{\bf
u})^{T})$ is the symmetric part of a gradient, while the demagnetization
energy vanishes on divergence-free magnetizations. Additionally, we note that
the exchange energy is one of the smallest energy contributions and its
contribution further decreases at increasing length scales. Although the
exchange constant could be affected by temperature, composition, and presence
of defects, its order of magnitude $A\approx 10^{-11}$Jm-1 does not
significantly affect large scale micromagnetic simulations. Consequently, we
do not vary the exchange constant in our calculations.

Figure 4: Growth of a spike-domain microstructure during a typical
magnetization reversal. A spike-domain forms around a non-magnetic inclusion
in the computation domain. This spike grows gradually on decreasing the
external field (a-c). At the coercivity value, the spike-domain grows abruptly
resulting in magnetization reversal (d-e). The hysteresis loop plots the
average magnetization in the computation domain $\overline{\mathrm{m}}_{1}$ as
a function of the applied field $\mathrm{H_{ext}}$.

We first explain from a theoretical viewpoint why the Néel spike attached to a
defect of the type we have chosen is a particularly potent perturbation. A key
observation comes from symmetry and holds also in the more general case of a
geometrically nonlinear magnetoelastic free energy (see Section 6 of [37]).
The observation is that domain walls involving a jump in the magnetization
${\bf m}^{+}-{\bf m}^{-}\neq 0$, where ${\bf m}^{+}$ and ${\bf m}^{-}$
minimize the magnetocrystalline anisotropy energy density and satisfy the
divergence free condition $({\bf m}^{+}-{\bf m}^{-})\cdot{\bf n}=0$ at an
interface with normal ${\bf n}$, have the property that they give rise to
strains that are perfectly mechanically compatible in the sense that ${\bf
E}_{0}({\bf m}^{+})-{\bf E}_{0}({\bf m}^{-})=\frac{1}{2}({\bf a}\otimes{\bf
n}+{\bf n}\otimes{\bf a})$ for some vector ${\bf a}$. The latter is the jump
condition implying the existence of a continuous displacement across the
interface. These conditions hold also for typical domain wall models with
remote boundary conditions given by $({\bf m}^{+},{\bf E}_{0}({\bf m}^{+}))$
and $({\bf m}^{-},{\bf E}_{0}({\bf m}^{-}))$. This argument applies not only
to materials with cubic symmetry but also to many lower symmetry systems (see
[37] for the precise statement of this result). .

To understand the relation between this symmetry argument and the Néel spike,
we substitute the form of $\mathbb{C}$ and ${\bf E}_{0}({\bf m})$ for cubic
materials into Eq. 1. The latter is,

${\bf E}_{0}({\bf m})=\frac{3}{2}\left(\begin{array}[]{ccc}\lambda_{100}({\rm
m_{1}^{2}}-\frac{1}{3})&\lambda_{111}{\rm m_{1}m_{2}}&\lambda_{111}{\rm
m_{1}m_{3}}\\\ \lambda_{111}{\rm m_{1}m_{2}}&\lambda_{100}({\rm
m_{2}^{2}}-\frac{1}{3})&\lambda_{111}{\rm m_{2}m_{3}}\\\ \lambda_{111}{\rm
m_{1}m_{3}}&\lambda_{111}{\rm m_{2}m_{3}}&\lambda_{100}({\rm
m_{3}^{2}}-\frac{1}{3})\end{array}\right)$ (4)

in the orthonormal cubic basis. We consider only the case $\kappa_{1}>0$
corresponding to $\langle 100\rangle$ easy axes, for which we have the most
data above. (The case $\kappa_{1}<0$ corresponding to $\langle 111\rangle$
easy axes is handled similarly.) Without loss of generality we also divide the
whole energy by the dimensionless number $\kappa_{1}/\mu_{0}m_{s}^{2}$. In
addition, we define a new displacement $\hat{{\bf u}}({\bf x})$ by ${{\bf
u}}({\bf x})=\lambda_{100}\hat{{\bf u}}({\bf x})$ with corresponding strain
tensor ${\bf E}=\lambda_{100}\hat{{\bf E}}$. Minimization of energy using
${\bf u}$ is equivalent to the same using $\hat{{\bf u}}$, and this
equivalence applies also to local minimization or the relative height of an
energy barrier. These changes give the explicit nondimensional form of the
free energy,

$\displaystyle\int_{\Omega^{\prime}}\mathrm{\frac{A}{\ell^{2}\kappa_{1}}\rm{m_{i,j}m_{i,j}}+(\mathit{\mathrm{{m}_{1}^{2}{m}_{2}^{2}}+\mathrm{{m}_{2}^{2}{m}_{3}^{2}}+\mathrm{{m}_{3}^{2}{m}_{1}^{2}}})}$
$\displaystyle+\frac{2c_{44}{\color[rgb]{0,0,0}\lambda_{100}^{2}}}{\kappa_{1}}\left[\left(\hat{\epsilon}_{12}-\frac{3}{2}\frac{\lambda_{111}}{\lambda_{100}}{\rm
m_{1}m_{2}}\right)^{2}+\left(\hat{\epsilon}_{13}-\frac{3}{2}\frac{\lambda_{111}}{\lambda_{100}}{\rm
m_{1}m_{3}}\right)^{2}+\left(\hat{\epsilon}_{23}-\frac{3}{2}\frac{\lambda_{111}}{\lambda_{100}}{\rm
m_{2}m_{3}}\right)^{2}\right]$
$\displaystyle+\left(\frac{(c_{11}-c_{12})\lambda_{100}^{2}}{2\kappa_{1}}\right)\left[\left(\hat{\epsilon}_{11}-\frac{3}{2}\left({\rm
m}_{1}^{2}-\frac{1}{3}\right)\right)^{2}+\left(\hat{\epsilon}_{22}-\frac{3}{2}\left({\rm
m}_{2}^{2}-\frac{1}{3}\right)\right)^{2}+\left(\hat{\epsilon}_{33}-\frac{3}{2}\left({\rm
m}_{3}^{2}-\frac{1}{3}\right)\right)^{2}\right]$
$\displaystyle+\frac{c_{12}\lambda_{100}^{2}}{2\kappa_{1}}\bigg{(}\hat{\epsilon}_{11}+\hat{\epsilon}_{22}+\hat{\epsilon}_{33}\bigg{)}^{2}$
$\displaystyle-\frac{\mu_{0}m_{s}}{\kappa_{1}}(H_{\mathrm{e}1}{\rm
m}_{1}+H_{\mathrm{e}2}{\rm m}_{2}+H_{\mathrm{e}3}{\rm m}_{3})$
$\displaystyle-\frac{\mu_{0}m_{s}}{2\kappa_{1}}(\mathbf{\mathit{H}_{\mathrm{d1}}{{\rm
m}_{\mathrm{1}}}+\mathit{H}_{\mathrm{d2}}{{\rm
m}_{\mathrm{2}}}+\mathit{H}_{\mathrm{d3}}{{\rm
m}_{\mathrm{3}}}})\,\mathrm{d{\bf x}},$ (5)

where we have used the magnetostatic equation to include the demagnetization
term in the integrand. Also,
$\hat{\epsilon}_{ij}=\frac{1}{2}(\hat{u}_{i,j}+\hat{u}_{j,i})$ and for
simplicity we have dropped the prime on ${\bf x}$. However, $\Omega^{\prime}$
remains dimensionless.

Figure 5: A 3D visualization of the needle domain. (a) This needle domain
forms around an embedded defect in our micromagnetic computations. (b) The
inset figure shows the needle domain with distinct surfaces (marked by surface
normal vectors $\mathbf{n}_{1},\mathbf{n}_{2}$). The surface $\mathbf{n}_{1}$
approximately satisfies the strain compatibility condition, however, the
surface $\mathbf{n}_{2}$ is not compatible and contributes to non-zero elastic
energy in the system. The magnetization is denoted by $\mathbf{m}$. Figure 6:
Schematic of the metastability of the Néel spikes. (a) With large applied
field the spikes are collapsed onto the nonmagnetic defect, stabilized by
(diffuse) domain wall and demagnetization energy. (b) As the applied field is
decreased the spikes grow, so as to decrease magnetostrictive and
demagnetization energies. (c) The associated energy barrier is breached at
small negative applied fields, and the spikes run to the boundary. The state
(c) is transient, and precedes the full switch to the lower energy
$\langle\bar{1}00\rangle$ magnetization; the final state is like (a) except
the spike domains are pointing NE and SW (see Fig. 4e). (d) An example
simulation of the spike domain at its transient domain state during
magnetization reversal. Labels ‘AA’, ‘BB’, ‘CC’ mark lengths across which we
examine the strain distribution in Fig. 7.

Now we make an observation about the structure of Eq. 5 relating to the
symmetry argument given above. With $\kappa_{1}>0$ as assumed, the two
magnetizations ${\bf m}^{-}=(1,0,0)$ and ${\bf m}^{+}=(0,1,0)$ satisfy $({\bf
m}^{+}-{\bf m}^{-})\cdot{\bf n}=0$, where ${\bf n}=\frac{1}{\sqrt{2}}(1,1,0)$.
Therefore, by the symmetry argument described above, ${\bf E}_{0}({\bf
m}^{+})-{\bf E}_{0}({\bf m}^{-})=\frac{1}{2}({\bf a}\otimes{\bf n}+{\bf
n}\otimes{\bf a})$ and it is indeed verified that that this holds with ${\bf
a}=\frac{3\sqrt{2}\lambda_{100}}{2}(-1,1,0)$. The values of the strains are
given in Fig. 6(c). Importantly, these choices of ${\bf m}^{+},{\bf
E}_{0}({\bf m}^{+})$ and ${\bf m}^{-},{\bf E}_{0}({\bf m}^{-})$ make the
magnetocrystalline anisotropy energy and all three magnetostrictive terms in
Eq. 5 vanish, and also give a locally divergence-free magnetization. Recalling
that the simulations are 3D (see Fig. 5), note that Fig. 6(c) can be confined
to a slab between two surfaces parallel to the plane of the page and, if ${\bf
m}({\bf x})={\bf m}^{-}$ is assigned outside the slab, then these surfaces are
also pole-free. These facts support the potency (i.e., low energy barrier)
provided by the Néel spikes.

At the transient state, the magnetic poles on the tips of the spike domains
are far apart and the geometric features of the spike-domain have evolved to
compatible interfaces. Similar symmetry arguments can be made for transient
state microstructures with multiple Néel spikes and/or other defect
geometries. In these cases the magnetic material constants dominate the energy
barrier, and the defect geometries appear to play a negligible role.

Thus we arrive at the following scenario typical of nucleation. When the
applied field is large the magnetocrystalline anisotropy energy of a large
spike is disfavored, and spike is small and collapsed near the defect, due to
the dominating influence of (diffuse) domain wall energy at small scales. As
the applied field is decreased, the magnetocrystalline anisotropy and
demagnetization energies favor the growth of the spike. A local energy maximum
is reached, beyond which an energy decreasing path is possible, leading to
complete reversal of the magnetization.

Figure 7: The strain and elastic energy distribution across the needle domain.
(a-c) The strain distribution $\hat{\epsilon}_{ij}$ across lengths ‘AA’, ‘BB’,
‘CC’ (as marked in Fig. 6(d)) in the computational domain. The strain
distribution from the simulation, $\hat{\epsilon}_{11}\approx 1$,
$\hat{\epsilon}_{22}\approx\hat{\epsilon}_{33}\approx-0.5$, is consistent with
the theoretical analysis in Fig. 6(a). The strains at the domain walls near
the spike domain vary indicating the diffuseness of the walls. The shear
strains are mostly zero throughout the computational domain. (d-f) A plot of
the three polynomials $f(\hat{\mathbf{E}},\mathbf{m})$ that accompany
coefficients $\frac{2c_{44}\lambda_{100}^{2}}{\kappa_{1}}$,
$\frac{(c_{11}-c_{12})\lambda_{100}^{2}}{2\kappa_{1}}$, and
$\frac{c_{12}\lambda_{100}^{2}}{2\kappa_{1}}$ in Eq. 5. In line with our
hypothesis the polynomials accompanying
$\frac{2c_{44}\lambda_{100}^{2}}{\kappa_{1}}$ and
$\frac{c_{12}\lambda_{100}^{2}}{2\kappa_{1}}$ make negligible energy
contributions for magnetization reversal at the transient state of the spike
domain, and the polynomial accompanying
$\frac{(c_{11}-c_{12})\lambda_{100}^{2}}{2\kappa_{1}}$ serves as an energy
barrier for magnetization reversal at the transient domain state.

The parabolic form of the locus of lowest coercivity points in Figs. 3(a) and
2(a) can now be understood heuristically. The terms involving the very small
nondimensional exchange constant and very large multiplier of the
nondimensional demagnetization energy are the only terms that involve
$\nabla{\bf m}$. These are expected to compete within the diffuse domain walls
present in these simulations, but otherwise lead only to a near divergence-
free magnetization, as suggested by the inequality of arithmetic-geometric
means in the form

$\left(\sqrt{\frac{A}{\ell^{2}\kappa_{1}}}{\rm
m_{i,j}}-\sqrt{\frac{\mu_{0}}{6\kappa_{1}}}|\nabla\zeta_{{\bf
m}}|\delta_{ij}\right)\left(\sqrt{\frac{A}{\ell^{2}\kappa_{1}}}{\rm
m_{i,j}}-\sqrt{\frac{\mu_{0}}{6\kappa_{1}}}|\nabla\zeta_{{\bf
m}}|\delta_{ij}\right)\geq 0,$ (6)

which implies, using the magnetostatic equations ${\bf H}_{{\bf
d}}=-\nabla\zeta_{{\bf m}}$ and $\nabla\cdot(-\nabla\zeta_{{\bf m}}+m_{s}{\bf
m})=0$, that

$\frac{A}{\ell^{2}\kappa_{1}}\rm{m_{i,j}m_{i,j}}-{\frac{\mu_{0}m_{s}}{2\kappa_{1}}}{\bf
H}_{{\bf d}}\cdot{\bf m}\geq
2\sqrt{\frac{A\mu_{0}}{6\ell^{2}\kappa_{1}^{2}}}({\rm div}\,{\bf
m})|\nabla\zeta_{{\bf m}}|.$ (7)

In the geometry being considered, the heuristic argument for the energy
barrier suggests the shear strains in this geometry play a minor role, and
this is supported by typical measurements of the shear strains across the
spike, Fig. 6(e), taken from the simulations. This observation, together with
the fact that the magnetization is near the easy axes $[100]$ or $[010]$ in
region of the spikes (so that $m_{i}m_{j}\approx 0,\ i\neq j$), indicates that
the constant $c_{44}{\color[rgb]{0,0,0}\lambda_{100}}^{2}/\kappa_{1}$ will
play only a minor role in determining the energy barrier.

Granted these approximations, the height of the energy barrier must then be
affected mainly by the remaining dimensionless constants

$\frac{(c_{11}-c_{12})\lambda_{100}^{2}}{2\kappa_{1}}\quad{\rm
and}\quad\frac{c_{12}\lambda_{100}^{2}}{2\kappa_{1}}.$ (8)

Again, from the simulations (Fig. 6(e)) tr$\ {\varepsilon}$ is quite small in
comparison to the multiplier of $(c_{11}-c_{12})\lambda_{100}^{2}$, and so we
expect the latter, ${(c_{11}-c_{12})\lambda_{100}^{2}}/{2\kappa_{1}}$, to be
the key dimensionless constant that governs the height of the barrier. In
simulations, the term containing this constant accounted for up to 80% of the
total magnetoelastic energy.

The value of this dimensionless constant giving the lowest energy barrier
indicated by the simulations is

$\frac{(c_{11}-c_{12})\lambda_{100}^{2}}{2\kappa_{1}}=81.0056$ (9)

As we have noticed previously [19], 78.5% Ni Permalloy does not exactly fall
on the computed parabola given by Eq. 9. This may indicate that there are
opportunities for lowering the coercivity of permalloy. It is possible that
favorable heat treatments of permalloy do just that by modifying the material
constants present in Eq. 9. An alternative possibility is that, while Eq. 9
may be the primary dimensionless constant affecting coercivity, both constants
given in Eq. 8 may play a role, indicating that the fine tuning of elastic
moduli according to their relative roles in the two constants Eq. 8 may be a
route to even lower coercivity.

## Discussion

At present, the conventional approach to develop soft magnets is to reduce the
magnetocrystalline anisotropy value to zero. Consequently the search for soft
magnets is concentrated around the $\kappa_{1}\to 0$ region. Our theory-guided
prediction suggests that in addition to the coercivity well at $\kappa_{1}\to
0$ region, there exists other regions along
$\frac{(c_{11}-c_{12})\lambda_{100}^{2}}{2\kappa_{1}}\approx 81$ at which
coercivity is small. This analytical relation between material constants gives
greater freedom for alloy development and increases the parameter space to
discover novel soft magnets. In summary, our findings show that the
magnetostriction constants, in addition to the magnetocrystalline anisotropy
constant, play an important role in governing magnetic coercivity. This was
the case in study 1 when coercivity was minimum along a parabolic relation
given by $\frac{(c_{11}-c_{12})\lambda_{100}^{2}}{2\kappa_{1}}\approx 81$ with
$\lambda_{111}=0$ in magnetic alloys. In study 2, we identified a 3D surface
topology on which the coercivity is small. Below, we discuss some limitations
of our findings and then present its potential impact to the magnetic alloy
development program.

Two features of this work limit the conclusions we can draw about the
fundamental relationship between magnetic material constants. First, we
compute coercivities assuming a simple defect structure (i.e., two spike-
domains formed around a non-magnetic inclusion) in an oblate ellipsoid. While
defect geometry has been shown to have surprisingly small effect on the
coercivity values [18], it is not known whether and how different defect
geometries, and crystallographic texture of the magnetic alloy would affect
our predictions of the parabolic relation for minimum coercivity. Second, the
presence of mechanical loads, such as residual strains or boundary loads,
could affect the delicate balance between the magnetic material constants. The
idealized stress-free conditions in our model is subject to the shortcomings
associated with the presence of microstructural inhomogeneities and surface
conditions in bulk materials. Whether accounting for these inhomogeneities
would yield comparable results in experiments is an open question. With these
reservations in mind we next discuss the impact of our findings to the alloy
development program.

A key feature of our work is that we demonstrate that the magnetostriction
constants play an important role in governing magnetic coercivity. These
findings contrast with prior research in which the magnetocrystalline
anisotropy constant has been regarded as the $\mathit{only}$ material
parameter that governs magnetic hysteresis, and the contribution from
magnetostriction constants have been largely neglected. Consequently, the
commonly accepted norm in the literature is that magnetic alloys with small
magnetocrystalline anisotropy constant have small coercivity, and magnetic
alloys with large magnetocrystalline anisotropy constant have large
coercivity. However, by accounting for both magnetocrystalline anisotropy and
magnetostriction constants, we show that magnetic alloys despite their large
$\kappa_{1}$ values have small coercivities at specific combinations of the
magnetostriction constants.

Another significant feature of our work is that we identify a relationship
$\frac{(c_{11}-c_{12})\lambda_{100}^{2}}{2\kappa_{1}}\approx 81$ between
magnetic material constants at which the coercivity is small. This generic
formula serves as a theoretical guide to the alloy development program, by
suggesting alternative combinations of material constants—beyond
$\kappa_{1}=0$—to develop soft magnets. While this relation is based on
$\langle 100\rangle$ easy axes and might differ for alloys with different easy
axes, our work is intended to guide the search for soft magnetism in ordinary
ferromagnets with large $m_{s}$. In future work, we intend to investigate
other easy axes, elastic stiffness constants, to further lower coercivity in
magnetic alloys [36]. With the recent advances in atomic-scale engineering the
compositions of magnetic alloys can be tuned atom-by-atom [38, 39]. For these
experimental approaches, our prediction of the material constant formula could
serve as a guiding principle, to engineer magnetic alloy compositions to small
hysteresis. Overall, the fundamental relationship between material constants
provides initial steps to experimentalists to discover soft magnets with high
magnetocrystalline anisotropy constants.

In conclusion, the present findings contribute to a more nuanced understanding
of how material constants, such as magnetocrystalline anisotropy and
magnetostriction constants, affect magnetic hysteresis. Specifically,
magnetoelastic interactions have been regarded to play a negligible role in
lowering magnetic coercivity. Given the current findings, we quantitatively
demonstrate that the delicate balance between magnetocrystalline anisotropy,
magnetostriction constants and the spike domain microstructure (localized
disturbance) is necessary to lower magnetic coercivity. We propose a
mathematical relationship between material constants
$\frac{(c_{11}-c_{12})\lambda_{100}^{2}}{2\kappa_{1}}\approx 81$ at which
minimum coercivity can be achieved in a material with $(100)$ easy axes. More
generally, our findings serve as a theoretical guide to discover novel
combinations of material constants that lower coercivity in magnetic alloys.

## Methods

### Micromagnetics

In our coercivity tool, we use micromagnetics theory that describes the total
free energy as a function of magnetization $\mathbf{m}$, strain $\mathbf{E}$,
and magnetostatic field $\mathbf{H_{d}}$:

$\displaystyle\mathcal{\psi}$
$\displaystyle=\int_{\Omega}\\{\nabla\mathbf{m}\cdot\mathrm{A\nabla\mathbf{m}+\kappa_{1}(\mathrm{m_{1}^{2}m_{2}^{2}}+\mathrm{m_{2}^{2}m_{3}^{2}}+\mathrm{m_{3}^{2}m_{1}^{2}})}+\frac{1}{2}[\mathbf{E}-\mathbf{E_{0}\mathrm{(\mathbf{m})}}]\cdot\mathbb{C}[\mathbf{E}-\mathbf{E_{0}\mathrm{(\mathbf{m})}}]-\mu_{0}m_{s}\mathbf{H_{\mathrm{ext}}\cdot
m}\\}\mathrm{d\mathbf{x}}$
$\displaystyle+\int_{\mathbb{R}^{3}}\frac{\mu_{0}}{2}\left|\mathbf{H_{d}}\right|^{2}\mathrm{d\mathbf{x}}.$
(10)

The form of Eq. 10 is identical to that used in our previous work, in which we
detail the meaning of the specific terms, material constants, and
normalizations [18]. For the present work, we note that the exchange energy
$\nabla\mathbf{m}\cdot\mathrm{A\nabla\mathbf{m}}$ penalizes spatial gradients
of magnetization. More generally, the anisotropy energy for cubic alloys would
have contributions from higher order energy terms (e.g.,
$\kappa_{2}(\mathrm{m_{1}^{2}m_{2}^{2}m_{3}^{2}})$), and these additional
anisotropy coefficients are likely in general to contribute to our coercivity
calculations. Although including these higher order anisotropy coefficients,
e.g. $\kappa_{2}$, could affect the easy axes of magnetization of the
material. For example at $\kappa_{1}\approx\frac{-\kappa_{2}}{2}$ other
crystallographic directions, such as $[111],[110]$, and a family of irrational
directions, would have similarly small magnetocrystalline anisotropy energy.
This warrants a systematic investigation in a future study, especially with
$\kappa_{1},\kappa_{2},\kappa_{3}$ near actual measured values. However, we do
not think that these higher order terms would significantly affect our
coercivity calculations for the following reasons: First, in our computations,
we model an oblate ellipsoid (pancake-shaped) that supports an in-plane
magnetization. This ellipsoid geometry and the defect shape penalizes out-of-
plane magnetization and thus $\mathrm{m_{3}}\approx 0$ in our calculations of
the magnetic hysteresis. Consequently, the energy contribution from the higher
order anisotropy term, $\kappa_{2}(\mathrm{m_{1}^{2}m_{2}^{2}m_{3}^{2}})$, is
negligible in our computations. Second, in the nondimensional form of the free
energy (see Eq. (5)), the energy contribution from the higher-order anisotropy
term would scale as
$\frac{\kappa_{2}}{\kappa_{1}}(\mathrm{m_{1}^{2}m_{2}^{2}m_{3}^{2}})$ with
$|\mathbf{m}|=1$. This sixth-order energy term is expected not to change the
coercivity calculations significantly. We propose to investigate the precise
role of the higher-order anisotropy terms in a future study, however, as a
first step we investigate coercivity as a function of magnetocrystalline
anisotropy $\kappa_{1}$ and magnetostriction constant $\lambda_{100}$. The
elastic energy
$\frac{1}{2}[\mathbf{E}-\mathbf{E_{0}\mathrm{(\mathbf{m})}}]\cdot\mathbb{C}[\mathbf{E}-\mathbf{E_{0}\mathrm{(\mathbf{m})}}]$
penalizes mechanical deformation away from the preferred strains, and the
external energy, $\mu_{0}\mathbf{H_{\mathrm{ext}}\cdot m}$ accounts for the
mutual interaction between magnetization moment and the applied field.
Finally, the magnetostatic energy
$\frac{\mu_{0}}{2}\left|\mathbf{H_{d}}\right|^{2}$ computed in all of space
$\mathbb{R}^{3}$ penalizes the stray fields generated by the magnetic body in
its surroundings.

We compute the evolution of the magnetization using an energy minimization
technique, the generalized Landau-Lifshitz-Ginzburg equation, see Fig. 4:

$\frac{\partial\mathbf{m}}{\partial\tau}=-\mathbf{m}\times\mathcal{H}-\alpha\mathbf{m}\times(\mathbf{m}\times\mathcal{H}).$
(11)

Here,
$\mathcal{H}=-\frac{1}{\mu_{0}m_{s}^{2}}\frac{\delta\Psi}{\delta\mathbf{m}}$
is the effective field, $\tau=\gamma m_{s}t$ is teh dimensionless time step,
$\gamma$ is the gyromagnetic ratio, and $\alpha$ is the damping constant. We
numerically solve Eq. 11 using the Gauss Siedel projection method [20], and
identify equilibrium states when the magnetization evolution converges,
$\left|\mathbf{m^{\mathrm{n+1}}}-\mathrm{\mathbf{m}^{n}}\right|<10^{-9}$. At
each iteration we compute the magnetostatic field
$\mathbf{H_{d}}=-\nabla\zeta_{\mathbf{m}}$ and the strain $\mathbf{E}$ by
solving their respective equilibrium equations:

$\displaystyle\nabla\cdot(\mathbf{\mathbf{-\nabla\zeta_{\mathbf{m}}}+\mathit{m_{s}}m})$
$\displaystyle=0\qquad\mathrm{on}\thinspace\mathbb{R}^{3}$ (12)
$\displaystyle\nabla\cdot\mathbb{C}(\mathbf{E-E}_{0})$
$\displaystyle=0\qquad\mathrm{on}\thinspace\mathcal{E}.$ (13)

The magnetostatic equilibrium condition arises from the Maxwell equations,
namely
$\nabla\times\mathbf{H_{d}}=0\to\mathbf{H_{d}}=-\nabla\zeta_{\mathbf{m}}$ and
$\nabla\cdot\mathbf{B}=\nabla\cdot(\mathbf{H_{d}+\mathit{m_{s}}m})=0$.

In our calculations, we model a finite sized domain $\Omega$ centered around a
non-magnetic defect $\Omega_{d}$, see inset images in Fig. 4. This domain is
several times smaller than the actual size of the ellipsoid $\mathcal{E},$ see
Fig. 1. We define the total demagnetization field as a sum of the local
$\mathbf{\widetilde{H}}(\mathbf{x})$ and non-local contributions
$\mathbf{\overline{H}}$. The local contribution
$\mathbf{\widetilde{H}}(\mathbf{x})$ varies spatially and accounts for the
magnetostatic fields generated from defects and other imperfections inside the
body. We calculate this local contribution by solving
$\nabla\cdot(\mathbf{\widetilde{H}+\mathit{m_{s}}\widetilde{m}})=0$ on
$\Omega$. The non-local contribution is computed as
$\mathbf{\overline{H}}=-\mathbf{N\mathit{m_{s}}\overline{m}}$. Here,
$\mathbf{N}$ is the demagnetization factor matrix that is a tabulated
geometric property of the ellipsoid. We note, from the tabulated values in
Ref. [40], the demagnetization factors for an oblate ellipsoid (pancake-
shaped) are $\mathrm{N_{11}=N_{22}=0,N_{33}=1}$. The $\mathbf{\overline{m}}$
is the constant magnetization which is defined such that
$\int_{\Omega}\mathbf{m}(\mathbf{x})\mathrm{d}\mathbf{x}=0$. This
decomposition simplifies our computational complexity, because we now model a
local domain $\Omega$ that is much smaller than modeling a domain in
$\mathbb{R}^{3}$, and yet account for the demagnetization contributions from
both the body geometry and the local defects. This decomposition is justified
in the appendix of our previous paper [18]. Both the magnetostatic and
mechanical equilibrium conditions in Eqs. 12–13 are solved in Fourier space,
see Refs. [18, 21] for further details.

### Numerical calculations

In the present work, we calibrate our micromagnetics model for the FeNi alloy
with the following material constants:
$A=10^{-11}\mathrm{Jm^{-1}},m_{s}=10^{6}\mathrm{Am^{-2}},\mu_{0}=1.3\times
10^{-6}\mathrm{NA^{-2}},c_{11}=240.8\mathrm{GPa},c_{12}=89.2\mathrm{GPa},c_{44}=75.8\mathrm{GPa}$.
The values of $\kappa_{1},\lambda_{100},\lambda_{111}$ are systematically
varied as detailed in the results section. Here, note that magnetic alloys
with positive and negative magnetocrystalline anisotropy constants, have their
easy axes along $\\{100\\}$ and $\\{111\\}$ crystallographic directions,
respectively. To accommodate this change of easy axes, we transform the energy
potential from a cubic basis, and further details of this transformation are
described in the appendix of [18].

We compute the coercivity of magnetic ellipsoids as a function of the
magnetocrystalline anisotropy and magnetostriction material constants. In each
computation, we model a domain of size $64\times 64\times 24$ containing a
defect with $16\times 16\times 6$ grid points. We choose a grid size such that
the domain walls span across 3-4 unit cells. We initialize the computational
domain with a uniform magnetization $\mathrm{\mathbf{m}=\mathrm{m_{1}}}$, and
force the magnetization inside the defect to be zero throughout the
computation, $\mathrm{\left|\mathbf{m}\right|=\mathrm{0}}$, see Fig. 4. We
apply a large external field along the easy axes,
$\mathbf{H}_{\mathrm{ext}}>>0$, and decrease it gradually in steps of
$\delta\mathbf{H}=0.25\mathrm{e_{1}}\mathrm{Oe}$. As we decrease the applied
field a spike domain forms around the defect, see Fig. 1. This spike domain
grows in size as the applied field is further lowered until a critical field
value—known as the coercivity $\mathrm{H_{\mathrm{c}}}$—at which the
magnetization vector reverses. We use this approach to predict the coercivity
of the magnetic alloys at each combination of the magnetocrystalline
anisotropy and magnetostriction constants.

Specifically, in Study 1 and Study 2 we carry out a total of, $n=2,163$ and
$n=605$, computations respectively. In these computations, we systematically
vary the material constants in the parameter space range of
$0\leq\kappa_{1}\leq 2000\mathrm{Jm^{-3}}$,
$-2000\mu\epsilon\leq\lambda_{100}\leq 2000\mu\epsilon$, and
$0\leq\lambda_{111}\leq 600\mu\epsilon$, respectively. Our investigation shows
that the minimum coercivity is attained for a parabolic relation
$\kappa_{1}\propto\lambda_{100}^{2}$, and a total of over 2,500 computations
are necessary to confirm this relationship.

Data Availability. The authors declare that the data supporting the findings
of this study are available within the paper and its supplementary information
files. Furthermore, additional data that support the findings of this study
are available from the corresponding author upon request.

Acknowledgment. The authors acknowledge the Center for Advanced Research
Computing at the University of Southern California and the Minnesota
Supercomputing Institute at the University of Minnesota for providing
resources that contributed to the research results reported within this paper.
The authors would like to thank Anjanroop Singh (University of Minnesota) for
help in checking some of the calculations. A.R.B acknowledges the support of a
Provost Assistant Professor Fellowship, Gabilan WiSE fellowship, and USC’s
start-up funds. R.D.J acknowledges the support of a Vannevar Bush Faculty
Fellowship. The authors thank NSF (DMREF-1629026) and ONR (N00014-18-1-2766)
for partial support of this work.

Author contribution. ARB and RDJ conceptualized the project, designed the
methodology, and procured funding. ARB worked on model development,
theoretical analysis, and visualization of data. RDJ worked on the theoretical
calculations. Both authors were involved with the writing of the paper.

Competing Interests. The authors declare that there are no competing
interests.

## References

  * [1] Albert Fert “Nobel Lecture: Origin, development, and future of spintronics” In _Reviews of modern physics_ 80.4 APS, 2008, pp. 1517
  * [2] Josefina M Silveyra, Enzo Ferrara, Dale L Huber and Todd C Monson “Soft magnetic materials for a sustainable and electrified world” In _Science_ 362.6413 American Association for the Advancement of Science, 2018
  * [3] JP Teter, AE Clark, M Wun-Fogle and OD McMasters “Magnetostriction and hysteresis for Mn substitutions in (TbxDy1-x)(MnyFe1-y)${}_{1}.95$” In _IEEE transactions on magnetics_ 26.5 IEEE, 1990, pp. 1748–1750
  * [4] AE Clark, JP Teter and OD McMasters “Magnetostriction “jumps” in twinned Tb0.3Dy0.7Fe1.9” In _Journal of applied physics_ 63.8 American Institute of Physics, 1988, pp. 3910–3912
  * [5] Sen Yang et al. “Large magnetostriction from morphotropic phase boundary in ferromagnets” In _Physical review letters_ 104.19 APS, 2010, pp. 197201
  * [6] Richard Bergstrom Jr et al. “Morphotropic Phase Boundaries in Ferromagnets: Tb1-xDyxFe2 Alloys” In _Physical review letters_ 111.1 APS, 2013, pp. 017203
  * [7] Cheng-Chao Hu et al. “Room-temperature ultrasensitive magnetoelastic responses near the magnetic-ordering tricritical region” In _Journal of Applied Physics_ 130.6 AIP Publishing LLC, 2021, pp. 063901
  * [8] RM Bozorth “The permalloy problem” In _Reviews of Modern Physics_ 25.1 APS, 1953, pp. 42
  * [9] M Takahashi, S Nishimaki and T Wakiyama “Magnetocrystalline anisotropy and magnetostriction of Fe-Si-Al (Sendust) single crystals” In _Journal of magnetism and magnetic materials_ 66.1 Elsevier, 1987, pp. 55–62
  * [10] R Tickle and Richard D James “Magnetic and magnetomechanical properties of Ni2MnGa” In _Journal of Magnetism and Magnetic Materials_ 195.3 Elsevier, 1999, pp. 627–638
  * [11] Jayasimha Atulasimha and Alison B Flatau “A review of magnetostrictive iron–gallium alloys” In _Smart Materials and Structures_ 20.4 IOP Publishing, 2011, pp. 043001
  * [12] Arthur E Clark, Marilyn Wun-Fogle, James B Restorff and Thomas A Lograsso “Magnetostrictive properties of Galfenol alloys under compressive stress” In _Materials transactions_ 43.5 The Japan Institute of MetalsMaterials, 2002, pp. 881–886
  * [13] William Fuller Brown “Micromagnetics” interscience publishers, 1963
  * [14] A Aharoni and S Shtrikman “Magnetization curve of the infinite cylinder” In _Physical Review_ 109.5 APS, 1958, pp. 1522
  * [15] William Fuller Brown “Magnetostatic principles in ferromagnetism” North-Holland Publishing Company, 1962
  * [16] R Schäfer and S Schinnerling “Tomography of basic magnetic domain patterns in ironlike bulk material” In _Physical Review B_ 101.21 APS, 2020, pp. 214430
  * [17] Alex Hubert and Rudolf Schäfer “Magnetic domains: the analysis of magnetic microstructures” Springer Science & Business Media, 2008
  * [18] Ananya Renuka Balakrishna and Richard D James “A tool to predict coercivity in magnetic materials” In _Acta Materialia_ 208 Elsevier, 2021, pp. 116697
  * [19] Ananya Renuka Balakrishna and Richard D James “A solution to the permalloy problem—A micromagnetic analysis with magnetostriction” In _Applied Physics Letters_ 118.21 AIP Publishing LLC, 2021, pp. 212404
  * [20] Xiao-Ping Wang, Carlos J Garcıa-Cervera and E Weinan “A Gauss–Seidel projection method for micromagnetics simulations” In _Journal of Computational Physics_ 171.1 Elsevier, 2001, pp. 357–372
  * [21] JX Zhang and LQ Chen “Phase-field microelasticity theory and micromagnetic simulations of domain structures in giant magnetostrictive materials” In _Acta Materialia_ 53.9 Elsevier, 2005, pp. 2845–2855
  * [22] Jie Wang and Jianwei Zhang “A real-space phase field model for the domain evolution of ferromagnetic materials” In _International Journal of Solids and Structures_ 50.22-23 Elsevier, 2013, pp. 3597–3609
  * [23] Felix Otto and Thomas Viehmann “Domain branching in uniaxial ferromagnets: asymptotic behavior of the energy” In _Calculus of variations and partial differential equations_ 38.1 Springer, 2010, pp. 135–181
  * [24] Felix Otto and Jutta Steiner “The concertina pattern” In _Calculus of variations and partial differential equations_ 39.1 Springer, 2010, pp. 139–181
  * [25] Lukas Döring, Radu Ignat and Felix Otto “A reduced model for domain walls in soft ferromagnetic films at the cross-over from symmetric to asymmetric wall types” In _Journal of the European Mathematical Society_ 16.7, 2014, pp. 1377–1422
  * [26] Eleonora Cinti and Felix Otto “Interpolation inequalities in pattern formation” In _Journal of Functional Analysis_ 271.11 Elsevier, 2016, pp. 3348–3392
  * [27] William Fuller Brown “Magnetoelastic interactions” Springer, 1966
  * [28] Nicolas Pilet “The relation between magnetic hysteresis and the micromagnetic state explored by quantitative magnetic force microscopy”, 2006
  * [29] Hans Knüpfer, Robert V Kohn and Felix Otto “Nucleation barriers for the cubic-to-tetragonal phase transformation” In _Communications on pure and applied mathematics_ 66.6 Wiley Online Library, 2013, pp. 867–904
  * [30] Zhiyong Zhang, Richard D James and Stefan Müller “Energy barriers and hysteresis in martensitic phase transformations” In _Acta Materialia_ 57.15 Elsevier, 2009, pp. 4332–4352
  * [31] Barbara Zwicknagl “Microstructures in low-hysteresis shape memory alloys: scaling regimes and optimal needle shapes” In _Archive for Rational Mechanics and Analysis_ 213.2 Springer, 2014, pp. 355–421
  * [32] Jun Cui et al. “Combinatorial search of thermoelastic shape-memory alloys with extremely small hysteresis width” In _Nature materials_ 5.4 Nature Publishing Group, 2006, pp. 286–290
  * [33] Som V Thomas et al. “Nanocrystallites via Direct Melt Spinning of Fe77Ni5.5Co5.5Zr7B4Cu for Enhanced Magnetic Softness” In _physica status solidi (a)_ 217.8 Wiley Online Library, 2020, pp. 1900680
  * [34] Michael E McHenry, Matthew A Willard and David E Laughlin “Amorphous and nanocrystalline materials for applications as soft magnets” In _Progress in materials Science_ 44.4 Elsevier, 1999, pp. 291–433
  * [35] Ivan Soldatov, Petru Andrei and Rudolf Schaefer “Inverted Hysteresis, Magnetic Domains, and Hysterons” In _IEEE Magnetics Letters_ 11 IEEE, 2020, pp. 1–5
  * [36] Negar Ahani, Ashwin Daware and Ananya Renuka Balakrishna “Magnetoelastic interactions reduce hysteresis in soft magnets” In _Manuscript in preparation_ , 2021
  * [37] Richard D James and Kevin F Hane “Martensitic transformations and shape-memory materials” In _Acta materialia_ 48.1 Elsevier, 2000, pp. 197–222
  * [38] S Ouazi et al. “Atomic-scale engineering of magnetic anisotropy of nanostructures through interfaces and interlines” In _Nature communications_ 3.1 Nature Publishing Group, 2012, pp. 1–9
  * [39] EN Kablov et al. “Bifurcation of magnetic anisotropy caused by small addition of Sm in (Nd1-xSmxDy)(FeCo)B magnetic alloy” In _Journal of Applied Physics_ 117.24 AIP Publishing LLC, 2015, pp. 243903
  * [40] JA Osborn “Demagnetizing factors of the general ellipsoid” In _Physical review_ 67.11-12 APS, 1945, pp. 351