# Deep-Learned Event Variables for Collider Phenomenology

Doojin Kim<EMAIL_ADDRESS>Mitchell Institute for Fundamental Physics and
Astronomy, Department of Physics and Astronomy, Texas A&M University, College
Station, TX 77843, USA  Kyoungchul Kong<EMAIL_ADDRESS>Department of Physics
and Astronomy, University of Kansas, Lawrence, KS 66045, USA  Konstantin T.
Matchev<EMAIL_ADDRESS>Institute for Fundamental Theory, Physics Department,
University of Florida, Gainesville, FL 32611, USA  Myeonghun Park
<EMAIL_ADDRESS>Institute of Convergence Fundamental Studies,
Seoultech, Seoul, 01811, Korea  School of Physics, KIAS, Seoul 02455, Korea
Center for Theoretical Physics of the Universe, Institute for Basic Science,
Daejeon 34126 Korea  Prasanth Shyamsundar<EMAIL_ADDRESS>Fermilab Quantum
Institute, Fermi National Accelerator Laboratory, Batavia, IL 60510, USA

###### Abstract

The choice of optimal event variables is crucial for achieving the maximal
sensitivity of experimental analyses. Over time, physicists have derived
suitable kinematic variables for many typical event topologies in collider
physics. Here we introduce a deep learning technique to design good event
variables, which are sensitive over a wide range of values for the unknown
model parameters. We demonstrate that the neural networks trained with our
technique on some simple event topologies are able to reproduce standard event
variables like invariant mass, transverse mass, and stransverse mass. The
method is automatable, completely general, and can be used to derive
sensitive, previously unknown, event variables for other, more complex event
topologies.

††preprint: FERMILAB-PUB-21-244-QIS††preprint: MI-TH-2111

Introduction. Data in collider physics is very high-dimensional, which brings
a number of challenges for the analysis, encapsulated in “the curse of
dimensionality” [1]. Mapping the raw data to reconstructed objects involves
initial dimensionality reduction in several stages, including track
reconstruction, calorimeter clustering, jet reconstruction, etc. Subsequently,
the kinematics of the reconstructed objects is used to define suitable
analysis variables, adapted to the specific channel and targeted event
topology. Each such step is essentially a human-engineered feature-extraction
process from complicated data to a handful of physically meaningful
quantities. While some information loss is unavoidable, physics principles and
symmetries help keep it to a minimum.

In this letter, we shall focus on the last stage of this dimensionality
reduction chain, namely, the optimal construction of kinematic variables,
which is essential to expedite the discovery of new physics and/or to improve
the precision of parameter measurements. By now, the experimentalist’s toolbox
contains a large number of kinematic variables, which have been thoroughly
tested in analyses with real data (see [2, 3, 4, 5] for reviews). The latest
important addition to this set are the so-called “singularity variables” [6,
7, 8, 9, 10], which are applicable to missing energy events — the harbingers
of dark matter production at colliders. In the machine learning era, a myriad
of algorithms have been invented or adopted to tackle various tasks that arise
in the analysis of collider data, e.g., signal–background discrimination (see
[11] for a continuously updated complete review of the literature). _Under the
hood_ , the machines trained in these techniques could learn to construct
useful features from the low-level event description, because they are
relevant to the task at hand. But it is difficult to interpret what exactly
the machines have learned in the process [12, 13]. Furthermore, it is rarely
studied whether the human-engineered features are indeed the best event
variables for certain purposes, and whether machines can outperform theorists
at constructing event variables.

These two issues, explainability and optimality, are precisely the two
questions which we shall address in this letter. We shall introduce a new
technique for training neural networks to directly _output_ useful features or
event variables (which offer sensitivity over a range of unknown parameter
values). This allows for explainability of the machine’s output by comparison
against known features in the data. At the same time, it is important to
verify that the variables obtained using our technique are indeed the optimal
choice, and we will test this by directly comparing them against the human-
engineered variables that are known to be optimal for their respective event
topologies. Once we have validated our training procedure in this way, we
could extend it to more complex event topologies and derive novel kinematic
variables in interesting and difficult scenarios.

Understanding how and what a neural network (NN) learns is a difficult task.
Here we shall consider relatively simple physics examples that are
nevertheless highly non-trivial from a machine learning point of view: (1)
visible two-body decay (to two visible daughter particles); (2) semi-invisible
two-body decay (to one visible and one invisible daughter particle); (3) semi-
invisible two-body decays of pair-produced particles. It is known that the
relevant variables in those situations are the invariant mass $m$, the
transverse mass $m_{T}$ [14, 15] and the stransverse mass $m_{T2}$ [16],
respectively. We will demonstrate that in each case, the NN can be trained to
learn the corresponding physics variable in the reduced latent space. The
method can be readily generalized to more complex cases to derive deep-
learned, next-generation event variables.

Methodology. Let $X$ represent the high-dimensional input features from a
collision event, e.g., the 4-momenta of the reconstructed physics objects. Let
$V(X)$ be a low-dimensional event variable constructed from $X$. In this work,
we shall model the function $V$ using a neural network, where for notational
convenience, the dependence of $V$ on the architecture and weights of the
network will not be explicitly indicated. We imagine that $V$ retains the
relevant physics information and will be the centerpiece of an experimental
study of a theory model with a set of unknown parameters $\Theta$. The goal is
to train the NN encoding the function $V$ to be “useful” over a wide range of
values for $\Theta$. For this purpose, we will need to train with events
generated from a range of $\Theta$ values. Note that this is a departure from
the traditional approach in particle physics, where training is done for
specific study points with fixed values of $\Theta$. In addition, we will have
to quantify the usefulness of a given event variable $V(X)$, as explained in
the remainder of this section.

Intuition from Information Theory. Each event $X$ carries some information
about the underlying model parameter values from which it was produced. Some
of this information could be lost when reducing the dimensionality of the data
from $d_{X}\equiv\dim(X)$ to $d_{V}\equiv\dim(V)$, as a consequence of the
data processing inequality [17]. Good event variables minimize this
information loss, and efficiently retain the information about the underlying
parameter values $\Theta$ [18, 19]. This is precisely why the invariant mass
$m$, the transverse mass $m_{T}$ and the stransverse mass $m_{T2}$ have been
widely used in particle physics for mass parameter measurements and for new
physics searches.

The mutual information of $V$ and $\Theta$ is given by

$I(V;\Theta)\equiv\int\limits_{\mathcal{V}}dv\int\limits_{\Omega}d\theta~{}p_{(V,\Theta)}(v,\theta)\,\ln\\!{\left[\frac{p_{(V,\Theta)}(v,\theta)}{p_{V}(v)\,p_{\Theta}(\theta)}\right]},$
(1)

where $p_{V}$ and $p_{\Theta}$ are the probability distribution functions of
$V$ and $\Theta$, respectively, and $p_{(V,\Theta)}$ is their joint
distribution. $\mathcal{V}$ and $\Omega$ are the domains of $V$ and $\Theta$,
respectively. One can think of $p_{\Theta}$ as the prior distribution of
$\Theta$.111For convenience, we will adopt the Bayesian interpretation of
probability in the presentation of this work. The distributions
$p_{(V,\Theta)}$ and $p_{V}$ can then be derived from $p_{\Theta}$ and the
conditional distribution $p_{V|\Theta}(v|\theta)$.

The mutual information $I(V;\Theta)$ quantifies the amount of information
contained in $V$ about $\Theta$. Therefore, a good event variable $V$ should
have relatively high values of $I(V;\Theta)$. From Eq. (1), one can see that
$I(V;\Theta)$ is nothing but the Kullback–Leibler (KL) divergence from (a) the
factorized distribution $p_{V}\otimes p_{\Theta}$ to (b) the joint
distribution $p_{(V,\Theta)}$. The KL divergence, in turn, is a measure of how
distinguishable the two distributions (a) and (b) are.

These observations lead to the following strategy schematically outlined in
FIG. 1: train the event variable network so that the distributions
$p_{V}\otimes p_{\Theta}$ and $p_{(V,\Theta)}$ are highly distinguishable. An
auxiliary classifier network can then be used for evaluating the
distinguishability of the two distributions. The basic blueprint of our
training technique will be described next.

Training Data Generation. In order to generate the training data, we start
with the two distributions $p_{\Theta}$ and $p_{X|\Theta}$. The specific
choice of a prior distribution $p_{\Theta}$ is not crucial — as long as it
allows to sample $\theta$ over a sufficiently wide range (the one in which we
want the event variable $V$ to be sensitive) any function will do, and one is
further free to impose theoretical prejudice like fine tuning, etc.

$p_{X|\Theta}$ is the distribution of the event $X$ conditional on $\Theta$.
General purpose event generators can be used to sample from this distribution.
The overall distribution of $X$, namely $p_{X}$, is given by

$p_{X}(x)=\int\limits_{\Omega}d\theta~{}p_{\Theta}(\theta)~{}p_{X|\Theta}(x|\theta).$
(2)

Our training data consists of two classes, whose generation is illustrated in
the left (green) block of FIG. 1. Each training data point is given by a
2-tuple $(X,\Theta)$ along with the class label
$y_{\mathrm{target}}\in\\{0,1\\}$ of the data point. Under class 0, $X$ and
$\Theta$ are independent of each other and their joint distribution is given
by $p_{X}\otimes p_{\Theta}$. This is accomplished by simply replacing the
true value of $\Theta$ used to generate $X$ with a fake one for the datapoints
in class 0. Under class 1, the joint distribution of $(X,\Theta)$ is given by

$p_{(X,\Theta)}(x,\theta)=p_{X|\Theta}(x|\theta)~{}p_{\Theta}(\theta).$ (3)

Figure 1: A schematic diagram of the training strategy for the artificial
event variable $V$. The left (green) and right (blue) blocks depict the
generation of training data and the composite neural network layout,
respectively.

Event Variable Training. As shown in the right (blue) block FIG. 1, we then
set up a composite network for classifying the data points $(X,\Theta)$ into
the two classes. The composite network consists of two parts. First, an event
variable network (EVN) takes the high-dimensional collider event information
$X$ as input and returns a low-dimensional $V(X)$ as output. As indicated,
this network parameterizes the artificial event variable function $V(X)$,
which is precisely what we are interested in training. The output layer of the
EVN network does not use an activation function (or, equivalently, uses the
identity activation). Since $d_{V}\ll d_{X}$, the main task of the EVN network
is to perform the needed dimensionality reduction. However, to ensure that
this retains the maximal amount of information, we introduce an auxiliary
classifier network which takes the event variable $V(X)$ and the parameters
$\Theta$ as input and returns a one-dimensional output, $y(V,\Theta)\in[0,1]$.
Note that the input received by the auxiliary network is distributed as
$p_{V}\otimes p_{\Theta}$ under class 0, and as $p_{(V,\Theta)}$ under class
1.

The information bottleneck [18] $V(X)$ created by the EVN module is optimized
by simply training the composite network as a classifier for the input data
$(X,\Theta)$, using the class labels $y_{\mathrm{target}}$ as the supervisory
signal.

Experiments. The EVN module in the network architecture from FIG. 1 reduces
the original $d_{X}$-dimensional features to a $d_{V}$-dimensional subspace of
event variables, which by construction are guaranteed to be highly sensitive
to the theory model parameters $\Theta$, but without any explicit dependence
on them. Such variables have been greatly valued in collider phenomenology,
and a significant number have been proposed and used in experimental analyses.
As a proof of principle, we shall now demonstrate how our approach is able to
reproduce the known kinematic variables in a few simple but non-trivial
examples. Here we shall only consider one variable at a time, i.e., $d_{V}=1$,
postponing the case of $d_{V}>1$ to future work [20].

Example 1: Fully visible two-body decay. First we consider the fully visible
decay of a parent particle $A$ into two massless visible daughter particles,
$A\rightarrow b\,c$. The parameter $\Theta$ in this example is the mass
$m_{A}$ of the mother particle $A$. The event $X$ is specified by the
4-momenta of the daughter particles $p_{b}$ and $p_{c}$, leading to $d_{X}=8$.

The prior $p_{\Theta}$ is chosen to sample $m_{A}$ uniformly in the range
$[100,500]$ GeV. For each sampled value of $m_{A}$, we generate an event as
follows. A generic boost for the parent particle $A$ is obtained by
isotropically picking the direction for its momentum and uniformly sampling
its lab frame energy in the range $[m_{A},1500\text{ GeV}]$. Subsequently, $A$
is decayed on-shell into two massless particles (isotropically in its own rest
frame), so that the input data consists of the lab-frame final state 4-momenta
$\\{p_{b},p_{c}\\}\equiv\\{E_{b},\vec{p}_{b},E_{c},\vec{p}_{c}\\}$.

All the neural networks used in this work were implemented in TensorFlow [21].
For the event variable network, we used a sequential fully connected
architecture with 5 hidden layers. The hidden layers, in order, have 128, 64,
64, 64, and 32 nodes, all of them using ReLU as the activation function. The
output layer has one node with no activation function. The classifier network
is a fully connected network with 3 hidden layers (16 nodes each, ReLU
activation). The output layer has one node with sigmoid activation. These two
networks were combined as shown in the right (blue) block of FIG. 1 and
trained with 2.5 million events total (50–50 split between classes 0 and 1),
out of which 20% was set aside for validation. The network was trained for 20
epochs with a mini-batch size of 50, using the Adam optimizer and the binary
crossentropy loss function.

For the event topology considered in this example, it is known that the event
variable most sensitive to the value of $m_{A}$ is the invariant mass of the
daughter particles

$m_{bc}=\sqrt{(E_{b}+E_{c})^{2}-(\vec{p}_{b}+\vec{p}_{c})^{2}~{}}\,,$ (4)

as well as any variable that is in one-to-one correspondence with it. In order
to test whether our artificial event variable $V$ learned by the NN correlates
with $m_{bc}$, we show a heatmap of the joint distribution of $(V,m_{bc})$ in
the upper-left panel of FIG. 2. Here, and in what follows, the heatmap is
generated using a separate test dataset with $10^{5}$ events. In the plot we
also show two nonparametric correlation coefficients, namely Kendall’s $\tau$
coefficient [22] and Spearman’s rank correlation coefficients $r_{s}$ [23]. A
value of $\pm 1$ for them would indicate one-to-one correspondence. Our
results depict an almost perfect correspondence between $V$ and $m_{bc}$.
Here, and in what follows, we append an overall minus sign to $V$ if needed,
to make the correlations positive and the plots in FIG. 2 intuitive.

Figure 2:  Top row: correlation plots between the artificial variable $V$ and
a relevant human-engineered variable for each of the three examples considered
in the text ($m_{bc}$, $m_{T}$ and $m_{T2}$ from left to right). Bottom row:
unit-normalized distributions of the corresponding artificial variable $V$ for
different mass inputs.

In practice, the artificial variable can be used to compare the data against
templates simulated for different values of $\Theta$. To illustrate this
usage, in the lower-left panel of FIG. 2 we show unit-normalized distributions
of the deep-learned variable $V$ for several different values of
$m_{A}=\\{200,280,320,400\\}$ GeV. It is seen that the distributions are
highly sensitive to the parameter choice $m_{A}$ and, if needed, $V$ can be
calibrated so that the peak location directly corresponds to $m_{A}$. The
observed spread around the peak values in the histogram, as well as the less
than perfect correspondence between $V$ and $m_{bc}$, are due to limitations
in the NN architecture and training.

Example 2: Semi-visible two body decay. Next we consider the semi-visible two-
body decay of a particle $A$ into a massless visible particle $b$ and a
possibly massive invisible particle $C$, $A\rightarrow b\,C$, where $A$ is
singly produced (with zero transverse momentum). The parameter $\Theta$ is
two-dimensional: $(m_{A},m_{C})$. The event $X$ is specified by the 4-momentum
$p_{b}=\\{E_{b},\vec{p}_{bT},p_{bz}\\}$ of $b$ and the missing transverse
momentum, leading to $d_{X}=6$.

We generate $(m_{A},m_{C})$ by uniformly sampling $(m_{A},\delta_{m})$ in the
region defined by $100\text{ GeV}\leq m_{A}\leq 500\text{ GeV}$ and
$0\leq\delta_{m}\leq m_{A}$, where $\delta_{m}\equiv(m_{A}-m_{C})^{2}/m_{A}$.
This choice of prior ensures that the relevant mass difference parameter in
this event topology $\mu\equiv(m_{A}^{2}-m_{C}^{2})/m_{A}$ is adequately
sampled in the range $[0,500]$ GeV. For the given value of $m_{A}$, we
generate an event as follows. The parent particle $A$ is boosted along the
beam axis $\pm z$ (with equal probability) to an energy chosen uniformly in
the range $[m_{A},1500\text{ GeV}]$. The particle $A$ is decayed on-shell into
$b$ and $C$, isotropically in its own rest frame. The details of network
architectures and training are the same as in Example 1.

The relevant variable for this event topology is the transverse mass $m_{T}$,
which in our setup is given by

$m_{T}(\tilde{m}_{C})\equiv p_{bT}+\sqrt{p_{bT}^{2}+\tilde{m}_{C}^{2}~{}},$
(5)

where the choice of mass ansatz $\tilde{m}_{C}$ for the mass of the invisible
particle $C$ does not affect the rank ordering of the events. For concreteness
in what follows we shall use $\tilde{m}_{C}=0$. The corresponding heatmap of
the joint distribution $(V,m_{T})$ and unit-normalized distributions of the
variable $V$ for several choices of $m_{A}=\\{200,280,320,400\\}$ GeV and
$m_{C}=100$ GeV are shown in the middle panels of FIG. 2. Once again, we
observe an almost perfect correlation between $V$ and $m_{T}$, and a high
sensitivity of the $V$ distributions to the input masses.

Example 3: Symmetric semi-visible two body decays. Finally, we consider the
exclusive production at a hadron collider of two equal-mass parent particles
$A_{1}$ and $A_{2}$ which decay semi-visibly as
$A_{1}A_{2}\rightarrow(b_{1}\,C_{1})\,(b_{2}\,C_{2})$. The parameter $\Theta$
is given by $(m_{A},m_{C})$, and the event $X$ is described by the 4-momemta
of $b_{1}$ and $b_{2}$, and the missing transverse momentum, leading to
$d_{X}=10$.

The masses $(m_{A},m_{C})$ are generated as in Example 2. In order to avoid
fine turning the network to the details of a particular collider, we uniformly
sampled the invariant mass $m_{A_{1}A_{2}}$ of the $A_{1}A_{2}$ system in the
range $[2m_{A},1500\text{ GeV}]$ and the lab-frame energy of the $A_{1}A_{2}$
system in the range $[m_{A_{1}A_{2}},2500\text{ GeV}]$. The direction of the
system was chosen to be along $\pm z$ with equal probability. The direction of
$A_{1}$ is chosen isotropically in the rest frame of the $A_{1}A_{2}$ system.
$A_{1}$ and $A_{2}$ are both decayed on-shell, isotropically in their
respective rest frames. The details of network architectures and training are
the same as in Example 1.

The straightforward generalization of the idea of the transverse mass to the
considered event topology leads to the stransverse mass variable
$m_{T2}(\tilde{m}_{C})$ [16]. In the upper-right panel of FIG. 2 we show a
heatmap of the joint distribution of $(V,m_{T2}(0))$, which reveals reasonably
good, but not perfect correlation, implying that the artificial event variable
encapsulates information beyond $m_{T2}$. This could have been expected for
the following two reasons: 1) unlike the previous two examples of singular
variables with sharp features in their distributions, $m_{T2}$ does not belong
to the class of singular variables [10]; 2) $m_{T2}$ only uses a subset of the
available kinematic information, namely the transverse momentum components. In
contrast, the artificial kinematic variable can use all of the available
information, and in a more optimal way. The lower-right panel of FIG. 2
displays unit-normalized distributions of the artificial variable for several
choices of $m_{A}$ and fixed $m_{C}=100$ GeV, again demonstrating the
sensitivity of $V$ to the mass spectrum.

Discussion and outlook. We proposed a new deep learning technique pictorially
summarized in FIG. 1 which allows the construction of event variables from a
set of training data produced from a given event topology. The novel component
is the simultaneous training for varying parameters $\Theta$ which allows the
algorithm to capture the underlying phase space structure irrespective of the
benchmark study point. This is the first such method for constructing event
variables with neural networks and can be applied to other, more challenging
event topologies in particle physics and beyond. In future applications of the
method one could enlarge the dimensionality of the latent space to $d_{V}>1$
and supplement the training data with additional features like tagging and
timing information, etc. By manipulating the specifics of the generation of
the training data, one can control what underlying physics effects are
available for the machine to learn from, and what physical parameters the
machine-learned variable will be sensitive to. Our method opens the door to
new investigations on intepretability and explainability by incorporating
modern representation learning approaches like contrastive learning [24].

Code and Data Availability. The code and data that support the findings of
this study are openly available at the following URL:
https://gitlab.com/prasanthcakewalk/code-and-data-availability/ under the
directory named arXiv_2105.xxxxx.

## Acknowledgements

We are indebted to the late Luc Pape for great insights and inspiration. This
work is supported in parts by US DOE DE-SC0021447 and DE-FG02-13ER41976. MP is
supported by Basic Science Research Program through the National Research
Foundation of Korea Research Grant No. NRF-2021R1A2C4002551. PS is partially
supported by the U.S. Department of Energy, Office of Science, Office of High
Energy Physics QuantISED program under the grants “HEP Machine Learning and
Optimization Go Quantum”, Award Number 0000240323, and “DOE QuantiSED
Consortium QCCFP-QMLQCF”, Award Number DE-SC0019219.

This manuscript has been authored by Fermi Research Alliance, LLC under
Contract No. DEAC02-07CH11359 with the U.S. Department of Energy, Office of
Science, Office of High Energy Physics.

## References

  * [1] R. Bellman, “Dynamic Programming”, Princeton University Press, 1957.
  * [2] T. Han, “Collider phenomenology: Basic knowledge and techniques,” [arXiv:hep-ph/0508097 [hep-ph]].
  * [3] A. J. Barr and C. G. Lester, “A Review of the Mass Measurement Techniques proposed for the Large Hadron Collider,” J. Phys. G 37, 123001 (2010) [arXiv:1004.2732 [hep-ph]].
  * [4] A. J. Barr, T. J. Khoo, P. Konar, K. Kong, C. G. Lester, K. T. Matchev and M. Park, “Guide to transverse projections and mass-constraining variables,” Phys. Rev. D 84, 095031 (2011) [arXiv:1105.2977 [hep-ph]].
  * [5] K. T. Matchev, F. Moortgat and L. Pape, “Dreaming Awake: Disentangling the Underlying Physics in Case of a SUSY-like Discovery at the LHC,” J. Phys. G 46, no.11, 115002 (2019) [arXiv:1902.11267 [hep-ph]].
  * [6] I. W. Kim, “Algebraic Singularity Method for Mass Measurement with Missing Energy,” Phys. Rev. Lett. 104, 081601 (2010) [arXiv:0910.1149 [hep-ph]].
  * [7] A. Rujula and A. Galindo, “Measuring the W-Boson mass at a hadron collider: a study of phase-space singularity methods,” JHEP 08, 023 (2011) [arXiv:1106.0396 [hep-ph]].
  * [8] A. De Rujula and A. Galindo, “Singular ways to search for the Higgs boson,” JHEP 06, 091 (2012) [arXiv:1202.2552 [hep-ph]].
  * [9] D. Kim, K. T. Matchev and P. Shyamsundar, “Kinematic Focus Point Method for Particle Mass Measurements in Missing Energy Events,” JHEP 10, 154 (2019) [arXiv:1906.02821 [hep-ph]].
  * [10] K. T. Matchev and P. Shyamsundar, “Singularity Variables for Missing Energy Event Kinematics,” JHEP 04, 027 (2020) [arXiv:1911.01913 [hep-ph]].
  * [11] M. Feickert and B. Nachman, “A Living Review of Machine Learning for Particle Physics,” [arXiv:2102.02770 [hep-ph]].
  * [12] S. Chang, T. Cohen and B. Ostdiek, “What is the Machine Learning?,” Phys. Rev. D 97, no.5, 056009 (2018) [arXiv:1709.10106 [hep-ph]].
  * [13] T. Faucett, J. Thaler and D. Whiteson, “Mapping Machine-Learned Physics into a Human-Readable Space,” Phys. Rev. D 103, no.3, 036020 (2021) [arXiv:2010.11998 [hep-ph]].
  * [14] V. D. Barger, A. D. Martin and R. J. N. Phillips, “Perpendicular $\nu_{e}$ Mass From $W$ Decay,” Z. Phys. C 21, 99 (1983).
  * [15] J. Smith, W. L. van Neerven and J. A. M. Vermaseren, “The Transverse Mass and Width of the $W$ Boson,” Phys. Rev. Lett. 50, 1738 (1983).
  * [16] C. G. Lester and D. J. Summers, “Measuring masses of semiinvisibly decaying particles pair produced at hadron colliders,” Phys. Lett. B 463, 99-103 (1999) [arXiv:hep-ph/9906349 [hep-ph]].
  * [17] N. Beaudry, “An intuitive proof of the data processing inequality,” Quantum Information & Computation, 12 (5–6): 432–441 (2012) [arXiv:1107.0740 [quant-ph]].
  * [18] N. Tishby, F. C. Pereira and W. Bialek, “The information bottleneck method”, Proc. of the 37-th Annual Allerton Conference on Communication, Control and Computing, page 368-377, (1999), [arXiv:physics/0004057 [physics.data-an]].
  * [19] R. Iten, T. Metger, H. Wilming, L. del Rio, R. Renner, “Discovering physical concepts with neural networks,” Phys. Rev. Lett. 124, 010508 (2020), [arXiv:1807.10300 [quant-ph]].
  * [20] D. Kim, K. Kong, K. T. Matchev, M. Park and P. Shyamsundar, in preparation.
  * [21] Martín Abadi, et. al,“Large-Scale Machine Learning on Heterogeneous Systems,” https://www.tensorflow.org/.
  * [22] Maurice G. Kendall, “The treatment of ties in ranking problems”, Biometrika Vol. 33, No. 3, pp. 239-251. 1945.
  * [23] Zwillinger, D. and Kokoska, S. “CRC Standard Probability and Statistics Tables and Formulae”. Chapman & Hall: New York. 2000. Section 14.7.
  * [24] P. Le-Khac, G. Healy and A. Smeaton, “Contrastive Representation Learning: A Framework and Review,” IEEE Access, vol. 8, pp. 193907-193934, 2020, [arXiv:2010.05113 [cs.LG]].


# Time-reversal symmetry breaking from lattice dislocations in superconductors

Clara N. Breiø  Niels Bohr Institute, University of Copenhagen, DK-2200
Copenhagen, Denmark  Andreas Kreisel  Niels Bohr Institute, University of
Copenhagen, DK-2200 Copenhagen, Denmark  Mercè Roig  Niels Bohr Institute,
University of Copenhagen, DK-2200 Copenhagen, Denmark  P. J. Hirschfeld
Department of Physics, University of Florida, Gainesville, Florida 32611, USA
Brian M. Andersen  Niels Bohr Institute, University of Copenhagen, DK-2200
Copenhagen, Denmark

###### Abstract

Spontaneous generation of time-reversal symmetry breaking in unconventional
superconductors is currently a topic of considerable interest. While chiral
superconducting order is often assumed to be the source of such signatures,
they can sometimes also arise from nonmagnetic disorder. Here we perform a
theoretical study of the impact of dislocations on the superconducting order
parameter within a microscopic one-band model which, in the homogeneous case,
features either extended $s$-wave, $d$-wave, or $s+id$-wave superconductivity
depending on the electron concentration. We find that the dislocations
minimize their impact on the superconducting condensate by inducing localized
supercurrents pinned by the dislocations, even well outside the $s+id$ regime.
We map out the parameter and density dependence of the induced currents. From
these results we conclude that quite generically unconventional
superconductors hosting dislocations tend to break time-reversal symmetry
locally.

##  I Introduction

A growing class of unconventional superconductors have been discovered to
exhibit evidence of spontaneous time-reversal symmetry breaking (TRSB) below
their critical temperatures $T_{c}$ [1, 2]. This is evidenced, for example, by
enhanced muon-spin relaxation ($\mu$SR) and polar Kerr effect measurements
reporting a change in the optical polar Kerr angle below $T_{c}$. This points
to internal magnetic fields spontaneously generated by the superconducting
state itself. A natural explanation for the occurrence of such fields would be
either in terms of a nonunitary triplet pairing state or multi-component
superconducting condensates entering a TRSB by a complex superposition of the
order. For superconducting instabilities condensing in two-dimensional (2D)
irreducible representations of the associated crystal point group, the latter
interpretation appears particularly appealing. Complex superpositions of two
symmetry-distinct order parameters may exhibit persistent supercurrents at
material edges, dislocations, or around various defect sites. The multi-
component scenario for the origin of TRSB below $T_{c}$ has been extensively
discussed for e.g. Sr2RuO4 [3, 4], and heavy-fermion compounds UPt3 [5, 6],
UTe2 [7, 8], PrOs4Sb12 [9, 10], and also for Ba1-xKxFe2As2 [11, 12], LaNiC2
[13], and more recently additionally for the kagome superconductors $A$V3Sb5
($A$: K, Rb, Cs) [14, 15, 16]. At present, however, the status of the precise
superconducting ground state remains controversial for several of these
materials (e.g. Sr2RuO4 and UTe2), partly due to specific heat data featuring
only a single transition that does not split under uniaxial strain.

Thus, it is important to pursue further possibilities for the origin of TRSB
generated inside the superconducting phase. In particular, one might ask what
mechanisms exist for TRSB for superconducting condensates composed of a single
component order parameter? The answer to this question necessarily involves
the presence of spatial inhomogeneities. For example, it is well-known that
nonmagnetic impurities can generate magnetic moments on nearby sites below
$T_{c}$ [17, 18, 19, 20, 21, 22, 21, 23, 24, 25, 26, 27, 28, 29, 30] in the
presence of electronic correlations. In addition, one can envision disorder-
induced slowing-down of magnetic fluctuations of preexisting magnetic
impurities, thereby lowering the characteristic fluctuation frequencies into
the muon time window. Another possibility for TRSB includes the generation of
localized orbital loops of supercurrents. The latter is well-known to arise
near nonmagnetic disorder sites in complex TRSB multi-component condensates
[31, 32, 33, 34, 35, 36, 37, 38], but was recently shown to be also present in
the strong disorder limit of single-component superconductors [39, 40].
Likewise, one might expect dislocations and grain boundaries to similarly
operate as seeds of localized loop currents. Recently, the latter scenario was
proposed by Willa et al.,[41] as a possible explanation for reconciling the
existence of TRSB and a single specific heat transition in Sr2RuO4. Indeed,
for this material it is well-known that dislocations are prevalent in many
samples [42].

Here, motivated by the above-mentioned developments we perform a detailed
theoretical study of the possibility of generation of localized supercurrents
by dislocations in single-component order parameter superconductors. For
concreteness, we focus on a one-band model with nearest neighbor attractive
interactions. This model supports either $B_{1g}$ ($d$-wave) or $A_{1g}$
(extended $s$-wave) spin-singlet superconductivity in the homogeneous case,
depending on filling. In addition, there exists a coexistence regime where the
complex combination $s+id$ is favored. We start from a single-band picture
without any material-specific details, except the assumption that the orbital
states are of $d_{x^{2}-y^{2}}$ symmetry. To obtain the modifications of the
hopping integrals due to a dislocation, we remove a line of $l$ atoms and
simulate the equilibrium positions assuming parabolic potentials between the
atoms in a molecular dynamics setup. The hopping integrals are then the
expectation value of the kinetic energy in atomic-like $d_{x^{2}-y^{2}}$
orbitals separated by the distances determined in this manner. From
selfconsistent real-space studies and associated calculations of the free
energy, we find that the dislocations favor the local generation of complex
pair potential order in a rather broad range of parameters and electron
concentrations. This implies that TRSB dislocation-bound orbital currents are
spontaneously stabilized. Our results highlight the role of spatial
inhomogeneity in generating TRSB and associated magnetic signals within the
superconducting phase.

##  II Model and Method

###  II.1 Dislocations and lattice relaxation

The normal state Hamiltonian of a system of real space points with hopping
matrix elements $t_{ij}$ is given by

$H_{0}=\sum_{i,j,\sigma}(t_{ij}-\mu\delta_{ij})c^{\dagger}_{i\sigma}c_{j\sigma},$
(1)

where the operator $c^{\dagger}_{i\sigma}$ creates an electron at real-space
point $i$ with spin $\sigma$. For a homogeneous system, we set $t_{ij}=t=-1$
if $i$ and $j$ are nearest neighbors (NN). In the presence of a dislocation,
the displacements of lattice points modifies the hopping integrals which we
capture in the following way. Given the distance vector
$\mathbf{r}_{ij}=(x_{ij},y_{ij})$, we calculate the hopping from the
expectation value $t(\mathbf{r})$ of the kinetic energy
$-1/(2m)\mathbf{\nabla}^{2}$ for overlapping atomic $d_{x^{2}-y^{2}}$
wavefunctions. A map of the hopping integral is shown in Fig. 1(a) where we
used the value $|t([1,0])|=1$ to fix the effective mass $m$. Black circle
shows $r_{\mathrm{cut}}=(1+\sqrt{2})/2$ defining the maximum distance allowed
for NN bonds. Hoppings for $|\mathbf{r}_{ij}|>r_{\mathrm{cut}}$ are simply set
to zero. Note that the implementation of a dislocation reduces the total
number of sites in an $N\times N$ square lattice to $N^{2}-l$.

Figure 1: (a) Hopping integral as function of connection vector $\mathbf{r}$
between sites. Black circle shows $r_{\mathrm{cut}}=(1+\sqrt{2})/2$ defining
the maximum distance allowed for NN bonds. Grey dots indicate site positions
of the homogeneous lattice. (b) An example of real space points (black dots)
describing a $l=5$ dislocation together with the value of the hopping integral
(colored bonds). The inset shows a zoom-in of the region close to the
dislocation.

The positions of the lattice points are obtained using a molecular dynamics
simulation with periodic boundary conditions where the sites of a square
lattice are bound to NN and NNN sites via unstretched strings of equal
strength. Next, a finite number $l$ of lattice sites along a line are removed
to form a dislocation pair and bonds to the “new” NN and NNN along that line
are established. The relaxation under the harmonic potential with periodic
boundary conditions is then done to find all positions $\mathbf{r}_{i}$ of the
lattice points. This procedure yields symmetric positions close to the
dislocation and reproduces the strain expected from continuum theory far away
from the dislocation. An example of the lattice positions (black dots) for
$l=5$ is shown in Fig. 1(b) where the value of the hopping integrals are
indicated by the color of the corresponding bond.

###  II.2 Selfconsistent solutions of superconductivity

The study presented here is based on the following one-band BCS mean-field
Hamiltonian

$\displaystyle H$ $\displaystyle=H_{0}-\sum_{\langle
i,j\rangle}(\Delta_{ij}(c_{i\downarrow}^{\dagger}c_{j\uparrow}^{\dagger}-c_{i\uparrow}^{\dagger}c_{j\downarrow}^{\dagger})+\mbox{H.c.}),$
(2)

where $\Delta_{ij}=V(\langle c_{j\uparrow}c_{i\downarrow}\rangle-\langle
c_{j\downarrow}c_{i\uparrow}\rangle)$ is the attractive ($V>0$) singlet
pairing on NN bonds also defined by $|\mathbf{r}_{ij}|<r_{cut}$. The
corresponding Bogoliubov-de Gennes (BdG) equations are solved selfconsistently
for systems of size $N\times N$ with a dislocation of length $l$ as described
above. The convergence criterion is set to
$|\mathcal{H}^{s}_{\alpha\beta}-\mathcal{H}^{s-1}_{\alpha\beta}|<10^{-8}\,|t|\,\,\forall\,\,\alpha,\beta$,
where $\mathcal{H}$ is the BdG Hamiltonian and $s$ denotes the iteration
counter. All results are computed with open boundary conditions (OBC) to avoid
super-cell interference effects and we set $V=0.75\,|t|$ and $N\times
N=39\times 39$ unless otherwise specified. Additionally, the average electron
density will be fixed by adjusting the chemical potential, $\mu$, in the
selfconsistent iterations. Results presented here are obtained for the hole-
doped region of the phase diagram ($\langle n\rangle<1$), but applies also to
the electron-doped region due to particle-hole symmetry at half-filling
($\mu=0$) with the symmetry $\langle n\rangle\rightarrow 2-\langle n\rangle$.

To study the possibility of spontaneous TRSB induced by the atomic
dislocations and manifested by generation of orbital currents, the bond-
current densities of all converged results are computed as

$\displaystyle\langle j_{ij}\rangle$ $\displaystyle=-i\frac{e}{\hbar
a^{2}}\sum_{\sigma}t_{ij}\langle
c^{\dagger}_{i\sigma}c_{j\sigma}-c^{\dagger}_{j\sigma}c_{i\sigma}\rangle$ (3)
$\displaystyle=2\frac{e}{\hbar
a^{2}}\sum_{n}t_{ij}\mathrm{Im}[u^{*}_{in\downarrow}u_{jn\downarrow}f(E_{n})-v^{*}_{in\uparrow}v_{jn\uparrow}f(E_{n})],$

where $e$ denotes the electron charge, $u_{in\downarrow},\,v_{in\uparrow}$ are
the usual eigenvector components from the Bogoliubov transformation and
$f(E_{n})$ is the Fermi-Dirac distribution function. All results are obtained
with the temperature fixed at $T=0.001\,|t|$. The selfconsistency ensures
current conservation as confirmed by incoming and outgoing currents being
equal in magnitude at all lattice sites [43]. We project the bond-current
densities to local densities using

$\displaystyle\mathbf{j}_{i}=\frac{1}{2}\sum_{\langle i,j\rangle}\frac{\langle
j_{ij}\rangle}{|\mathbf{r}_{ij}|}\begin{pmatrix}x_{ij}\\\
y_{ij}\end{pmatrix}.$ (4)

###  II.3 Free energy calculation

Introducing spatial inhomogeneities in the selfconsistent solution of the
model vastly extends the parameter space in which the free energy should be
minimized. As the selfconsistency only ensures results to be at stationary
points in the free energy landscape, an enhancement of parameter space also
increases the risk of obtaining converged results which do not represent the
global minimum. Thus, to ensure that the spontaneous TRSB and generation of
local loop currents are indeed energetically favorable, we compare the free
energy of converged configurations with and without TRSB. The latter is
obtained by restricting the superconducting order parameters to take real
values at all bonds throughout the system. In this context, it should be noted
that while Eq. (2) is written in the grand canonical ensemble we compare
solutions with fixed particle number. Thus, the relevant quantity to consider
is the free energy rather than the grand potential i.e.,
$F=\Omega+\mu\hat{N}$, where $\Omega$ is the grand potential and
$\hat{N}=\sum_{i\sigma}\langle c^{\dagger}_{i\sigma}c_{i\sigma}\rangle$. In
the zero temperature limit, this can be obtained as $F=\langle
H\rangle+\mu\hat{N}$ yielding simply

$\displaystyle F$ $\displaystyle=\sum_{\langle
i,j\rangle,\,\sigma}t_{ij}\langle
c^{\dagger}_{i\sigma}c_{j\sigma}\rangle-\sum_{\langle
i,j\rangle}\frac{|\Delta_{ij}|^{2}}{V}\,,$ (5)

where the expectation value is calculated using the eigenvalues $E_{n}$ and
eigenvector components $u_{in\downarrow},\,v_{in\uparrow}$ of the Hamiltonian
Eq. (2). Note that the constant term originating from the mean-field
decoupling of the interactions which is not included in Eq. (2) cancels the
contribution from the expectation value of one of the superconducting pairing
terms.

##  III Results

###  III.1 Homogeneous phase diagram

To establish a baseline for the results obtained for systems with dislocation
defects, the phase diagram for the homogeneous system as a function of
electron density is shown in Fig. 2. The magnitude of the $A_{1g}$ (extended
$s$-wave, red line) and $B_{1g}$ ($d_{x^{2}-y^{2}}$, blue line) order
parameters are shown on the left axis while the site-averaged free energy
difference, $\overline{\delta
F}=(F_{\Delta\in\mathbb{C}}-F_{\Delta\in\mathbb{R}})/N^{2}$, is indicated on
the right axis (black line). Here $F$ is computed according to Eq. (5) and the
subscript indicates whether $\Delta_{ij}$ are free to take on any value
($\Delta\in\mathbb{C}$) or they are restricted to real values preventing
currents from developing ($\Delta\in\mathbb{R}$). For electron densities close
to half-filling, $\langle n\rangle=1$, the $A_{1g}$ solution is strongly
disfavored due to the near-alignment of the Fermi surface and line nodes, and
therefore the selfconsistent solution exhibits $B_{1g}$ symmetry. On the
contrary, for small Fermi pockets ($\langle n\rangle<0.36$), a full gap can
develop in the $A_{1g}$ channel which is more favorable than the $B_{1g}$
channel hosting four point nodes at the Fermi surface for all electron
densities [44, 45, 38].

Figure 2: Phase diagram showing on the left axis the amplitude of the
homogeneous $A_{1g}$ (extended $s$-wave, red line) and $B_{1g}$
($d_{x^{2}-y^{2}}$, blue line) order parameters as a function of electron
density. The free energy difference per site, $\overline{\delta
F}=(F_{\Delta\in\mathbb{C}}-F_{\Delta\in\mathbb{R}})/N^{2}$, is displayed on
the right axis (black line). Results are computed on a $201\times 201$ square
lattice with periodic boundary conditions.

In the intermediate region ($0.36<\langle n\rangle<0.48$) a complex
superposition of the two channels is favored, yielding a TRSB solution of the
form $A_{1g}\pm iB_{1g}$. A possible nematic solution $A_{1g}\pm B_{1g}$
exhibits line nodes at $k_{x}=\pm\pi/2$, while the TRSB solution reduces the
line nodes to point nodes at $k_{x}=\pm k_{y}=\pm\pi/2$, resulting in a gain
in condensation energy. This energy gain is reflected by the free energy
difference shown in black in Fig. 2. Note by the definition stated above, that
a negative value of $\overline{\delta F}$ signifies a TRSB ground state. Pure
$A_{1g}/B_{1g}$ regions have $\overline{\delta F}=0$, since the two solutions
($\Delta\in\mathbb{C}$ and $\Delta\in\mathbb{R}$) can only differ by a change
in the global superconducting phase.

###  III.2 Loop currents and pairing modulations

Several earlier studies investigated the emergence of orbital currents induced
by non-magnetic defects in TRSB superconductors employing both
phenomenological Ginzburg-Landau approaches and microscopic mean-field methods
[33, 46, 31, 47, 34, 39, 40]. In particular, the studies reveal that
superconductors with $A_{1g}+iB_{1g}$ pairing symmetries generate spontaneous
supercurrents in the vicinity of various types of non-magnetic potential
scatterers, including point-like and spatially extended impurities [40, 33],
system boundaries [31] as well as defects exhibiting non-trivial spatial
structures [34]. While the dislocations introduced in this work do not act as
simple potential scatterers, it is reasonable to assume that the spatial
inhomogeneities also generate orbital localized currents in the
$A_{1g}+iB_{1g}$ phase due to local fluctuations in the superconducting order
parameter analogous to the results presented in previous studies [33, 46, 31,
47, 34, 39, 40]. Indeed, the computed selfconsistent results for systems in
the homogeneous TRSB density range including a dislocation in the center
region confirm the formation of local loop currents surrounding the defect. We
note that, while the specific current pattern and magnitude depend on the
dislocation length, the formation of currents is robust for all dislocations
studied here.

Figure 3: (a) Site-resolved magnitude of the current density. Localized loop
currents have developed near the dislocation with a small, long-range, tail
required by current conservation. Zoom-in show the local current direction
indicated by arrows. Arrow lengths have been re-scaled for visual clarity.
(b,c) Modulations of $\Delta_{ij}$ for the converged solution with (b) and
without (c) currents. Bond thickness indicate $|\Delta_{ij}|$ and bond color
display modulations of the bond phase, $\phi_{ij}$. Here $\lfloor.\rfloor$ is
the floor-function such that the nearest multiple of $\pi$ is subtracted and
only the deviations in $\phi_{ij}$ away from pure $B_{1g}$ structure is shown.
All results are computed with an $l=5$ dislocation at $\langle n\rangle=0.64$.

While the generation of defect-induced supercurrents in the $A_{1g}+iB_{1g}$
phase might be expected, the emergence of currents in the pure $A_{1g}$ or
$B_{1g}$ phases would require spontaneous TRSB occurring solely due to the
lattice inhomogeneity associated with the presence of a dislocation.
Surprisingly, we find that this scenario is indeed realized for a range of
densities well within the homogeneous $B_{1g}$ phase. This is unlike the
property of single impurities which are unable to seed orbital currents in
this region. Figure 3(a) displays the induced orbital currents for result
obtained at $\langle n\rangle=0.64$ and $l=5$. The supercurrents are
predominantly localized in the region extending $\sim 5$ lattice sites away
from the dislocation, but current conservation requires a small, yet finite,
long-range tail, see left panel of Fig. 3(a). The lack of corner/edge currents
emphasises that the near-homogeneous regions remain in the $B_{1g}$ phase and,
as such, the apparent breaking of TRS is a purely local and dislocation-
induced effect. The zoom-in shown in the right panel displays the current
direction, demonstrating the orbital loops in the supercurrent pattern. The
magnitudes are indicated by arrow color, while the arrow lengths have been re-
scaled for visual clarity.

Since supercurrents arise from modulations in the phase and amplitude of the
order parameter, we display the bond order parameters corresponding to Fig.
3(a) in Fig. 3(b). For comparison, Fig. 3(c) shows the converged result for
identical parameters with the $\Delta_{ij}\in\mathbb{R}$ restriction applied
(i.e. no currents allowed). The line thickness indicate the amplitudes of the
bond pairings while the color displays the phase-deviation away from
$\phi_{ij}=0,\pm\pi$, i.e. the deviation of the superconducting phase away
from a pure $B_{1g}$ phase with e.g. $\phi_{ij}=0$ on $x$-bonds and
$\phi_{ij}=\pi$ on $y$-bonds. The restriction applied to the result shown in
Fig. 3(c) obviously hinders phase modulations and the spatial inhomogeneity
associated with the dislocation mainly serve as a pair-breaking effect as
evident from the smaller gap amplitude on the bonds near the ends. In
contrast, the current-carrying solution shown in Fig. 3(b) features clear
deviations from the pure $B_{1g}$ phase up to $18\%$ alongside the amplitude
fluctuations. Explicit comparison of the free energy difference per site in
this particular case yields indeed a negative $\overline{\delta F}=-2.4\cdot
10^{-6}|t|$. This is much smaller than the typical reduction of the free
energy when a homogeneous $A_{1g}+iB_{1g}$ state emerges; again a signature of
the local nature of TRSB in this scenario.

###  III.3 Dislocation phase diagram and re-emerging supercurrents

Figure 4: Left axis: Site-averaged local current density magnitude (dashed
lines, circle markers). A main TRSB dome extends from the homogeneous TRSB
region (marked in grey) with re-emerging domes extending to $\langle
n\rangle\sim 0.7$. Right axis: Site-averaged free energy gain (solid lines,
triangle markers) associated with the spontaneous TRSB. All TRSB solutions are
energetically favored or, in rare cases, degenerate with the time-reversal
symmetric solution. Colors indicate different dislocation lengths, $l$. Note
that the total number of sites used in the averages is given by $N^{2}-l$.

Given the rather unexpected finding of energetically favorable loop currents
far inside the homogeneous $B_{1g}$ phase, it is desirable to perform a more
extensive study of this phenomenon as a function of dislocation length and
electron density. Thus, in Fig. 4 the left axis shows site-averaged
supercurrent (dashed lines, circle markers) for dislocation lengths ranging
from $l=1$ to $l=10$ as a function of density. As a direct reference to the
phase diagram in Fig. 2, the homogeneous TRSB region has been marked in grey.
It is evident that all dislocation lengths induce currents within this region
as anticipated from the previous studies discussed above. We note, however,
that the dome emerging from the homogeneous TRSB region is significantly
enhanced with the upper phase boundary shifted from $\langle n\rangle=0.48$ to
$\langle n\rangle\sim 0.53$.

More interestingly, the generation of loop currents is clearly extended far
beyond this main dome with re-emerging TRSB domes ranging all the way to
$\langle n\rangle\sim 0.7$. While the density ranges and current strengths of
the TRSB domes varies for different dislocation lengths, the multi-dome
feature is universal across all systems investigated here, where some show an
additional shoulder at the edge of the main dome (e.g. $l=9$, marked in
orange) and other show clearly distinguishable domes (e.g. $l=2$, marked in
purple, or $l=5$, marked in cyan).

Figure 4 also shows the free energy gain (solid lines, triangle markers)
associated with the generation of loop currents on the right axis verifying
that the TRSB solution is energetically preferred with the exception of a few
points with degenerate time-reversal symmetric and TRSB solutions down to the
convergence criterion, see e.g. $l=7$ (light green) at $\langle
n\rangle=0.68$. Note that the individual contributions from the two terms in
Eq. (5) varies between loss and gain for different parameters and only the
total free energy difference reveal a consistent gain across all TRSB
solutions.

##  IV Discussion and conclusions

It is relevant to highlight a fundamental difference between the spatial
inhomogeneity introduced by dislocation defects and multiple impurities acting
as dopants. Refs. [39] and [40] both found qualitatively similar results for
strongly disordered systems ($>20\%$ dopant atoms) where spontaneous loop
currents are formed well within the homogeneous $B_{1g}$ phase. However, such
strong disorder breaks all spatial symmetries and, as argued in Ref. [40],
unavoidably introduce finite components of the other pairing symmetry
channels. In contrast, the dislocation defects studied in this work preserves
a global $C_{2}$-symmetry ($C_{4}$-symmetry for $l=1$) around
$\mathbf{r}=(0,0)$ and, from a naive perspective, it is reasonable to assume
that disorder with a higher symmetry should have less impact on the preferred
symmetry of $\Delta_{ij}$. However, as evident from the results presented
here, the disorder does not need to be strong from a symmetry-breaking
perspective to obtain TRSB ground states.

This is quite important, because many unconventional superconducting crystals
are rather clean by many standard criteria, yet undoubtedly contain
significant concentrations of such linear defects oriented by strain, which
may be preferentially along high symmetry axes. For example, the highest
quality Sr2RuO4 crystals are famously extremely clean; tiny amounts of
nonmagnetic disorder (e.g. 0.15% Ti4+ substituting for Ru4+) suppress $T_{c}$
to zero [48, 49, 50]. Both Kerr and $\mu$SR measurements indicate a small TRSB
signal at the superconducting transition, yet the “standard” explanations in
terms of two-component order parameters are seriously challenged by the
absence of any splitting of the specific heat transition with or without
strain [51, 52], suggesting alternative explanations for the weak TRSB. One
such explanation was provided by Willa et al.[41], who assumed a close
proximity of two pair components, $d$\- and $g$-wave such that a TRSB state
was nucleated near an edge dislocation. Here we have demonstrated a similar
scenario within a microscopic model, and also shown that the degeneracy need
not be particularly close.

Additionally, apart from a few bonds located at the dislocation edges, all
bonds exhibit near-homogeneous values of the hopping integrals leading to the
expectation that a near-homogeneous $B_{1g}$ phase should develop. To examine
this in detail, we performed a systematic investigation of the amount of
lattice distortion necessary to reproduce TRSB domes within the homogeneous
$B_{1g}$ phase and found within additional self-consistent calculations that
even small modifications in the hopping integrals down to $5\%$ must be
included to qualitatively reproduce the results shown in Fig. 4. To
quantitatively reproduce the results, all relaxation effects must be included.
Furthermore, we have also found that the particular density ranges and current
magnitudes shown in Fig. 4 varies somewhat with system size and boundary
conditions (e.g. employing periodic boundaries). We stress, however, that the
generic re-emergence of TRSB ground states within the homogeneous $B_{1g}$
phase persists for a vast range of investigated system sizes and boundary
conditions. These combined findings suggest that the weak long-range
perturbations and currents induced by the dislocation are essential. In this
context, it is worth mentioning that typical dislocations in crystals may
involve many hundreds of atoms, which presumably enhances these long-range
effects.

In summary, we have performed a theoretical study of the preferred
superconducting order near dislocations of a microscopic model that
incorporates two symmetry-distinct superconducting order parameters. It is
found that dislocations tend to nucleate orbital currents and thereby locally
break time-reversal symmetry spontaneously. This takes place even in regions
of the phase diagram where the corresponding homogeneous superconducting
condensate is single component order. This effect may be important for the
understanding of experiments reporting TRSB at $T_{c}$ in unconventional
superconductors.

## Acknowledgements

We acknowledge useful discussions with C. Hicks and M. Pal. C.N.B. and A.K.
acknowledge support by the Danish National Committee for Research
Infrastructure (NUFI) through the ESS-Lighthouse Q-MAT. P.J.H. acknowledges
support from NSF-DMR-1849751. M.R. acknowledges support from the Novo Nordisk
Foundation grant NNF20OC0060019.

## References

  * Ghosh  _et al._ [2021] S. K. Ghosh, M. Smidman, T. Shang, J. F. Annett, A. D. Hillier, J. Quintanilla, and H. Yuan, Recent progress on superconductors with time-reversal symmetry breaking, J. Phys. Condens. Matter 33, 033001 (2021).
  * Røising  _et al._ [2022] H. S. Røising, G. Wagner, M. Roig, A. T. Rømer, and B. M. Andersen, Heat capacity double transitions in time-reversal symmetry broken superconductors, Phys. Rev. B 106, 174518 (2022).
  * Maeno  _et al._ [1994] Y. Maeno, H. Hashimoto, K. Yoshida, S. Nishizaki, T. Fujita, J. G. Bednorz, and F. Lichtenberg, Superconductivity in a layered perovskite without copper, Nature (London) 372, 532 (1994).
  * Mackenzie  _et al._ [2017] A. P. Mackenzie, T. Scaffidi, C. W. Hicks, and Y. Maeno, Even odder after twenty-three years: the superconducting order parameter puzzle of Sr2RuO4, npj Quantum Mater. 2, 40 (2017).
  * Stewart [1984] G. R. Stewart, Heavy-fermion systems, Rev. Mod. Phys. 56, 755 (1984).
  * Joynt and Taillefer [2002] R. Joynt and L. Taillefer, The superconducting phases of ${\mathrm{UPt}}_{3}$, Rev. Mod. Phys. 74, 235 (2002).
  * Ran  _et al._ [2019] S. Ran, C. Eckberg, Q.-P. Ding, Y. Furukawa, T. Metz, S. R. Saha, I.-L. Liu, M. Zic, H. Kim, J. Paglione, and N. P. Butch, Nearly ferromagnetic spin-triplet superconductivity, Science 365, 684 (2019).
  * Aoki  _et al._ [2022] D. Aoki, J.-P. Brison, J. Flouquet, K. Ishida, G. Knebel, Y. Tokunaga, and Y. Yanase, Unconventional superconductivity in UTe2, J. Phys.: Condens. Matter 34, 243002 (2022).
  * Bauer  _et al._ [2002] E. D. Bauer, N. A. Frederick, P.-C. Ho, V. S. Zapf, and M. B. Maple, Superconductivity and heavy fermion behavior in ${\mathrm{PrOs}}_{4}{\mathrm{Sb}}_{12}$, Phys. Rev. B 65, 100506(R) (2002).
  * Aoki  _et al._ [2003] Y. Aoki, A. Tsuchiya, T. Kanayama, S. R. Saha, H. Sugawara, H. Sato, W. Higemoto, A. Koda, K. Ohishi, K. Nishiyama, and R. Kadono, Time-Reversal Symmetry-Breaking Superconductivity in Heavy-Fermion ${\mathrm{P}\mathrm{r}\mathrm{O}\mathrm{s}}_{4}{\mathrm{S}\mathrm{b}}_{12}$ Detected by Muon-Spin Relaxation, Phys. Rev. Lett. 91, 067003 (2003).
  * Rotter  _et al._ [2008] M. Rotter, M. Tegel, and D. Johrendt, Superconductivity at 38 K in the Iron Arsenide $({\mathrm{Ba}}_{1-x}{\mathrm{K}}_{x}){\mathrm{Fe}}_{2}{\mathrm{As}}_{2}$, Phys. Rev. Lett. 101, 107006 (2008).
  * Böker  _et al._ [2017] J. Böker, P. A. Volkov, K. B. Efetov, and I. Eremin, $s+is$ superconductivity with incipient bands: Doping dependence and STM signatures, Phys. Rev. B 96, 014517 (2017).
  * Hillier  _et al._ [2009] A. D. Hillier, J. Quintanilla, and R. Cywinski, Evidence for Time-Reversal Symmetry Breaking in the Noncentrosymmetric Superconductor ${\mathrm{LaNiC}}_{2}$, Phys. Rev. Lett. 102, 117007 (2009).
  * Guguchia  _et al._ [2023] Z. Guguchia, C. Mielke, D. Das, R. Gupta, J.-X. Yin, H. Liu, Q. Yin, M. H. Christensen, Z. Tu, C. Gong, N. Shumiya, M. S. Hossain, T. Gamsakhurdashvili, M. Elender, P. Dai, A. Amato, Y. Shi, H. C. Lei, R. M. Fernandes, M. Z. Hasan, H. Luetkens, and R. Khasanov, Tunable unconventional kagome superconductivity in charge ordered RbV3Sb5 and KV3Sb5, Nature Communications 14, 153 (2023).
  * Mielke  _et al._ [2022] C. Mielke, D. Das, J.-X. Yin, H. Liu, R. Gupta, Y.-X. Jiang, M. Medarde, X. Wu, H. C. Lei, J. Chang, P. Dai, Q. Si, H. Miao, R. Thomale, T. Neupert, Y. Shi, R. Khasanov, M. Z. Hasan, H. Luetkens, and Z. Guguchia, Time-reversal symmetry-breaking charge order in a kagome superconductor, Nature 602, 245 (2022).
  * Rømer  _et al._ [2022] A. T. Rømer, S. Bhattacharyya, R. Valentí, M. H. Christensen, and B. M. Andersen, Superconductivity from repulsive interactions on the kagome lattice, Phys. Rev. B 106, 174514 (2022).
  * Tsuchiura  _et al._ [2001] H. Tsuchiura, Y. Tanaka, M. Ogata, and S. Kashiwaya, Local magnetic moments around a nonmagnetic impurity in the two-dimensional t-J model, Phys. Rev. B 64, 140501 (2001).
  * Wang and Lee [2002] Z. Wang and P. A. Lee, Local moment formation in the superconducting state of a doped mott insulator, Phys. Rev. Lett. 89, 217002 (2002).
  * Zhu  _et al._ [2002] J.-X. Zhu, I. Martin, and A. R. Bishop, Spin and Charge Order around Vortices and Impurities in High-${T}_{c}$ Superconductors, Phys. Rev. Lett. 89, 067003 (2002).
  * Chen and Ting [2004] Y. Chen and C. S. Ting, States of Local Moment Induced by Nonmagnetic Impurities in Cuprate Superconductors, Phys. Rev. Lett. 92, 077203 (2004).
  * Andersen  _et al._ [2007] B. M. Andersen, P. J. Hirschfeld, A. P. Kampf, and M. Schmid, Disorder-Induced Static Antiferromagnetism in Cuprate Superconductors, Phys. Rev. Lett. 99, 147002 (2007).
  * Harter  _et al._ [2007] J. W. Harter, B. M. Andersen, J. Bobroff, M. Gabay, and P. J. Hirschfeld, Antiferromagnetic correlations and impurity broadening of NMR linewidths in cuprate superconductors, Phys. Rev. B 75, 054520 (2007).
  * Andersen  _et al._ [2010] B. M. Andersen, S. Graser, and P. J. Hirschfeld, Disorder-Induced Freezing of Dynamical Spin Fluctuations in Underdoped Cuprate Superconductors, Phys. Rev. Lett. 105, 147002 (2010).
  * Schmid  _et al._ [2010] M. Schmid, B. M. Andersen, A. P. Kampf, and P. J. Hirschfeld, d-Wave superconductivity as a catalyst for antiferromagnetism in underdoped cuprates, New Journal of Physics 12, 053043 (2010).
  * Gastiasoro  _et al._ [2013] M. N. Gastiasoro, P. J. Hirschfeld, and B. M. Andersen, Impurity states and cooperative magnetic order in Fe-based superconductors, Phys. Rev. B 88, 220509 (2013).
  * Gastiasoro  _et al._ [2014] M. N. Gastiasoro, I. Paul, Y. Wang, P. J. Hirschfeld, and B. M. Andersen, Emergent Defect States as a Source of Resistivity Anisotropy in the Nematic Phase of Iron Pnictides, Phys. Rev. Lett. 113, 127001 (2014).
  * Gastiasoro and Andersen [2015] M. N. Gastiasoro and B. M. Andersen, Local Magnetization Nucleated by Non-magnetic Impurities in Fe-based Superconductors, Journal of Superconductivity and Novel Magnetism 28, 1321 (2015).
  * Martiny  _et al._ [2015] J. H. J. Martiny, M. N. Gastiasoro, I. Vekhter, and B. M. Andersen, Impurity-induced antiferromagnetic order in Pauli-limited nodal superconductors: Application to heavy-fermion ${\mathrm{CeCoIn}}_{5}$, Phys. Rev. B 92, 224510 (2015).
  * Gastiasoro  _et al._ [2016] M. N. Gastiasoro, F. Bernardini, and B. M. Andersen, Unconventional Disorder Effects in Correlated Superconductors, Phys. Rev. Lett. 117, 257002 (2016).
  * Martiny  _et al._ [2019] J. H. J. Martiny, A. Kreisel, and B. M. Andersen, Theoretical study of impurity-induced magnetism in FeSe, Phys. Rev. B 99, 014509 (2019).
  * Lee  _et al._ [2009] W.-C. Lee, S.-C. Zhang, and C. Wu, Pairing state with a time-reversal symmetry breaking in FeAs-based superconductors, Phys. Rev. Lett. 102, 217002 (2009).
  * Garaud and Babaev [2014] J. Garaud and E. Babaev, Domain Walls and Their Experimental Signatures in $s+is$ Superconductors, Phys. Rev. Lett. 112, 017003 (2014).
  * Maiti  _et al._ [2015] S. Maiti, M. Sigrist, and A. Chubukov, Spontaneous currents in a superconductor with $s+is$ symmetry, Phys. Rev. B 91, 161102 (2015).
  * Lin  _et al._ [2016] S.-Z. Lin, S. Maiti, and A. Chubukov, Distinguishing between $s+id$ and $s+is$ pairing symmetries in multiband superconductors through spontaneous magnetization pattern induced by a defect, Phys. Rev. B 94, 064519 (2016).
  * Silaev  _et al._ [2017] M. Silaev, J. Garaud, and E. Babaev, Phase diagram of dirty two-band superconductors and observability of impurity-induced $s+is$ state, Phys. Rev. B 95, 024517 (2017).
  * Garaud  _et al._ [2016] J. Garaud, M. Silaev, and E. Babaev, Thermoelectric Signatures of Time-Reversal Symmetry Breaking States in Multiband Superconductors, Phys. Rev. Lett. 116, 097002 (2016).
  * Benfenati  _et al._ [2020] A. Benfenati, M. Barkman, T. Winyard, A. Wormald, M. Speight, and E. Babaev, Magnetic signatures of domain walls in $s+is$ and $s+id$ superconductors: Observability and what that can tell us about the superconducting order parameter, Phys. Rev. B 101, 054507 (2020).
  * Roig  _et al._ [2022] M. Roig, A. T. Rømer, P. J. Hirschfeld, and B. M. Andersen, Revisiting superconductivity in the extended one-band hubbard model: Pairing via spin and charge fluctuations, Phys. Rev. B 106, 214530 (2022).
  * Li  _et al._ [2021a] Z.-X. Li, S. A. Kivelson, and D.-H. Lee, Superconductor-to-metal transition in overdoped cuprates, npj Quantum Materials 6, 36 (2021a).
  * Breiø  _et al._ [2022] C. N. Breiø, P. J. Hirschfeld, and B. M. Andersen, Supercurrents and spontaneous time-reversal symmetry breaking by nonmagnetic disorder in unconventional superconductors, Phys. Rev. B 105, 014504 (2022).
  * Willa  _et al._ [2021] R. Willa, M. Hecker, R. M. Fernandes, and J. Schmalian, Inhomogeneous time-reversal symmetry breaking in ${\mathrm{Sr}}_{2}{\mathrm{RuO}}_{4}$, Phys. Rev. B 104, 024511 (2021).
  * Ying  _et al._ [2013] Y. A. Ying, N. E. Staley, Y. Xin, K. Sun, X. Cai, D. Fobes, T. J. Liu, Z. Q. Mao, and Y. Liu, Enhanced spin-triplet superconductivity near dislocations in Sr2RuO4, Nature Communications 4, 2596 (2013).
  * Andersen  _et al._ [2006] B. M. Andersen, I. V. Bobkova, P. J. Hirschfeld, and Y. S. Barash, $0-\pi$ Transitions in Josephson Junctions with Antiferromagnetic Interlayers, Phys. Rev. Lett. 96, 117005 (2006).
  * Kreisel  _et al._ [2017] A. Kreisel, A. T. Rømer, P. J. Hirschfeld, and B. M. Andersen, Superconducting Phase Diagram of the Paramagnetic One-Band Hubbard Model, Journal of Superconductivity and Novel Magnetism 30, 85 (2017).
  * Rømer  _et al._ [2015] A. T. Rømer, A. Kreisel, I. Eremin, M. A. Malakhov, T. A. Maier, P. J. Hirschfeld, and B. M. Andersen, Pairing symmetry of the one-band Hubbard model in the paramagnetic weak-coupling limit: A numerical RPA study, Phys. Rev. B 92, 104505 (2015).
  * Graf  _et al._ [2000] M. J. Graf, A. V. Balatsky, and J. A. . Sauls, Local time-reversal-symmetry breaking in ${d}_{{x}^{2}-{y}^{2}}$ superconductors, Phys. Rev. B 61, 3255 (2000).
  * Seibold  _et al._ [2015] G. Seibold, L. Benfatto, C. Castellani, and J. Lorenzana, Amplitude, density, and current correlations of strongly disordered superconductors, Phys. Rev. B 92, 064512 (2015).
  * Mackenzie  _et al._ [1998] A. P. Mackenzie, R. K. W. Haselwimmer, A. W. Tyler, G. G. Lonzarich, Y. Mori, S. Nishizaki, and Y. Maeno, Extremely Strong Dependence of Superconductivity on Disorder in ${\mathrm{Sr}}_{2}{\mathrm{RuO}}_{4}$, Phys. Rev. Lett. 80, 161 (1998).
  * Kikugawa  _et al._ [2003] N. Kikugawa, A. Peter Mackenzie, and Y. Maeno, Effects of In-Plane Impurity Substitution in Sr2RuO4, Journal of the Physical Society of Japan 72, 237 (2003).
  * Kikugawa  _et al._ [2004] N. Kikugawa, A. P. Mackenzie, C. Bergemann, R. A. Borzi, S. A. Grigera, and Y. Maeno, Rigid-band shift of the Fermi level in the strongly correlated metal: ${\mathrm{Sr}}_{2-y}{\mathrm{La}}_{y}{\mathrm{RuO}}_{4}$, Phys. Rev. B 70, 060508 (2004).
  * Deguchi  _et al._ [2004] K. Deguchi, Z. Q. Mao, and Y. Maeno, Determination of the Superconducting Gap Structure in All Bands of the Spin-Triplet Superconductor Sr2RuO4, J. Phys. Soc. Jpn. 73, 1313 (2004).
  * Li  _et al._ [2021b] Y.-S. Li, N. Kikugawa, D. A. Sokolov, F. Jerzembeck, A. S. Gibbs, Y. Maeno, C. W. Hicks, J. Schmalian, M. Nicklas, and A. P. Mackenzie, High-sensitivity heat-capacity measurements on Sr2RuO4 under uniaxial pressure, Proc. Natl. Acad. Sci. USA 118, e2020492118 (2021b).