[bc]Daniel Nogradi

# Mesonic decay constant and mass ratios and the conformal window

Hee Sok Chung

###### Abstract

Two particular ratios related to mesons are proposed for the study of the
conformal window in $SU(3)$ gauge theory and fundamental fermions. Lattice and
other studies indicate that the lower end, $N_{f}^{*}$, is at around 7 - 13
flavors which is a wide range without a clear consensus. Here we propose the
decay constant to mass ratios of mesons, $f_{PS,V}/m_{V}$, as a proxy since
below the conformal window lattice studies have shown that they are largely
$N_{f}$-independent while at the upper end of the conformal window they are
vanishing. The drop from the non-zero constant value to zero at $N_{f}=16.5$
might be indicative of $N_{f}^{*}$. We compute $f_{V}/m_{V}$ to N3LO and
$f_{PS}/m_{V}$ to NNLO order in (p)NRQCD. The results are unambiguously
reliable just below $N_{f}=16.5$, hence the results are expanded á la Banks-
Zaks in $\varepsilon=16.5-N_{f}$. The convergence properties of the series and
matching with the non-perturbative infinite volume, continuum and chiral
extrapolated lattice results at $N_{f}=10$ suggest that the perturbative
results might be reliable down to $N_{f}=12$. A sudden drop is observed at
$N_{f}=12$ and $N_{f}=13$ in $f_{V}/m_{V}$ and $f_{PS}/m_{V}$, respectively.

##  1 Introduction and summary

The lower end of the conformal window, usually denoted by $N_{f}^{*}$, for a
given gauge group and fermion representation, has been an elusive object of
study [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]
. Naively, one would think the lattice approach would be an ideal way to study
it because once the infinite volume, continuum and chiral limits are taken at
each flavor $N_{f}$ below the loss of asymptotic freedom, the results would be
unambiguous with a similarly unambiguous conclusion about $N_{f}^{*}$. However
it became clear that systematic effects close to the conformal window are
significantly larger than at lower flavor numbers where the models are very
similar to QCD. As a result there is not a clear consensus for $SU(3)$ and
fundamental fermions, lattice results and non ab initio methods estimate
$N_{f}^{*}$ to be somewhere in the range $7-13$ which is rather broad.

In this contribution a new approach is proposed [22]: by matching the fully
controlled lattice results obtained for low fermion numbers [23, 24, 25] and
fully controlled perturbative results obtained close to but below $N_{f}=16.5$
where asymptotic freedom is still present. The latter calculation will be the
focus of our contribution. The former are available in the range $2\leq
N_{f}\leq 10$ and the task is then to study how far down the perturbative
results can be trusted from $N_{f}=16.5$ and if they can be matched to the
last non-perturbatively obtained point at $N_{f}=10$. The particular
quantities to be investigated are dimensionless and finite ratios related to
bound states; the decay constant to mass ratios of mesons. More precisely the
ratios $f_{V}/m_{V}$ and $f_{PS}/m_{V}$ will be investigated in (p)NRQCD [26,
27, 28, 29, 30] which is the appropriate framework inside the conformal
window. These are readily available in the range $2\leq N_{f}\leq 10$ from
past lattice studies either directly, as for $f_{PS}/m_{V}$, or indirectly by
using the KSRF relations [31, 32] for $f_{V}/m_{V}$.

The setup of the perturbative calculation is as follows. We start from a CFT
close to but below $N_{f}=16.5$ which is weakly coupled as shown by Banks-Zaks
[33]. All particles are of course massless and correlation functions fall off
algebraically. A flavor singlet mass term is introduced which leads to bound
states whose masses and decay constants are proportional to $m^{\alpha}$ with
the same exponent $\alpha=1/(1+\gamma)$ related to the mass anomalous
dimension $\gamma$ [34, 35]. The constant of proportionality can be computed
perturbatively following the (p)NRQCD prescription. In the NRQCD language all
$N_{f}$ flavors are “heavy” and there are no “light” flavors and we need to
keep the purely perturbative terms only. The ratio of decay constants and
meson masses are then obtained as a series in the coupling with coefficients
depending on $N_{f}$. In the final step both the coupling and any explicit
$N_{f}$ dependence is expanded in $\varepsilon=16.5-N_{f}$ leading to constant
coefficients. The final series obtained in this way contains both powers of
$\varepsilon$ and its logarithm.

As always with a perturbative result its reliability or convergence properties
are non-trivial. For the case of $f_{V}/m_{V}$ we have N3LO results and a
comparison of the NNLO and N3LO results show that it might be reliable down to
$N_{f}=12$. We assign a theoretical error by taking the difference between the
last two available orders. The lattice result for $f_{V}$ is not available
directly, only for $f_{PS}$ but we utilise the KSRF relation to estimate
$f_{V}=\sqrt{2}f_{PS}$ on the range $2\leq N_{F}\leq 10$. Curiously, the
perturbative result at $N_{f}=12$ is compatible with the last non-perturbative
result at $N_{f}=10$ within errors. Assuming a monotonous behavior the
following picture emerges: $f_{V}/m_{V}$ is constant outside the conformal
window and drops sharply at around $N_{f}=12$ finally reaching zero at
$N_{f}=16.5$. The sudden drop might be indicative of the lower end of the
conformal window.

A similar analysis for the other ratio, $f_{PS}/m_{V}$, could not be fully
carried out because the perturbative result is only available to NNLO order
for $f_{PS}$. Nonetheless assuming the convergence properties are similar to
$f_{V}/m_{V}$ we are able to conclude similarly that the perturbative result
might be reliable down to $N_{f}=12$. A sudden drop in the ratio
$f_{PS}/m_{V}$ seems to occur at around $N_{f}=13$.

The perturbative calculation can be viewed in one of two ways. First, as
alluded to above, it may be thought of as perturbative (p)NRQCD without any
terms containing $\Lambda_{QCD}$ explicitly. Or, it is also instructive to
view it as (p)NRQED with more diagrams due to the non-abelian nature of the
interaction. If viewed this way the bound states in question are analogous to
the positronium. This view is useful because it is easy to see that all decay
constants and meson masses will be proportional to the fermion mass, just as
the decay constants and masses of positronium are proportional to the electron
mass. There is no other scale in QED than the electron mass, and there is no
other scale than the fermion mass in our calculation either since we started
from a CFT. The constant of proportionality in QED is well-known to contain
powers of $\alpha$ as well as $\log(\alpha)$ which in our case will lead to
$\varepsilon$ and $\log(\varepsilon)$.

##  2 Perturbative results

As with any perturbative calculation a running scale $\mu$ is introduced and
since we start from a CFT the natural scale is $\mu=m$, the mass of the
fermions111One could choose any scale $\mu>m$. All results will be given in
the $\overline{\rm MS\kern-0.5pt}\kern 0.5pt$ scheme and the renormalized
coupling will be denoted by $g^{2}(\mu)/(16\pi^{2})=g^{2}(m)/(16\pi^{2})=a$.
The decay constants and meson masses are expanded in $a$ leading to [36, 37,
38, 39, 40]

$\displaystyle f_{PS,V}$ $\displaystyle=$ $\displaystyle
b_{0}\,m\,a^{3/2}\left(1+b_{10}a+b_{11}a\log a+b_{20}a^{2}+b_{21}a^{2}\log
a+b_{22}a^{2}\log^{2}a+O(a^{3})\right)$ $\displaystyle m_{V}$ $\displaystyle=$
$\displaystyle c_{0}\,m\left(1+c_{20}a^{2}+c_{30}a^{3}+c_{31}a^{3}\log
a+O(a^{4})\right)\;.$ (1)

The leading term in the mass, $c_{0}=2$, just follows from having a free
fermion and anti-fermion pair whereas the first correction $c_{20}<0$ is
familiar from the quantum mechanical binding energy in a Coulomb potential.
Further radiative corrections are systematically obtained using (p)NRQCD. The
leading expression, $b_{0}$, for the decay constants is proportional to the
ground state wave function at the origin. The explicit form of the
corrections, NNLO for $f_{PS}$ and N3LO for $f_{V},m_{V}$ can be found in
[22].

In the ratio $m$ drops out and the massless limit takes the running coupling
to the fixed point $a(m)\to a_{*}$ which can be expanded in
$\varepsilon=16.5-N_{f}$ and is known to 5-loops [41, 42, 43, 44, 45, 46],

$\displaystyle
a_{*}=\varepsilon\,(e_{0}+e_{1}\,\varepsilon+e_{2}\,\varepsilon^{2}+e_{3}\,\varepsilon^{3}+\ldots)\,,$
(2)

with some coefficients $e_{i}$. Combining this expansion with (2) leads to the
final results,

$\displaystyle\frac{f_{V}}{m_{V}}$ $\displaystyle=$
$\displaystyle\varepsilon^{3/2}C_{0}\left(1+\sum_{n=1}^{3}\sum_{k=0}^{n}C_{nk}\,\varepsilon^{n}\log^{k}\varepsilon+O(\varepsilon^{4})\right)$
(3) $\displaystyle C_{0}$ $\displaystyle=$ $\displaystyle 0.005826678$
$\displaystyle C_{10}$ $\displaystyle=$ $\displaystyle 0.4487893\qquad
C_{11}=-0.2056075$ $\displaystyle C_{20}$ $\displaystyle=$ $\displaystyle
0.2444502\qquad C_{21}=-0.1624891\quad C_{22}=0.03522870$ $\displaystyle
C_{30}$ $\displaystyle=$ $\displaystyle 0.10604(3)\;\;\,\quad
C_{31}=-0.1128420\quad C_{32}=0.03695458\quad C_{33}=-0.005633665$

for the vector case and

$\displaystyle\frac{f_{PS}}{m_{V}}$ $\displaystyle=$
$\displaystyle\varepsilon^{3/2}C_{0}\left(1+\sum_{n=1}^{2}\sum_{k=0}^{n}D_{nk}\,\varepsilon^{n}\log^{k}\varepsilon+O(\varepsilon^{3})\right)$
(4) $\displaystyle D_{10}$ $\displaystyle=$ $\displaystyle 0.4654041\quad
D_{11}=-0.2056075$ $\displaystyle D_{20}$ $\displaystyle=$ $\displaystyle
0.2845697\quad D_{21}=-0.1737620\quad D_{22}=0.03528692$

for the pseudo-scalar case.

Note that the coefficients in (3) and (4) are well-behaved as the order grows,
in contrast to the generally factorially growing perturbative coefficients
such as (2) and (2). Furthermore, since $f_{PS,V}/m_{V}$ are finite, RG-
invariant, scheme independent physical quantities, all coefficients in (3) and
(4) are scheme independent as well222The coefficient $C_{30}$ is only
available in numerical form with some uncertainty at the moment.. These
coefficients are the main results of our work.

Figure 1: Left: the $f_{V}/m_{V}$ ratio in increasing perturbative order. The
non-perturbative result from combined lattice calculations [23, 24, 25] and
the KSRF-relation is also shown. The smaller error band corresponds to the
uncertainty of the lattice calculation, the wider one combines this with a
conservative estimate of the uncertainty of the KSRF-relation itself. Right:
The corresponding results for $f_{PS}/m_{V}$.

##  3 Matching low $N_{f}$ and high $N_{f}$

The increasing orders for the two ratios are shown in figure 1. Clearly, the
deviation between the NNLO and N3LO results of $f_{V}/m_{V}$ for $N_{f}\geq
12$ is not substantial. Quantitatively, in the range $11.9\leq N_{f}\leq
12.1$, the deviation between the NNLO and N3LO results is at most 4%, or in
the range $11.5\leq N_{f}\leq 12.5$ at most 13%. We thus conclude that in the
region of interest, $N_{f}\sim 12$, the N3LO result is robust and reliable. We
take as an estimate of the neglected higher orders the difference between the
last two available orders.

In order to compare with the lattice results we would need $f_{V}$, however
only $f_{PS}$ was measured directly [23, 24, 25]. Here we use the KSRF
relation $f_{V}=\sqrt{2}f_{PS}$ originating in vector meson universality to
estimate $f_{V}$ and assign a $12\%$ uncertainty which holds in QCD. Somewhat
unexpectedly the perturbative result at $N_{f}=12$ matches the last non-
perturbative lattice result almost exactly; see left panel of figure 2.

A similar analysis unfortunately cannot be completed for $f_{PS}/m_{V}$
because $f_{PS}$ is only available to NNLO order. One may nonetheless assume
that the theoretical uncertainty is similar to that of $f_{V}$ making a
comparison with the lattice results at low $N_{f}$ feasible; see right panel
of 2.

In both cases a match between the low $N_{f}$ and high $N_{f}$ regions seems
plausible. Assuming a monotonous behavior and trusting the lattice results
below $N_{f}=10$ and the perturbative ones above $N_{f}=12$ or $13$ leaves
only a narrow range to be interpolated or more rigorously calculated in future
lattice work. The following picture seems to emerge: a mostly
$N_{f}$-independent flat curve drops sharply at around $N_{f}=12$ and
$N_{f}=13$ for the two ratios, respectively. This sudden change in behavior
might be indicative of $N_{f}^{*}$, the lower end of the conformal window.

Figure 2:  Non-perturbative lattice results in the range $2\leq N_{f}\leq 10$
and the perturbative ones where they seem reliable. Assuming a monotonous
behavior only a small range needs to be interpolated. Left: $f_{V}/m_{V}$, the
wider error bands for $2\leq N_{f}\leq 10$ includes the error from the usage
of the KSRF relation. The error of the perturbative curve is estimated from
the difference of the last two available orders. Right: $f_{PS}/m_{V}$, the
error of the perturbative result is estimated from that of $f_{V}/m_{V}$.

##  4 Conclusion and outlook

In this contribution we presented a new approach to shed light on the
emergence of conformal behavior from chirally broken dynamics as the flavor
number increases which combines both perturbative and non-perturbative
information. Well chosen dimensionless quantities were presented which can be
easily measured in lattice calculations in the chiral limit and which can also
be computed in perturbation theory, again in the massless limit. Current
lattice calculations are able to provide unambiguous results for low $N_{f}$
and the perturbative results are reliable at high $N_{f}\leq 16.5$. Curiously,
the two approaches seem to match at around $N_{f}=12,13$ where a sudden change
in behavior as a function of $N_{f}$ is observed.

The results can be improved in a number of ways. First, direct lattice results
at $N_{f}=12$ would be very useful. The difficulty is controlling all 3
sources of systematic effects, finite volume, finite lattice spacing and
finite mass which certainly would lead to very costly calculations. Second,
the dominant source of uncertainty of $f_{V}/m_{V}$ was the use of the KSRF
relation, which in principle could be eliminated once $f_{V}$ is measured
directly on the lattice. Third, currently the highest (p)NRQCD order for
$f_{PS}$ is NNLO which in principle could be extended to N3LO, similarly to
$f_{V}$. However going beyond N3LO order for any quantity does not seem
feasible in the near future since the 6-loop $\beta$-function would be needed
for that.

## References

  * [1]
  * [2] T. Appelquist, K. D. Lane and U. Mahanta, Phys. Rev. Lett. 61, 1553 (1988)
  * [3] A. G. Cohen and H. Georgi, Nucl. Phys. B 314, 7-24 (1989)
  * [4] F. Sannino and K. Tuominen, Phys. Rev. D 71, 051901 (2005) [arXiv:hep-ph/0405209 [hep-ph]].
  * [5] D. D. Dietrich and F. Sannino, Phys. Rev. D 75, 085018 (2007) [arXiv:hep-ph/0611341 [hep-ph]].
  * [6] T. Appelquist et al., Phys. Rev. Lett. 100, 171607 (2008), 102, 149902 (2009) [arXiv:0712.0609 [hep-ph]].
  * [7] A. Armoni, Nucl. Phys. B 826, 328-336 (2010) [arXiv:0907.4091 [hep-ph]].
  * [8] T. Appelquist, G. T. Fleming and E. T. Neil, Phys. Rev. D 79, 076010 (2009) [arXiv:0901.3766 [hep-ph]].
  * [9] Z. Fodor et al., Phys. Lett. B 681, 353-361 (2009) [arXiv:0907.4562 [hep-lat]].
  * [10] M. T. Frandsen, T. Pickup and M. Teper, Phys. Lett. B 695, 231-237 (2011) [arXiv:1007.1614 [hep-ph]].
  * [11] Z. Fodor et al., Phys. Lett. B 703, 348-358 (2011) [arXiv:1104.3124 [hep-lat]].
  * [12] A. Hasenfratz, Phys. Rev. Lett. 108, 061601 (2012) [arXiv:1106.5293 [hep-lat]].
  * [13] Y. Aoki et al., Phys. Rev. D 86, 054506 (2012) [arXiv:1207.3060 [hep-lat]].
  * [14] Z. Fodor et al., Phys. Rev. D 94, no.9, 091501 (2016) [arXiv:1607.06121 [hep-lat]].
  * [15] A. Hasenfratz and D. Schaich, JHEP 02, 132 (2018) [arXiv:1610.10004 [hep-lat]].
  * [16] D. Nogradi and A. Patella, Int. J. Mod. Phys. A 31, no.22, 1643003 (2016) [arXiv:1607.07638 [hep-lat]].
  * [17] Z. Fodor et al., Phys. Lett. B 779, 230-236 (2018) [arXiv:1710.09262 [hep-lat]].
  * [18] A. Hasenfratz, C. Rebbi and O. Witzel, Phys. Lett. B 798, 134937 (2019) [arXiv:1710.11578 [hep-lat]].
  * [19] B. S. Kim, D. K. Hong and J. W. Lee, Phys. Rev. D 101, no.5, 056008 (2020) [arXiv:2001.02690 [hep-ph]].
  * [20] J. W. Lee, Phys. Rev. D 103, no.7, 076006 (2021) [arXiv:2008.12223 [hep-ph]].
  * [21] K. Rummukainen and K. Tuominen, Universe 8, no.3, 188 (2022)
  * [22] H. S. Chung and D. Nogradi, Phys. Rev. D 107, no.7, 074039 (2023) [arXiv:2302.06411 [hep-ph]].
  * [23] D. Nogradi and L. Szikszai, JHEP 05, 197 (2019) [erratum: JHEP 06, 031 (2022)] [arXiv:1905.01909 [hep-lat]].
  * [24] D. Nogradi and L. Szikszai, PoS LATTICE2019, 237 (2019) [arXiv:1912.04114 [hep-lat]].
  * [25] A. Y. Kotov et al., JHEP 07, 202 (2021) [erratum: JHEP 06, 032 (2022)] [arXiv:2107.05996 [hep-lat]].
  * [26] W. E. Caswell and G. P. Lepage, Phys. Lett. B 167, 437-442 (1986)
  * [27] G. T. Bodwin et al., Phys. Rev. D 51, 1125-1171 (1995), 55, 5853 (1997) [arXiv:hep-ph/9407339 [hep-ph]].
  * [28] A. Pineda and J. Soto, Nucl. Phys. B Proc. Suppl. 64, 428-432 (1998) [arXiv:hep-ph/9707481 [hep-ph]].
  * [29] N. Brambilla, A. Pineda, J. Soto and A. Vairo, Nucl. Phys. B 566, 275 (2000) [arXiv:hep-ph/9907240 [hep-ph]].
  * [30] N. Brambilla, A. Pineda, J. Soto and A. Vairo, Rev. Mod. Phys. 77, 1423 (2005) [arXiv:hep-ph/0410047 [hep-ph]].
  * [31] K. Kawarabayashi and M. Suzuki, Phys. Rev. Lett. 16, 255 (1966)
  * [32] Riazuddin and Fayyazuddin, Phys. Rev. 147, 1071-1073 (1966)
  * [33] T. Banks and A. Zaks, Nucl. Phys. B 196, 189-204 (1982)
  * [34] L. Del Debbio and R. Zwicky, Phys. Rev. D 82, 014502 (2010) [arXiv:1005.2371 [hep-ph]].
  * [35] L. Del Debbio and R. Zwicky, Phys. Lett. B 700, 217-220 (2011) [arXiv:1009.2894 [hep-ph]].
  * [36] A. A. Penin et al., Nucl. Phys. B 699, 183-206 (2004), 829, 398-399 (2010) [arXiv:hep-ph/0406175 [hep-ph]].
  * [37] B. A. Kniehl et al., Phys. Lett. B 638, 209-213 (2006) [arXiv:hep-ph/0604072 [hep-ph]].
  * [38] G. T. Bodwin, H. S. Chung, J. Lee and C. Yu, Phys. Rev. D 79, 014007 (2009) [arXiv:0807.2634 [hep-ph]].
  * [39] A. A. Penin and N. Zerf, JHEP 04, 120 (2014) [arXiv:1401.7035 [hep-ph]].
  * [40] H. S. Chung, JHEP 12, 065 (2020) [arXiv:2007.01737 [hep-ph]].
  * [41] O. V. Tarasov, A. A. Vladimirov and A. Y. Zharkov, Phys. Lett. B 93, 429-432 (1980)
  * [42] S. A. Larin and J. A. M. Vermaseren, Phys. Lett. B 303, 334-336 (1993) [arXiv:hep-ph/9302208 [hep-ph]].
  * [43] T. van Ritbergen et al., Phys. Lett. B 400, 379-384 (1997) [arXiv:hep-ph/9701390 [hep-ph]].
  * [44] M. Czakon, Nucl. Phys. B 710, 485-498 (2005) [arXiv:hep-ph/0411261 [hep-ph]].
  * [45] P. A. Baikov et al., Phys. Rev. Lett. 118, no.8, 082002 (2017) [arXiv:1606.08659 [hep-ph]].
  * [46] F. Herzog, B. Ruijl, T. Ueda, J. A. M. Vermaseren and A. Vogt, JHEP 02, 090 (2017) [arXiv:1701.01404 [hep-ph]].


# Learning dislocation dynamics mobility laws from large-scale MD simulations

Nicolas Bertin<EMAIL_ADDRESS>Vasily V. Bulatov  Fei Zhou  Lawrence
Livermore National Laboratory, Livermore, CA, USA

###### Abstract

The computational method of discrete dislocation dynamics (DDD), used as a
coarse-grained model of true atomistic dynamics of lattice dislocations, has
become of powerful tool to study metal plasticity arising from the collective
behavior of dislocations. As a mesoscale approach, motion of dislocations in
the DDD model is prescribed via the mobility law; a function which specifies
how dislocation lines should respond to the driving force. However, the
development of traditional “hand-crafted” mobility laws can be a cumbersome
task and may involve detrimental simplifications. Here we introduce a machine-
learning (ML) framework to streamline the development of data-driven mobility
laws which are modeled as graph neural networks (GNN) trained on large-scale
Molecular Dynamics (MD) simulations of crystal plasticity. We illustrate our
approach on BCC tungsten and demonstrate that our GNN mobility implemented in
large-scale DDD simulations accurately reproduces the challenging
tension/compression asymmetry observed in ground-truth MD simulations while
correctly predicting the flow stress at lower straining rate conditions unseen
during training, thereby demonstrating the ability of our method to learn
relevant dislocation physics. Our DDD+ML approach opens new promising avenues
to improve fidelity of the DDD model and to incorporate more complex
dislocation motion behaviors in an automated way, providing a faithful proxy
for dislocation dynamics several orders of magnitude faster than ground-truth
MD simulations.

###### keywords:

Dislocation mobility , Dislocation dynamics , Graph neural networks , Machine
learning

††journal: arXiv

##  1 Introduction

In metals, plastic deformation is ordinarily defined by the motion and
interaction of dislocation lines through the lattice. While direct Molecular
Dynamics (MD) simulations of crystal plasticity are now within reach to
provide insights on the collective behavior of dislocations and metal strength
[1, 2, 3, 4], their scale still commands considerable computing resources.
Alternatively, the method of Discrete Dislocation Dynamics (DDD), used as an
expedient mesoscale proxy, has been widely regarded as a powerful tool to
connect the evolution of the dislocation microstructure to the macroscopic
response of crystals, thereby bridging the gap between mesoscale and continuum
approaches in metal plasticity [5, 6, 7, 8, 9, 10, 11]. In contrast to MD
which simulates “all the atoms”, only the dislocation microstructure is
considered in DDD, which is represented as a set of dislocation segments
inter-connected through nodes and generally evolved with an over-damped
equation of motion [12]

$\boldsymbol{V}_{i}=\frac{d\boldsymbol{r}_{i}}{dt}=\mathcal{M}\left[\boldsymbol{F}_{i}\left(\\{\boldsymbol{r}_{j}\\},\\{\boldsymbol{b}_{jk}\\}\right)\right]$
(1)

where $\boldsymbol{r}_{j}\in\mathcal{V}$ and
$\boldsymbol{b}_{jk}\in\mathcal{E}$ are the dislocation node positions and
Burgers vector connectivity, respectively defining the set of vertices and
edges of the dislocation graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$. In Eq.
(1), $\mathcal{M}$ is the mobility function that relates the nodal force
$\boldsymbol{F}_{i}$ exerted on the dislocation line at node $i$ to its
resulting nodal velocity $\boldsymbol{V}_{i}$. The mobility law is the core
ingredient of the DDD model which encapsulates the physics of the simulation
(analogous to the interatomic potential in MD), i.e. by which DDD
practitioners specify how dislocations should move in response to the driving
force.

In principle, the DDD mobility needs to be calibrated to a reference behavior
in order for simulations to produce meaningful predictions. Given the relative
difficulty to extract force-velocity relations from experiments, a widely used
strategy has been to calibrate mobility laws using MD simulations of
individual dislocations. Within this approach, a mobility law is typically
constructed by postulating a “hand-crafted” functional form based on physical
knowledge and intuition, and quantifying its parameters by running MD
calculations of isolated dislocations [13, 14, 15], e.g. by computing the
velocity of dislocations of distinct types (for instance, edge and screw)
under varying levels of stresses and temperatures [16, 17, 18, 19, 20],
extracting nucleation energy barriers and critical stresses for motion
controlled by the kink-pair mechanism [21, 22, 15], or extracting dislocation
drag by analysis of the thermal fluctuations spectrum [23].

However, while this conventional workflow has been widely adopted, it presents
several important drawbacks. First, constructing mobility laws can be a time-
consuming task, and differences in metals crystallography and dislocation
behaviors make it difficult to fully automate the process. Second, this
approach necessarily involves simplifications. In particular, it remains
unclear if and how behaviors extracted from straight, isolated dislocations in
MD can be transferred [24] to large-scale DDD simulations in which
dislocations form dense and intricate networks of curved lines [25]. Third, to
the authors’ knowledge, no work has yet verified how DDD predictions using
mobility laws constructed in the way described above would compare to direct,
large-scale MD simulations performed using the same interatomic potential used
in the calibration.

In this work, we introduce the DDD+ML framework to streamline the construction
of DDD mobility laws using a machine learning (ML) approach akin to the task
of developing interatomic potentials for MD simulations. For this, we propose
to define the mobility law in Eq. (1) as a graph neural network (GNN) operator
trained to learn the evolution of dislocation network trajectories produced by
large-scale MD simulations. Here, we therefore view DDD as a coarse-grained
model of direct MD simulations of crystal plasticity, which we regard as the
ground-truth that we wish to emulate. Thus, a faithful mobility law should
result in DDD+ML predictions that closely match those of direct MD simulations
performed under identical conditions. Using body-centered cubic (BCC) W as a
testbed material, we show that our trained GNN-based mobility implemented in
our DDD model can learn relevant dislocation physics and is able to reproduce
well the strength behavior and tension/compression asymmetry predicted in
large-scale MD simulations. We believe that our new GNN-based mobility
framework provides a novel, unbiased approach to accurately incorporate
dislocation physics in the DDD model while paving the way to automated
workflows for constructing mobility laws of arbitrary complexity.

##  2 Results

###  2.1 DDD+ML model

Figure 1: DDD+ML framework. (a) Schematic of two consecutive dislocation
networks $\mathcal{G}^{t}$ and $\mathcal{G}^{t+\Delta t}$ extracted from a MD
trajectory using DXA at times $t$ (solid lines) and $t+\Delta t$ (dashed
lines). As extracted, the configurations contain nodal positions
$\boldsymbol{r}=\\{\boldsymbol{r}_{i}\\}$ and Burgers vector connectivity
$\boldsymbol{b}=\\{\boldsymbol{b}_{ij}\\}$ but (i) contain no information
about the line forces driving their evolution, and (ii) nodal velocity vectors
generally cannot be unambiguously defined due to the nonisomorphism between
both networks. For instance, colored arrows show three possible ways to define
the ground-truth velocity $\boldsymbol{V}_{2}$ at node $2$. (b) Schematic of
the inference and training loop for the GNN-based mobility law. (c)
Illustration of the Nye’s tensor field-matching approach to circumvent the
ill-defined matching of velocities. As exemplified here, the procedure is
agnostic to line discretization and network topology.

Taking advantage of the graph nature of dislocation networks, we view the
construction of the DDD mobility law in Eq. (1) as a ML task defined on a
graph network, Fig. 1. Very recently, we successfully developed a machine-
learned mobility from a simplified dislocation-obstacle DDD system [26]. Here,
we propose to model the mobility function with a GNN trained to predict the
evolution of subsequent dislocation network configurations, e.g. as extracted
from large-scale (ground-truth) MD trajectories using the dislocation
extraction algorithm (DXA) [27, 28].

In practice however, the task is complex and presents several challenges that
need to be addressed. To illustrate the matter, consider two consecutive
dislocation graphs $\mathcal{G}^{t}$ and $\mathcal{G}^{t+\Delta t}$ at time
$t$ and $t+\Delta t$, respectively, as depicted in Fig. 1a. As implied by Eq.
(1), training a ML mobility law entails learning the relation between the
(input) nodal force and (output) nodal velocity vectors given the geometry and
topology of the dislocation network. Yet, by themselves dislocation networks
extracted with DXA (i) contain no information about the driving force on the
dislocation lines, and (ii) the dislocation velocities are only implicitly
defined by comparing dislocation line configurations attained at consecutive
times, and thus not explicitly contained in the data. This is because the
topology of dislocation networks constantly evolves (e.g. as a result of
dislocation intersections and core reactions) while the discretization of
dislocation links is essentially arbitrary, i.e. subsequent dislocation graphs
are generally nonisomorphic. In other words, there does not exist a one-to-one
correspondence between dislocation nodes in consecutive network configurations
and velocity vectors cannot be unambiguously defined. As such, learning on
nodal velocity vectors is not a well-defined task for ML.

We propose to address both these challenges of unknown ground-truth force and
velocity vectors in the following way. First, we estimate the nodal forces
$\boldsymbol{F}_{i}$ on nodes $i$ of the MD dislocation networks using the
force calculation procedure of the DDD model

$\boldsymbol{F}_{i}\sim\boldsymbol{F}_{i}^{\rm DDD}=\boldsymbol{F}_{i}^{\rm
app}+\boldsymbol{F}_{i}^{\rm lr}+\boldsymbol{F}_{i}^{\rm
sr}+\boldsymbol{F}_{i}^{\rm core}$ (2)

where the force $\boldsymbol{F}_{i}^{\rm DDD}$ calculated with DDD is the sum
of several contributions, namely the force due to the applied loading and
boundary conditions $\boldsymbol{F}_{i}^{\rm app}$, the long-range
$\boldsymbol{F}_{i}^{\rm lr}$ and short-range $\boldsymbol{F}_{i}^{\rm sr}$
elastic interactions, and the force associated with the dislocation core
$\boldsymbol{F}_{i}^{\rm core}$. Given the various approximations involved in
force calculations in the DDD model, this force is only an _estimate_ of the
true (ground-truth) but unknown force $\boldsymbol{F}_{i}^{\rm MD}$ driving
the dislocations in the MD simulations. Mathematically, we can write this
approximation exactly as

$\boldsymbol{F}_{i}^{\rm MD}=\boldsymbol{F}_{i}^{\rm
DDD}+\boldsymbol{F}_{i}^{\rm corr}$ (3)

where $\boldsymbol{F}_{i}^{\rm corr}$ is a correction term accounting for the
error made in approximating the true force with the DDD model. Assuming that
the applied and long-range components of the force are smoothly varying
functions that can be accurately computed with DDD, it follows that the
correction term in Eq. (3) is a rather local function of dislocation network,
which can itself be easily learned by a GNN model. Thus, we propose to define
a GNN-based mobility law $\mathcal{M}^{\rm GNN}_{\boldsymbol{\theta}}$ with
learnable parameters $\boldsymbol{\theta}$,

$\boldsymbol{V}_{i}=\mathcal{M}^{\rm
GNN}_{\boldsymbol{\theta}}\left(\boldsymbol{F}_{i}^{\rm
DDD},\\{N_{j}\\},\\{E_{jk}\\}\right),$ (4)

and let the model learn the resulting nodal velocity vector
$\boldsymbol{V}_{i}$ as a function of the estimated force
$\boldsymbol{F}_{i}^{\rm DDD}$ and local graph neighborhood of nodes $j$ and
segments $jk$ with attributes $N_{j}$ and $E_{jk}$, respectively (see
Methods). As defined, the mobility function in Eq. (4) thus fully subsumes the
calculation of the as-yet unspecified local force correction. Since in DDD
simulations forces are always computed using Eq. (2), this formulation further
ensures that the approach is self-consistent when employing the trained
mobility in rollout simulations.

We now turn our attention to the second issue, namely the absence of well-
defined nodal velocity vectors in the ground-truth MD data. Here we note that
this issue is essentially identical to the task of tracing dislocation motion
between consecutive MD dislocation snapshots, for which we have recently
introduced the sweep-tracing algorithm (STA) [3]. By using a dual line/field
representation of the dislocation networks, STA can “reconnect” successive
networks in a way that is independent of line discretization and agnostic to
network topology. Conceptually, the task is defined as an optimization problem
that seeks to minimize the distance between Nye’s tensor [29] field
representations of consecutive graphs [3]. Following the STA approach, the
training of the GNN-based mobility operator in Eq. (4) is achieved by
minimizing the following loss function

$\mathcal{L}(\boldsymbol{\theta})=\sum_{s}\sum_{\boldsymbol{g}}\sum_{kl}\Bigl{[}\alpha^{\boldsymbol{g}}_{kl}\left(\boldsymbol{r}^{t_{s}}+\boldsymbol{\mathcal{M}}_{\boldsymbol{\theta}}^{\rm
GNN}\Delta t_{s},\boldsymbol{b}^{t_{s}}\right)\\\
-\alpha^{\boldsymbol{g}}_{kl}\left(\boldsymbol{r}^{t_{s}+\Delta
t_{s}},\boldsymbol{b}^{t_{s}+\Delta t_{s}}\right)\Bigr{]}^{2}$ (5)

over the set of training examples $s$, where $\alpha_{kl}^{\boldsymbol{g}}$
are the components of the Nye’s tensor field numerically computed at grid
point $\boldsymbol{g}$ using the method introduced in [30]. During training,
the mobility law thus learns to predict nodal velocity vectors
$\boldsymbol{V}_{i}$ which, when applied to configurations at time $t_{s}$,
best match the tensor field representation of the ground-truth evolution of
the networks at time $t_{s}+\Delta t_{s}$, Fig. 1c. In other words, the field-
matching loss of Eq. (5) circumvents of the seemingly straightforward but ill-
defined matching of velocities.

We point out that the calculation of the Nye’s tensor field over networks is
only required to compute and back-propagate mobility parameters
$\boldsymbol{\theta}$ information during training. Once trained, only the GNN-
based mobility function in Eq. (4) needs to be evaluated at inference time,
Fig. 1b. Our GNN model and training protocol are implemented using PyTorch
[31].

###  2.2 Validation on DDD trajectories

To validate our approach, we first apply our method to DDD trajectories to
evaluate the ability of our framework to learn a simple “hand-crafted” DDD
mobility law. Although the approach was described for DXA configurations
extracted from MD in §2.1, it can be applied equally well to DDD
configurations. Doing so can be insightful since, in contrast to MD
simulations, we have complete control over traditional mobility laws in DDD
and can prescribe any desired ground-truth force-velocity relation, and then
examine how well our approach can recover it. Here we choose a generic linear
mobility function developed for BCC metals [32] and implemented in our DDD
code ParaDiS [11], which can be expressed as

$\boldsymbol{V}_{i}=\boldsymbol{\mathcal{M}}^{\rm
DDD}(M_{s},M_{e})\cdot\boldsymbol{F}_{i}^{\rm DDD}$ (6)

where $\boldsymbol{\mathcal{M}}^{\rm DDD}$ is a matrix parameterized by scalar
mobility coefficients $M_{s}=20~{}({\rm Pa}\cdot{\rm s})^{-1}$ and
$M_{e}=2600~{}({\rm Pa}\cdot{\rm s})^{-1}$ associated with the velocity of
pure screw and pure edge dislocations, respectively.

To generate the training data, we run 5 large-scale DDD simulations each
starting with 12 prismatic loops seeded at different random positions.
Crystals are compressed at a strain rate of $2\times 10^{8}$/s along the [001]
direction. Simulations are run to 0.1 strain during which the dislocation
networks are regularly saved at intervals of 0.25 ps, generating a total of
10,000 configurations for training. These configurations then follow the same
procedure as the one described for DXA snapshots in §2.1: first nodal forces
are computed using DDD, and then training is performed by minimizing the STA-
based loss function in Eq. (5).

Once trained, we inspect the best resulting GNN-based mobility law by
inferring the nodal velocities of pure edge and screw dislocations as a
function of shear stress. For this, we create small infinite, straight
dislocation configurations of edge and screw characters and compute the nodal
forces $\boldsymbol{F}_{i}^{\rm DDD}=\boldsymbol{F}_{i}^{\rm app}(\tau)$,
where $\tau$ is the probing shear stress value. This simple technique allows
us to slice through the otherwise “black-box” nature of the model.

Figure 2: Edge and screw dislocation velocities predicted by the GNN-based
mobility law trained on DDD trajectories. The GNN model does capture the
linear nature of the ground-truth DDD mobility law, with linear fits of the
stress-velocity curves that give mobility coefficients in good agreement with
ground-truth values.

Results in Fig. 2 show the velocity magnitudes as function of shear stress
predicted by the trained GNN-mobility law. The linear nature of the ground-
truth DDD mobility law is well captured by the GNN model. In addition, fitting
the stress-velocity curves with a linear regression yields mobility
coefficients of $\hat{M}_{e}\approx 2564~{}({\rm Pa}\cdot{\rm s})^{-1}$ and
$\hat{M}_{s}\approx 6~{}({\rm Pa}\cdot{\rm s})^{-1}$ for the edge and screw
dislocations, respectively, which are in good agreement with the ground-truth
values. While further training could lead to yet better accuracy, this simple
example demonstrates the general feasibility of our proposed framework.

###  2.3 Application to MD trajectories

We now apply our method to learn a mobility law from large-scale MD
simulations. To illustrate the approach, we focus on the case of BCC
plasticity which is notoriously complex due to the strong plastic anisotropy
[33, 34, 35], as manifested in the tension/compression asymmetry [36, 37, 3],
and the still debated issue of slip crystallography [38]. We choose BCC W as a
testbed material because it is nearly isotropic, which makes it convenient for
the use of isotropic elasticity formulations in the DDD calculations. MD
simulations are run using the EAM-style interatomic potential developed in
[39], which yields an anisotropic ratio $A=2C_{44}/(C_{11}-C_{12})\sim 1.15$
at 300 K.

We perform 2 large-scale simulations of $\sim 35$ millions atoms in which BCC
crystals initially seeded with prismatic dislocation loops are deformed at 300
K under uniaxial tension and compression along the [001] direction. Crystals
are deformed at a true strain rate of $2\times 10^{8}$/s until reaching a
strain of 1.0. Consistent with our previous results in BCC Ta [3, 40], the
simulations predict a strong tension/compression asymmetry, with flow stress
of $\sim 4.2$ GPa in compression and $\sim 2.8$ GPa in tension. During both
runs, DXA is used to extract the dislocation configurations at every 1 ps
interval of time. This results in a total of 10,000 DXA snapshots which are
then used to train a GNN-based mobility.

Figure 3: Edge and screw dislocation velocities predicted by the GNN-based
mobility law trained on large-scale MD trajectories of BCC W deformed under
[001] tension and compression. The GNN mobility law predicts a non-linear
force-velocity relation for the edge dislocation and a strong asymmetry in the
velocity of the screw dislocation sheared in the twinning (T) and anti-
twinning (AT) direction.

Results for edge and screw dislocation velocities as a function of shear
stress predicted by the best trained GNN-based mobility law are shown in Fig.
3. The model exhibits smooth and well-behaved functions. The predicted edge
dislocation velocity shows a linear stress-velocity relation up to $\sim 300$
MPa, followed by a non-linear regime reaching an asymptotic velocity of $\sim
2000$ m/s at a stress of 2 GPa. Qualitatively, this result is fully consistent
with the known behavior of the edge dislocation [15] and previous MD
calculations in BCC metals [16, 41].

Of particular interest are the predicted velocities for the screw dislocation.
When sheared along the twinning (T) direction (corresponding to a MRSSP angle
$\chi=-30^{o}$), the screw velocity response resembles that of thermally-
activated motion, with very slow velocity at stresses below an activation
threshold (finite temperature Peierls stress), and a transition to a drag-
controlled regime at higher stresses. In contrast, when sheared in the anti-
twinning (AT) direction (corresponding to a MRSSP angle $\chi=+30^{o}$), the
screw dislocation is predicted to be much slower on average. Here again these
results are fully consistent with theoretical models [42] and earlier
observations [3].

###  2.4 Large-scale DDD+ML simulations with trained mobility

To assess the validity of our approach and the quality of the learned model in
realistic conditions, we now examine predictions of large-scale DDD+ML
simulations using the GNN-based mobility law trained in §2.3, which is
implemented within the ParaDiS DDD code [11] via the C++ PyTorch interface. To
establish a one-to-one comparison, DDD simulations are performed under
identical conditions to those used in the ground-truth MD simulations. Cubic
DDD simulation boxes of side length $296b$, where $b=0.2743$ nm is the
magnitude of the Burgers vector, are initially seeded with randomly positioned
prismatic loops (using different seeds than in the ground-truth MD runs) and
then deformed under [001] tension and compression at a strain rate of $2\times
10^{8}$/s. We use values of the shear modulus $\mu=149.78$ GPa and Poisson’s
ratio $\nu=0.289$, corresponding to values of the elastic constants of the
ground-truth interatomic potential [39] at 300 K. Dislocation core forces are
computed from core energies extracted from the same interatomic potential
using the framework detailed in [43].

Figure 4: Comparison of ground-truth MD results with DDD predictions of [001]
tension and compression deformation of BCC W using the GNN-based mobility
trained from large-scale MD data. MD and DDD simulation volumes of identical
size are deformed under identical uniaxial straining rates of $2\times
10^{8}$/s. (a) Axial stress and (b) dislocation density evolution as a
function of strain.

Results of the DDD predictions compared to MD simulations are shown in Fig. 4.
Flow stress predictions are in remarkable agreement with the MD simulations,
Fig. 4a. Specifically, the asymmetry in the tension/compression response is
naturally captured by the GNN mobility law, demonstrating that the model has
correctly learned the distinct behavior of the dislocations depending on
shearing direction. In Fig. 4b however, we note a quantitative difference in
the predicted evolution of the dislocation densities. While DDD indeed
predicts a higher dislocation density in compression compared to tension, in
agreement with the MD data, DDD also systematically underpredicts the value of
the dislocation density by a factor of approximately two. Potential causes of
these discrepancies will be discussed in the next section.

##  3 Discussion

Overall, results presented in §2.4 clearly demonstrate the feasibility and
potential of the approach. These results are all the more significant that our
previous attempts at developing traditional “hand-crafted” mobility laws to
satisfactorily capture the tension/compression asymmetry in BCC metals have
proved surprisingly challenging, with mixed results when comparing large-scale
DDD predictions to reference MD simulations. Indeed, the strong anisotropic
plasticity of BCC metals, rooted in the T/AT asymmetry [34, 44, 45, 35],
translates into a complex motion of the dislocations as a function of the
local configuration and stress state. Here, we show that such a behavior can
be naturally captured by our ML approach. Through seeing a large number of
examples during training, the GNN-based model learns to predict accurate,
effective velocities, which in our example generally result in smooth stress-
velocity relations, e.g. see Fig. 3.

Of particular interest, the GNN mobility law also encodes information about
slip crystallography (i.e. direction of the velocity vector), without pre-
existing assumption about slip planes, which is an unsettled and not fully
understood issue in BCC metals [38]. This is in contrast with earlier DDD
mobility law frameworks that assume a fixed set of glide planes and only focus
on modeling the magnitude of the dislocation velocity vector [13, 15]. This
example illustrates the novelty and power of such a tool to encode complex
dislocation physics into the DDD model, without imposing constraints (e.g.
functional forms, slip planes, etc.) that may prove detrimental.

In our first application, we nevertheless observe discrepancies compared to
the ground-truth MD simulations. While the flow stress in compression is
remarkably captured within statistical fluctuations, flow stress under tension
is slightly underestimated. Yet, as noted previously the largest discrepancy
is found in the evolution of the dislocation densities, which are roughly
underestimated by a factor two in both DDD simulations.

Several reasons could be at the origin of these discrepancies. The first set
of reasons could be related to the DDD model itself. In DDD, dislocation
mobility is not the only ingredient by which the dislocation network is
evolved. Indeed, topological operations performed to handle core reactions
(e.g. junction formation and splitting) constitute another critical step of
the DDD method [11]. This seperate step is performed in turn with the
integration of the equation of motion, and core reactions are thus absent from
the learning procedure. It is thus possible that the lack of multiplication
observed in Fig. 4b is rather a consequence of the details of the treatment of
topological operations than due to inaccuracies of the mobility law itself. An
example of this issue was recently uncovered in [46].

The second reason could be a more general one, analogous to the issue of
uncertainty quantification and transferability of interatomic potentials [47].
For instance, it is possible that – if it were feasible – large-scale _ab-
initio_ simulations would produce ground-truth outcomes that differ
significantly from MD predictions using even the most accurate interatomic
potentials fitted to it. In other words, if the emerging behavior (e.g.
crystal plasticity) is sensitive to the smallest details (e.g. details of the
electronic structure), then the transferabilty of physics across the scales
may be a daunting task. By analogy a similar difficulty may exist in coarse-
graining lattice motion to a dislocation mobility law.

As a way to address this issue, recent work has focused on the selection of
the training data [48, 49, 50, 51, 52], which has been an everlasting
challenge in the realm of interatomic potentials. As a first application of
our DDD+ML framework, here we have used training data produced on a very
narrow set of loading conditions (unique high straining rate in tension and
compression). It is thus possible that a broader training dataset could help
reduce the discrepancies. Simultaneously, it is also possible that the
applicability of our trained mobility may be limited to a narrow range of
conditions as well.

To first examine the latter issue, we have run and compared results of
simulations performed at a lower straining rate of $2\times 10^{7}$/s under
compression. To ensure that volumes remain statically representative of bulk
plasticity, the MD crystal size was increased by a factor 4 to $\sim 140$
millions atoms, i.e. corresponding to a DDD box of side length $470b$. Results
for the stress/strain predictions are shown in Fig. 5. Although the mobility
law was only trained on $2\times 10^{8}$/s trajectories, the DDD+ML model is
seen to remarkably capture the flow stress at this lower strain rate. This
result is consistent with our previous attempt at employing GNN to evolve
simplified DDD systems [26], in which we similarly observed that the GNN model
trained on high flow stress (i.e. high strain rate) simulations was able to
correctly predict the evolution of simulations at lower strain rates not seen
during training. We believe one reason for this success could be that, despite
the system being driven at a high external stress in the training data, local
dislocation lines sample a much richer set of stress conditions. Indeed, the
stress acting along dislocations lines is not just the applied stress, but the
sum of the applied stress and the internal stress arising from interactions
with all other dislocations (e.g. short-range interactions). While we leave
further aspects of this issue out the scope of this paper, future work will
focus on the accuracy and generalizability of the approach to different
conditions.

Figure 5: Comparison of ground-truth MD results with DDD predictions of [001]
compression for straining rates of $2\times 10^{8}$/s and $2\times 10^{7}$/s.
Although training was only performed on $2\times 10^{8}$/s trajectories (red),
our DDD+ML method also correctly predicts the flow stress at the lower
straining rate of $2\times 10^{7}$/s (green). Comparison of the dislocation
density is shown in the inset for the straining rate of $2\times 10^{7}$/s.

Finally, we point out that, as coarse-grained model, DDD simulations are
typically several orders of magnitude faster than MD simulations under
identical conditions. In DDD+ML, the GNN-based mobility only comes with a
marginal cost overhead compared to traditional mobility laws, and the
simulation cost is still fully dominated by force calculations. For large-
scale examples presented in this work, MD simulations (at $2\times 10^{8}$/s)
typically require $\sim 100-1000$ GPU-hours (e.g. depending on the interatomic
potential), compared to $\sim 10-20$ CPU-hours needed for the same DDD
simulations. At lower straining rates the gains are only even larger, and DDD
can access volume sizes (e.g. several $\mu$m) that are out-of-reach from
direct MD.

To summarize, we have introduced a ML framework to construct DDD mobility laws
trained on large-scale MD data using a GNN model. Training is achieved by
matching the evolution of the continuum dislocation density tensor between
subsequent DXA configurations for which an estimate of the dislocation driving
forces computed with DDD is provided as input. This scheme allows us to
circumvent both the issues of unknown nodal force and velocity vectors in the
ground-truth MD data. By applying our approach to the complex case of BCC W,
we showed that DDD+ML simulations using a trained GNN mobility are able to
reproduce well the behavior of ground-truth MD simulations it was trained on,
and that of unseen loading conditions. We believe that our proposed approach
to construct data-driven mobility laws is a promising avenue that has the
potential to significantly improve the fidelity of the DDD model while
allowing to incorporate more complex dislocation motion behaviors which, if
desired, can be made unbiased of human intervention.

##  4 Methods

###  4.1 GNN mobility law

We model the mobility law in Eq. (4) with a message-passing GNN [53], which
has proved very powerful for predicting force fields and other materials
properties [54, 55, 56, 57, 58] and for simulating complex physics [59].
Following our recent work on applying GNN to DDD simulations [26], a
dislocation configuration is represented by a graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}=\\{N_{i}\\}$ is a
collection of dislocation node/vertex features (attributes), and
$\mathcal{E}=\\{E_{ij}\\}$ is a collection of dislocation segment/edge
features. We define the input features for each node $i$ as

$N_{i}=(n_{i},\boldsymbol{F}^{\rm DDD},\|\boldsymbol{F}^{\rm DDD}\|)$ (7)

where $n_{i}$ is a flag used to specify whether node $i$ is a discretization
or physical (junction) node. The input edge features on edge $ij$ are

$E_{ij}=(e_{ij},\boldsymbol{b}_{ij},\|\boldsymbol{b}_{ij}\|,\boldsymbol{r}_{j}-\boldsymbol{r}_{i},\|\boldsymbol{r}_{j}-\boldsymbol{r}_{i}\|)$
(8)

where $e_{ij}$ is a flag used to specify whether segment $ij$ is a glissile or
junction segment, and $\boldsymbol{r}_{j}-\boldsymbol{r}_{i}$ are the local
segments line vectors, naturally compatible with the use of periodic boundary
conditions. To satisfy Burgers vector conservation, we use a directed graph,
i.e. if $ij$ is an edge with Burgers vector $\boldsymbol{b}_{ij}$, then $ji$
is also an edge but with opposite Burgers vector $-\boldsymbol{b}_{ij}$.

Our GNN architecture follows [60] and is first composed of vertex and edge
encoders $\text{ENC}^{V}$, $\text{ENC}^{E}$ transforming concatenated input
features into a latent space

$v^{(0)}_{i}=\text{ENC}^{V}(N_{i}),\ e^{(0)}_{ij}=\text{ENC}^{E}(E_{ij}),$ (9)

followed by $K$ stacked message passing layers $f^{E(k)}$, $f^{V(k)}$ ($1\leq
k\leq K$) sequentially updating the latent vertex and edge variables

$\displaystyle e^{(k)}_{ij}$
$\displaystyle=f^{E(k)}(e^{(k-1)}_{ij},v^{(k-1)}_{i},v^{(k-1)}_{j}),$ (10)
$\displaystyle v^{(k)}_{i}$
$\displaystyle=f^{V(k)}(v^{(k-1)}_{i},\sum_{j}e^{(k)}_{ij}),$ (11)

and finally a node decoder DEC that translates the latent node variables
$v^{(K)}$ into the desired output properties, i.e. nodal velocity vectors:

$\boldsymbol{V}_{i}=\text{DEC}(v^{(K)}_{i}).$ (12)

Functions $\text{ENC}^{V}$, $\text{ENC}^{E}$, $f^{V}$, $f^{E}$, and DEC are
neural network operators built from multi-layer perceptrons with two hidden
layers, layer normalization [61], skip connections [62], and GELU activation
functions [63].

###  4.2 MD simulations

Large-scale MD simulations of BCC tungsten are performed with LAMMPS [64]
under the Kokkos GPU implementation [65], using the EAM-style interatomic
potential developed in [39]. Simulations are performed following the protocol
introduced in [1]. Periodic, orthorombic BCC perfect crystals are initially
seeded with twelve $1/2\langle 111\rangle\\{110\\}$ prismatic loops of the
vacancy type. The crystals are first equilibrated at the temperature of 300 K,
after which they are deformed at a constant true strain rate. Temperature and
uniaxial loading conditions are maintained during deformation using the
_langevin_ thermostat and the _nph_ barostat. For the $\sim 35$ millions atoms
simulations deformed at a rate of $2\times 10^{8}$/s, DXA [28] is executed
every 1 ps to save the detailed evolution of the dislocation networks.

###  4.3 Mobility workflow and training

DXA configurations produced in the MD simulations are first converted to the
ParaDiS format [11]. For consistency, during this operation the line networks
are also remeshed with discretization size of $\sim 10b$, corresponding to the
average segment size in used our DDD simulations.

The so-converted dislocation configurations $\\{\mathcal{G}^{t_{s}}\\}$ are
then fed to the ParaDiS code to compute nodal forces, Eq. (2). Applied forces
$\boldsymbol{F}_{i}^{app}$ are computed by integrating the Peach-Koehler force
$(\boldsymbol{\sigma}^{t_{s}}\cdot\boldsymbol{b}_{ij})\times\boldsymbol{t}_{ij}$
along the dislocation segments $ij$ with unit tangent $\boldsymbol{t}_{ij}$,
where $\boldsymbol{\sigma}^{t_{s}}$ is the instantaneous stress applied to
network $\mathcal{G}^{t_{s}}$ at time $t_{s}$ as recorded during the MD runs.
Long-range interaction forces $\boldsymbol{F}_{i}^{lr}$ are computed using
DDD-FFT approach introduced in [30] which can easily handle non-cubic,
deforming simulation boxes as produced by MD simulations. Short-range
interaction forces $\boldsymbol{F}_{i}^{sr}$ are computed using the non-
singular isotropic analytical formulation [66]. Core forces
$\boldsymbol{F}_{i}^{core}$ are computed from core energies extracted from the
ground-truth interatomic potential [39] using the framework developed in [43].

The networks containing nodal forces are then used as inputs to our training
procedure implemented within PyTorch [31]. To facilitate training, input
forces are rescaled by a factor $10^{9}~{}{\rm Pa}\cdot b^{2}$ so that their
average magnitude is on the order of unity. The loss function, Eq. (5), is
computed by evaluating the Nye’s tensor on a grid of $32^{3}$ voxels using a
fully vectorized implementation of the discrete-to-continuous method
introduced in [30].

We trained different GNN models to explore different sets of hyper-parameters.
We tested a combination of models with $K=\\{2,3\\}$ message-passing layers
each with latent space of size $L=\\{48,96\\}$, leading to 4 different trained
models. The models were trained using the Adamw optimizer with weight decay of
$1\times 10^{-5}$ [67]. Training was performed for 12 hours with batch size of
4 on a single NVidia V100 GPU. We find that the GNN model with $K=3$ and
$L=48$ offers the best trade-off between accuracy and complexity (78,768 total
parameters) while showing no sign of ovefitting. We thus selected it as the
best model for results presented in this work.

## Acknowledgement

NB and VB acknowledge support by the Laboratory Directed Research and
Development (LDRD) program (22-ERD-016) and by the ASC PEM program at Lawrence
Livermore National Laboratory (LLNL). FZ was supported by the Critical
Materials Institute, an Energy Innovation Hub funded by the U.S. Department of
Energy, Office of Energy Efficiency and Renewable Energy, and Advanced
Materials and Manufacturing Technologies Office. Computing support for this
work came from LLNL Institutional Computing Grand Challenge program. This work
was performed under the auspices of the U.S. Department of Energy by LLNL
under contract DE-AC52-07NA27344.

## References

  * [1] L. A. Zepeda-Ruiz, A. Stukowski, T. Oppelstrup, V. V. Bulatov, Probing the limits of metal plasticity with molecular dynamics simulations, Nature 550 (7677) (2017) 492–495.
  * [2] L. A. Zepeda-Ruiz, A. Stukowski, T. Oppelstrup, N. Bertin, N. R. Barton, R. Freitas, V. V. Bulatov, Atomistic insights into metal hardening, Nature materials 20 (3) (2021) 315–320.
  * [3] N. Bertin, L. Zepeda-Ruiz, V. Bulatov, Sweep-tracing algorithm: in silico slip crystallography and tension-compression asymmetry in bcc metals, Materials Theory 6 (1) (2022) 1–23.
  * [4] J. C. Stimac, N. Bertin, J. K. Mason, V. V. Bulatov, Energy storage under high-rate compression of single crystal tantalum, Acta Materialia 239 (2022) 118253\.
  * [5] L. P. Kubin, G. Canova, M. Condat, B. Devincre, V. Pontikis, Y. Bréchet, Dislocation microstructures and plastic flow: a 3d simulation, Solid state phenomena 23 (1992) 455–472.
  * [6] H. M. Zbib, M. Rhee, J. P. Hirth, On plastic deformation and the dynamics of 3d dislocations, International Journal of Mechanical Sciences 40 (2-3) (1998) 113–127.
  * [7] V. Bulatov, F. F. Abraham, L. Kubin, B. Devincre, S. Yip, Connecting atomistic and mesoscale simulations of crystal plasticity, Nature 391 (6668) (1998) 669–672.
  * [8] K. Schwarz, Simulation of dislocations on the mesoscopic scale. i. methods and examples, Journal of Applied Physics 85 (1) (1999) 108–119.
  * [9] N. Ghoniem, M, S.-H. Tong, L. Sun, Parametric dislocation dynamics: a thermodynamics-based approach to investigations of mesoscopic plastic deformation, Physical Review B 61 (2) (2000) 913.
  * [10] D. Weygand, L. Friedman, E. Van der Giessen, A. Needleman, Aspects of boundary-value problem solutions with three-dimensional dislocation dynamics, Modelling and Simulation in Materials Science and Engineering 10 (4) (2002) 437\.
  * [11] A. Arsenlis, W. Cai, M. Tang, M. Rhee, T. Oppelstrup, G. Hommes, T. G. Pierce, V. V. Bulatov, Enabling strain hardening simulations with dislocation dynamics, Modelling and Simulation in Materials Science and Engineering 15 (6) (2007) 553–595.
  * [12] V. Bulatov, W. Cai, Computer simulations of dislocations, Vol. 3, Oxford University Press on Demand, 2006.
  * [13] Z. Wang, I. Beyerlein, An atomistically-informed dislocation dynamics model for the plastic anisotropy and tension–compression asymmetry of bcc metals, International Journal of Plasticity 27 (10) (2011) 1471–1484.
  * [14] K. Srivastava, R. Gröger, D. Weygand, P. Gumbsch, Dislocation motion in tungsten: atomistic input to discrete dislocation simulations, International Journal of Plasticity 47 (2013) 126–142.
  * [15] G. Po, Y. Cui, D. Rivera, D. Cereceda, T. D. Swinburne, J. Marian, N. Ghoniem, A phenomenological dislocation mobility law for bcc metals, Acta Materialia 119 (2016) 123–135.
  * [16] J. Chang, W. Cai, V. V. Bulatov, S. Yip, Dislocation motion in bcc metals by molecular dynamics, Materials Science and Engineering: A 309 (2001) 160–163.
  * [17] D. L. Olmsted, L. G. Hector, W. Curtin, R. Clifton, Atomistic simulations of dislocation mobility in al, ni and al/mg alloys, Modelling and Simulation in Materials Science and Engineering 13 (3) (2005) 371.
  * [18] S. Queyreau, J. Marian, M. Gilbert, B. Wirth, Edge dislocation mobilities in bcc fe obtained by molecular dynamics, Physical Review B 84 (6) (2011) 064106\.
  * [19] D. Cereceda, A. Stukowski, M. Gilbert, S. Queyreau, L. Ventelon, M.-C. Marinica, J. Perlado, J. Marian, Assessment of interatomic potentials for atomistic analysis of static and dynamic properties of screw dislocations in w, Journal of Physics: Condensed Matter 25 (8) (2013) 085702.
  * [20] J. Cho, J.-F. Molinari, G. Anciaux, Mobility law of dislocations with several character angles and temperatures in fcc aluminum, International Journal of Plasticity 90 (2017) 66–75.
  * [21] G. Monnet, D. Terentyev, Structure and mobility of the 12¡ 111¿$\\{$112$\\}$ edge dislocation in bcc iron studied by molecular dynamics, Acta Materialia 57 (5) (2009) 1416–1426.
  * [22] K. Kang, V. V. Bulatov, W. Cai, Singular orientations and faceted motion of dislocations in body-centered cubic crystals, Proceedings of the National Academy of Sciences 109 (38) (2012) 15174–15178.
  * [23] P.-A. Geslin, D. Rodney, Thermal fluctuations of dislocations reveal the interplay between their core energy and long-range elasticity, Physical Review B 98 (17) (2018) 174115.
  * [24] N. Bertin, R. B. Sills, W. Cai, Frontiers in the simulation of dislocations, Annual Review of Materials Research 50 (2020) 437–464.
  * [25] R. B. Sills, N. Bertin, A. Aghaei, W. Cai, Dislocation networks and the microstructural origin of strain hardening, Physical review letters 121 (8) (2018) 085501.
  * [26] N. Bertin, F. Zhou, Accelerating discrete dislocation dynamics simulations with graph neural networks, Journal of Computational Physics 487 (2023) 112180.
  * [27] A. Stukowski, K. Albe, Extracting dislocations and non-dislocation crystal defects from atomistic simulation data, Modelling and Simulation in Materials Science and Engineering 18 (8) (2010) 085001.
  * [28] A. Stukowski, A triangulation-based method to identify dislocations in atomistic models, Journal of the Mechanics and Physics of Solids 70 (2014) 314–319.
  * [29] J. Nye, Some geometrical relations in dislocated crystals, Acta metallurgica 1 (2) (1953) 153–162.
  * [30] N. Bertin, Connecting discrete and continuum dislocation mechanics: A non-singular spectral framework, International Journal of Plasticity 122 (2019) 268–284.
  * [31] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, S. Chintala, Pytorch: An imperative style, high-performance deep learning library, in: Advances in Neural Information Processing Systems 32, Curran Associates, Inc., 2019, pp. 8024–8035.
  * [32] W. Cai, V. V. Bulatov, Mobility laws in dislocation dynamics simulations, Materials Science and Engineering: A 387 (2004) 277–281.
  * [33] J. Christian, Some surprising features of the plastic deformation of body-centered cubic metals and alloys, Metallurgical Transactions A 14 (7) (1983) 1237–1256.
  * [34] M. a.-S. Duesbery, V. Vitek, Plastic anisotropy in bcc transition metals, Acta Materialia 46 (5) (1998) 1481–1492.
  * [35] L. Dezerald, D. Rodney, E. Clouet, L. Ventelon, F. Willaime, Plastic anisotropy and dislocation trajectory in bcc metals, Nature Communications 7 (2016) 11695\.
  * [36] P. Sherwood, F. Guiu, H. C. Kim, P. L. Pratt, Plastic anisotropy of tantalum, niobium, and molybdenum, Canadian Journal of Physics 45 (2) (1967) 1075–1089.
  * [37] G. L. Webb, R. Gibala, T. E. Mitchell, Effect of normal stress on yield asymmetry in high purity tantalum crystals, Metallurgical Transactions 5 (7) (1974) 1581–1584.
  * [38] C. R. Weinberger, B. L. Boyce, C. C. Battaile, Slip planes in bcc transition metals, International Materials Reviews 58 (5) (2013) 296–314.
  * [39] N. Juslin, B. Wirth, Interatomic potentials for simulation of he bubble formation in w, Journal of nuclear materials 432 (1-3) (2013) 61–66.
  * [40] N. Bertin, R. Carson, V. V. Bulatov, J. Lind, M. Nelms, Crystal plasticity model of bcc metals from large-scale md simulations, Acta Materialia 260 (2023) 119336.
  * [41] Y. N. Osetsky, D. J. Bacon, An atomic-level model for studying the dynamics of edge dislocations in metals, Modelling and simulation in materials science and engineering 11 (4) (2003) 427.
  * [42] K. Edagawa, T. Suzuki, S. Takeuchi, Motion of a screw dislocation in a two-dimensional peierls potential, Physical Review B 55 (10) (1997) 6180.
  * [43] N. Bertin, W. Cai, S. Aubry, V. Bulatov, Core energies of dislocations in bcc metals, Physical Review Materials 5 (2) (2021) 025002.
  * [44] K. Ito, V. Vitek, Atomistic study of non-schmid effects in the plastic yielding of bcc metals, Philosophical Magazine A 81 (5) (2001) 1387–1407.
  * [45] V. Vitek, Core structure of screw dislocations in body-centred cubic metals: relation to symmetry and interatomic bonding, Philosophical Magazine 84 (3-5) (2004) 415–428.
  * [46] N. Bertin, W. Cai, S. Aubry, A. Arsenlis, V. V. Bulatov, Enhanced mobility of dislocation network nodes and its effect on dislocation multiplication and strain hardening, arXiv preprint arXiv:2210.14343 (2022).
  * [47] A. Chernatynskiy, S. R. Phillpot, R. LeSar, Uncertainty quantification in multiscale simulation of materials: A prospective, Annual Review of Materials Research 43 (2013) 157–182.
  * [48] W. Jeong, K. Lee, D. Yoo, D. Lee, S. Han, Toward reliable and transferable machine learning potentials: uniform training by overcoming sampling bias, The Journal of Physical Chemistry C 122 (39) (2018) 22790–22795.
  * [49] G. Imbalzano, A. Anelli, D. Giofré, S. Klees, J. Behler, M. Ceriotti, Automatic selection of atomic fingerprints and reference configurations for machine-learning potentials, The Journal of chemical physics 148 (24) (2018).
  * [50] E. V. Podryabinkin, E. V. Tikhonov, A. V. Shapeev, A. R. Oganov, Accelerating crystal structure prediction by machine-learning interatomic potentials with active learning, Physical Review B 99 (6) (2019) 064114.
  * [51] N. Bernstein, G. Csányi, V. L. Deringer, De novo exploration and self-guided learning of potential-energy surfaces, npj Computational Materials 5 (1) (2019) 99.
  * [52] D. Montes de Oca Zapiain, M. A. Wood, N. Lubbers, C. Z. Pereyra, A. P. Thompson, D. Perez, Training data selection for accuracy and transferability of interatomic potentials, npj Computational Materials 8 (1) (2022) 189.
  * [53] P. W. Battaglia, J. B. Hamrick, V. Bapst, A. Sanchez-Gonzalez, V. Zambaldi, M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, et al., Relational inductive biases, deep learning, and graph networks (2018).  arXiv:1806.01261.
  * [54] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, G. E. Dahl, Neural message passing for quantum chemistry, in: International conference on machine learning, PMLR, 2017, pp. 1263–1272.
  * [55] T. Xie, J. C. Grossman, Crystal Graph Convolutional Neural Networks for an Accurate and Interpretable Prediction of Material Properties, Physical Review Letters 120 (14) (2018) 145301.
  * [56] C. Chen, W. Ye, Y. Zuo, C. Zheng, S. P. Ong, Graph Networks as a Universal Machine Learning Framework for Molecules and Crystals, Chemistry of Materials 31 (2018) 3564.
  * [57] C. W. Park, M. Kornbluth, J. Vandermause, C. Wolverton, B. Kozinsky, J. P. Mailoa, Accurate and scalable graph neural network force field and molecular dynamics with direct force architecture, npj Computational Materials 7 (1) (2021) 73.
  * [58] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, B. Kozinsky, E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials, Nature Communications 13 (1) (2022) 2453.
  * [59] A. Sanchez-Gonzalez, J. Godwin, T. Pfaff, R. Ying, J. Leskovec, P. Battaglia, Learning to simulate complex physics with graph networks, in: International Conference on Machine Learning, PMLR, 2020, pp. 8459–8468.
  * [60] T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, P. W. Battaglia, Learning mesh-based simulation with graph networks (2020).  arXiv:2010.03409.
  * [61] J. L. Ba, J. R. Kiros, G. E. Hinton, Layer normalization (2016).  arXiv:1607.06450.
  * [62] K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition (2015).  arXiv:1512.03385.
  * [63] D. Hendrycks, K. Gimpel, Gaussian error linear units (gelus) (2016).  arXiv:1606.08415.
  * [64] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, et al., Lammps-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Computer Physics Communications 271 (2022) 108171.
  * [65] H. C. Edwards, C. R. Trott, D. Sunderland, Kokkos: Enabling manycore performance portability through polymorphic memory access patterns, Journal of parallel and distributed computing 74 (12) (2014) 3202–3216.
  * [66] W. Cai, A. Arsenlis, C. R. Weinberger, V. V. Bulatov, A non-singular continuum theory of dislocations, Journal of the Mechanics and Physics of Solids 54 (3) (2006) 561–587.
  * [67] I. Loshchilov, F. Hutter, Decoupled weight decay regularization (2017).  arXiv:1711.05101.


* (20) H. Tashiro, W. Chen, F. Ferrer and T. Vachaspati, Mon. Not. R. Astron. Soc. 445, L41(2014), arXiv:1310.4826 [astro-ph.CO].
  * (21) A. M. Taylor, I. Vovk and A. Neronov, Astron. Astrophys. 529, A144 (2011), arXiv:1101.0932 [astro-ph.HE].
  * (22) K. Dolag, M. Kachelriess, S. Ostapchenko and R. Tomas, Astrophys. J. Lett. 727, L4 (2011), arXiv:1009.1782 [astro-ph.HE].
  * (23) A. Neronov and I. Vovk, Science 328, 73 (2010), arXiv:1006.3504 [astro-ph.HE].
  * (24) F. Tavecchio, G. Ghisellini, L. Foschini, G. Bonnoli, G. Ghirlanda and P. Coppi, Mon. Not. Roy. Astron. Soc. 406, L70 (2010), arXiv:1004.1329 [astro-ph.CO].
  * (25) V. Skokov, A. Y. Illarionov, and V. Toneev, Int. J. Mod. Phys. A 24, 5925 (2009), arXiv:0907.1396 [nucl-th].
  * (26) W.-T. Deng and X.-G. Huang, Phys. Rev. C 85, 044907 (2012), arXiv:1201.5108 [nucl-th].
  * (27) D. E. Kharzeev, J. Liao, S. A. Voloshin and G. Wang, Prog. Part. Nucl. Phys. 88, 1 (2016), arXiv:1511.04050 [hep-ph].
  * (28) V. A. Miransky and I. A. Shovkovy, Phys. Rep. 576, 1 (2015), arXiv:1503.00732 [hep-ph].
  * (29) A. Bandyopadhyay, C. A. Islam, and M. G. Mustafa, Phys. Rev. D 94, 114034 (2016), arXiv:1602.06769 [hep-ph].
  * (30) A. Das, N. Haque, M. G. Mustafa, and P. K. Roy, Phys. Rev. D 99, 094022 (2019), arXiv:1903.03528 [hep-ph].
  * (31) R. Ghosh, B. Karmakar, and M. G. Mustafa, Phys. Rev. D 101, 056007 (2020), arXiv:1911.00744 [hep-ph].
  * (32) N. Sadooghi and F. Taghinavaz, Ann. Phys. (Amsterdam) 376, 218 (2017), arXiv:1601.04887 [hep-ph].
  * (33) S. Ghosh and V. Chandra, Phys. Rev. D 98, 076006 (2018), arXiv:1808.05176 [hep-ph].
  * (34) S. Ghosh, N. Chaudhuri, S. Sarkar, and P. Roy, Phys. Rev. D 101, 096002 (2020), arXiv:2004.09203 [nucl-th].
  * (35) The absorptive part of the photon polarization function includes the imaginary part of the symmetric tensor structure and the real part of the antisymmetric tensor structure. Note, however, that the latter vanishes at $\mu=0$ because of the charge conjugation symmetry.
  * (36) A. A. Sokolov and I. M. Ternov, Radiation from Relativistic Electrons, (American Institute of Physics, New York, 1986).
  * (37) X. Wang, I. A. Shovkovy, L. Yu, and M. Huang, Phys. Rev. D 102, 076010 (2020), arXiv:2006.16254 [hep-ph].
  * (38) K. Tuchin, Phys. Rev. C 87, 024912 (2013), arXiv:1206.0485 [hep-ph].
  * (39) G. Basar, D. E. Kharzeev, and V. Skokov, Phys. Rev. Lett. 109, 202303 (2012), arXiv:1206.1334 [hep-ph].
  * (40) H. U. Yee, Phys. Rev. D 88, 026001 (2013), [arXiv:1303.3571 [nucl-th]].
  * (41) K. Tuchin, Phys. Rev. C 91, 014902 (2015), [arXiv:1406.5097 [nucl-th]].
  * (42) B. G. Zakharov, Eur. Phys. J. C 76, 609 (2016), [arXiv:1609.04324 [nucl-th]].
  * (43) I. S. Gradshtein and I. M. Ryzhik, Tables of Integrals, Series, and Products (Academic Press, Orlando, 1980).
  * (44) J. Kapusta and C. Gale, in Finite-temperature field theory: Principles and applications, (Cambridge University Press, Cambridge, England, 2006), p. 330.
  * (45) X. Wang and I. Shovkovy, arXiv:2106.09029 [nucl-th].
  * (46) P. E. C. Ashby and J. P. Carbotte, Phys. Rev. B 87, 245131 (2013), arXiv:1305.0275 [cond-mat.mes-hall].
  * (47) Z. Long, Y. Wang, M. Erukhimova, M. Tokman, and A. Belyanin, Phys. Rev. Lett. 120, 037403 (2018), arXiv:1708.05498 [cond-mat.mes-hall].
  * (48) M. Stålhammar, J. Larana-Aragon, J. Knolle, and E. J. Bergholtz, Phys. Rev. B 102, 235134 (2020), arXiv:2009.07018 [cond-mat.mes-hall].
  * (49) R. Y. Chen, Z. G. Chen, X.-Y. Song, J. A. Schneeloch, G. D. Gu, F. Wang, and N. L. Wang, Phys. Rev. Lett. 115, 176404 (2015), arXiv:arXiv:1506.08676 [cond-mat.mes-hall].
  * (50) M. Abramowitz and I. A. Stegun, _Handbook of mathematical functions: with formulas, graphs, and mathematical tables._ (Dover Publications, New York, 1964).