# Axion-like Particles at Future Neutrino Experiments:
Closing the “Cosmological Triangle”

Vedran Brdar<EMAIL_ADDRESS>Fermi National Accelerator
Laboratory, Batavia, IL, 60510, USA  Northwestern University, Department of
Physics & Astronomy, 2145 Sheridan Road, Evanston, IL 60208, USA  Bhaskar
Dutta<EMAIL_ADDRESS>Mitchell Institute for Fundamental Physics and
Astronomy, Department of Physics and Astronomy, Texas A&M University, College
Station, TX 77843, USA  Wooyoung Jang<EMAIL_ADDRESS>Department of
Physics, University of Texas, Arlington, TX 76019, USA  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  Ian M. Shoemaker<EMAIL_ADDRESS>Center for Neutrino Physics,
Department of Physics, Virginia Tech, Blacksburg, VA 24061, USA  Zahra Tabrizi
<EMAIL_ADDRESS>Center for Neutrino Physics, Department of Physics, Virginia
Tech, Blacksburg, VA 24061, USA  Adrian Thompson<EMAIL_ADDRESS>Mitchell
Institute for Fundamental Physics and Astronomy, Department of Physics and
Astronomy, Texas A&M University, College Station, TX 77843, USA  Jaehoon Yu
<EMAIL_ADDRESS>Department of Physics, University of Texas, Arlington, TX
76019, USA

###### Abstract

Axion-like particles (ALPs) provide a promising direction in the search for
new physics, while a wide range of models incorporate ALPs. We point out that
future neutrino experiments, such as DUNE, possess competitive sensitivity to
ALP signals. The high-intensity proton beam impinging on a target can not only
produce copious amounts of neutrinos, but also cascade photons that are
created from charged particle showers stopping in the target. Therefore, ALPs
interacting with photons can be produced (often energetically) with high
intensity via the Primakoff effect and then leave their signatures at the near
detector through the inverse Primakoff scattering or decays to a photon pair.
Moreover, the high-capability near detectors allow for discrimination between
ALP signals and potential backgrounds, improving the signal sensitivity
further. We demonstrate that a DUNE-like detector can explore a wide range of
parameter space in ALP-photon coupling $g_{a\gamma}$ vs ALP mass $m_{a}$,
including some regions unconstrained by existing bounds; the “cosmological
triangle” will be fully explored and the sensitivity limits would reach up to
$m_{a}\sim 3-4$ GeV and down to $g_{a\gamma}\sim 10^{-8}~{}{\rm GeV}^{-1}$.

††preprint: MI-TH-2029††preprint: nuhep-th/20-13

Introduction. Axions not only address the strong CP problem Peccei and Quinn
(1977); Wilczek (1978); Weinberg (1978), but also provide an explanation for
the 27% of the Universe’s energy content constituting dark matter Preskill
_et al._ (1983); Abbott and Sikivie (1983); Dine and Fischler (1983); Duffy
and van Bibber (2009); Marsh (2016); Battaglieri  _et al._ (2017).
Investigation of the QCD axion has been extended to incorporate general axion-
like particles (ALPs) in a wide range of models. The experimental effort in
the search for ALPs in these many incarnations is vigorously made using their
couplings with Standard Model (SM) particles, primarily to photons, electrons,
and nucleons. These experiments include helioscopes: CAST Zioutas  _et al._
(1999); Anastassopoulos  _et al._ (2017); Irastorza  _et al._ (2013),
haloscopes: Abracadabra Kahn  _et al._ (2016); Salemi (2019), ADMX Asztalos
_et al._ (2001); Du  _et al._ (2018), CASPEr Jackson Kimball  _et al._ (2017),
HAYSTAC Brubaker  _et al._ (2017); Droster and van Bibber (2019), light-
shining-through-wall experiments: ALPSII Spector (2019), interferometry
Melissinos (2009); DeRocco and Hook (2018): ADBC Liu  _et al._ (2019), DANCE
Obata  _et al._ (2018), current and proposed accelerator-based experiments:
FASER Feng  _et al._ (2018), LDMX Berlin  _et al._ (2019); Åkesson  _et al._
(2018), NA62 Volpe (2019), NA64 Dusaev  _et al._ (2020); Banerjee  _et al._
(2020), SeaQuest Berlin  _et al._ (2018), SHiP Alekhin  _et al._ (2016),
hybrids of beam dump and helioscope approaches: PASSAT Bonivento  _et al._
(2020), reactor experiments: MINER, CONUS etc. Dent  _et al._ (2020a);
Aristizabal Sierra  _et al._ (2020), dark matter experiments: XENON Aprile
_et al._ (2020); Dent  _et al._ (2020b), SuperCDMS Aralis  _et al._ (2020),
PandaX Fu  _et al._ (2017) etc.

The coupling to photons is particularly interesting, as it allows the ALPs to
be copiously produced in accelerator-based neutrino experiments where, along
with the neutrino flux, a high-intensity photon flux is generated from intense
proton beams impinging on a target. Photons emerge from bremsstrahlung and
meson decays, and can convert to ALPs via the Primakoff process in the forward
region. The produced ALPs would then be observed at a detector via the inverse
Primakoff scattering process, or via decays into photon pairs. A complementary
analysis exploiting the axion-gluon coupling for QCD axions is done for the
DUNE near detector Kelly  _et al._ (2020).

In this Letter, we present the sensitivity of the ALP search at future
neutrino experiments utilizing their near-detector facilities, taking DUNE as
a concrete example, while the overall search strategies discussed here are
readily applicable to other similar experiments. The near-detector complex at
DUNE will include three detectors, among which we will employ liquid argon
(LAr) and gaseous argon (GAr) detectors for constraining parameter space of
ALPs.

The DUNE near detectors will be located 574 m downstream of the target hall
where the 120-GeV proton beam impinges on a a $1.5$-m long, segmented,
cylindrical graphite target with a $16$-mm diameter. The photon flux will be
created with much higher energy compared to many ongoing accelerator-based
neutrino experiments, e.g., COHERENT, CCM, JSNS2 etc. where $\mathcal{O}$(GeV)
proton beams are used. The higher energy photons would allow us to probe a
considerable range of ALP parameter space. We demonstrate that a DUNE-like
experiment can explore a wide range of the ALP mass [from $\mathcal{O}(1)$ GeV
to the massless limit] vs the ALP-photon coupling parameter space, covering
some regions where there exist no (laboratory-based) constraints, e.g., the
“cosmological triangle” and the region beyond the current beam-dump limits. An
accurate estimate in the photon flux is essential for a more accurate estimate
in the signal sensitivity. For this purpose, we simulate the photon flux using
the GEANT4 Agostinelli  _et al._ (2003) package as conventional event
generators are unable to handle secondary photon production such as cascade
photons whose contribution is significant Dutta  _et al._ (2020). We also
consider the SM background to estimate more realistic constraints on the
parameter space.

ALP production at the target. In order to investigate the ALP parameter space,
we will focus on a generic model where the ALP field (henceforth denoted by
$a$) can couple to a photon as described by interaction terms in the
Lagrangian of the form

$\mathcal{L}\supset-\dfrac{1}{4}g_{a\gamma}aF_{\mu\nu}\tilde{F}^{\mu\nu}\,,$
(1)

where $F_{\mu\nu}$ and $\tilde{F}_{\mu\nu}$ are the usual field strength
tensor of the SM photon and its dual and where $g_{a\gamma}$ parameterizes the
ALP-photon coupling in the unit of inverse energy. While this coupling appears
at the tree level in DFSZ-type models of the QCD axion, it also arises in
models of ALP dark matter and shows up naturally in string axiverse scenarios.

Figure 1: Tree-level ALP production through the Primakoff process (top left),
ALP scattering through the inverse Primakoff process (top right), and ALP
decays (bottom).

Given the above coupling, ALPs can be produced by the Primakoff scattering
process, $\gamma+A\to a+A$ with $A$ symbolizing the atomic system of interest,
as diagramatically displayed in Fig. 1. The production cross-section
$\sigma_{P}$ in ALP scattering angle
$d\Omega_{a}^{\prime}=\sin\theta_{a}^{\prime}d\theta_{a}^{\prime}d\phi_{a}^{\prime}$
is given by

$\frac{\partial^{2}\sigma_{P}}{\partial\theta_{a}^{\prime}\partial\phi_{a}^{\prime}}=\frac{g_{a\gamma}^{2}\alpha}{8\pi}\frac{p_{a}^{4}}{q^{4}}\sin^{3}\theta_{a}^{\prime}F^{2}(q)\,,$
(2)

where $\alpha$ and $m_{a}$ are the fine structure constant and the ALP mass,
respectively and where
$q^{2}=m_{a}^{2}-2E_{\gamma}(E_{a}-p_{a}\cos\theta_{a}^{\prime})$. Here, form
factor $F^{2}(q)$ encapsulates an important coherent enhancement of the atomic
number $Z^{2}$.

For the high energy photon flux in the target, the relevant coherence length
(and momentum transfer) will be at the nuclear scale. We therefore use the
nuclear form factor, for which we adopt the Helm parameterization

$F_{N}^{2}(q)=Z^{2}\bigg{(}\dfrac{3j_{1}(qR_{1})}{qR_{1}}\bigg{)}^{2}e^{-q^{2}s^{2}}\,,$
(3)

with $s=0.9$ fm and $R_{1}=\sqrt{(1.23A^{1/3}-0.6)^{2}+2.18}$ fm with $A$
being the atomic mass number. The results do not change if we replace the
nuclear form factor by the atomic form factor.

We next discuss how the ALP flux is connected to the photon flux at DUNE. Let
us work in the frame $(x,y,z)$ such that $z$ points along the beam axis and
each photon created in the target lies in the $xz$ plane with an angle
$\theta_{\gamma}$ with respect to the $z$ axis. Then, the conversion event on
an atom at $(0,0,0)$ involves a photon with unit vector
$\hat{\gamma}=(\sin\theta_{\gamma},0,\cos\theta_{\gamma})$ and an ALP, which
is generated at an angle $\theta_{a}^{\prime}$ with respect to the photon
direction according to the angular distribution generated by
$\partial^{2}\sigma_{P}/\partial\theta_{a}^{\prime}\partial\phi_{a}^{\prime}$
as in Eq. (2). The ALP direction in the frame
$(x^{\prime},y^{\prime},z^{\prime})$ where the photon momentum points along
the $z^{\prime}$ axis is

$\displaystyle\hat{a}^{\prime}$
$\displaystyle=(\cos\phi_{a}^{\prime}\sin\theta_{a}^{\prime},\sin\phi_{a}^{\prime}\sin\theta_{a}^{\prime},\cos\theta_{a}^{\prime})\,$
(4)

where $\phi_{a}^{\prime}$ is an azimuthal angle around the $z^{\prime}$ axis.
$\hat{a}^{\prime}$ can be transformed to $\hat{a}$ defined in the unprimed
coordinate frame by a rotation by $-\theta_{\gamma}$ about the $y^{\prime}$
axis:

$\displaystyle\hat{a}=($
$\displaystyle\cos\theta_{\gamma}\cos\phi_{a}^{\prime}\sin\theta_{a}^{\prime}-\cos\theta_{a}^{\prime}\sin\theta_{\gamma},\sin\theta_{a}^{\prime}\sin\phi_{a}^{\prime},$
$\displaystyle\cos\theta_{a}^{\prime}\cos\theta_{\gamma}+\cos\phi_{a}^{\prime}\sin\theta_{a}^{\prime}\sin\theta_{\gamma})\,,$
(5)

from which we find that the polar angle of ALP $\theta_{a}$ measured in the
unprimed frame is

$\theta_{a}=\arccos(\cos\theta_{a}^{\prime}\cos\theta_{\gamma}+\cos\phi_{a}^{\prime}\sin\theta_{a}^{\prime}\sin\theta_{\gamma})\,.$
(6)

Assuming that the face of the detector spans a circular aperture of radius
$r$, a distance $\ell$ from the target for simplicity, we see that the ALP
flux would enter the detector as far as $\theta_{a}$ is less than the detector
opening angle $\theta_{\rm det}=\arctan(r/\ell)$. We then convolve the
differential photon flux in $E_{\gamma}$ and $\theta_{\gamma}$ with the
differential Primakoff production cross-section in Eq. (2):

$\displaystyle\frac{d\Phi_{a}}{dE_{a}}=$
$\displaystyle\int\frac{\partial^{2}\Phi_{\gamma}}{\partial
E_{\gamma}\partial\theta_{\gamma}}\bigg{[}\dfrac{1}{\sigma_{P}+\sigma_{\gamma}}\frac{\partial^{2}\sigma_{P}}{\partial\theta_{a}^{\prime}\partial\phi_{a}^{\prime}}\bigg{]}\delta(E_{a}-E_{\gamma})$
(7) $\displaystyle\times$
$\displaystyle\Theta(\theta_{\text{det}}-\theta_{a})d\phi_{a}^{\prime}d\theta_{a}^{\prime}d\theta_{\gamma}\,.$

The factor of
$\frac{1}{\sigma_{P}+\sigma_{\gamma}}\frac{\partial^{2}\sigma_{P}}{\partial\theta_{a}^{\prime}\partial\phi_{a}^{\prime}}$
is the differential branching fraction of ALP production, which is used to
take into account the fraction of photons that convert into ALPs in the target
versus those that are absorbed through standard interactions such as
photoelectric absorption and pair production. One may find that above an MeV,
the total photon absorption cross-section $\sigma_{\gamma}\sim 1$ barn in
carbon.

Finally, we discuss the generation of photon flux before closing this section.
A large fraction of photons are produced via the decays of mesons such as
$\pi^{0}$ and $\eta$ and radiation off the incoming beam. These contributions
can be estimated using conventional event generators. There are additional
significant contributions, e.g., cascade photons from secondary particles such
as ionized electrons created while final-state particles lose their energy and
stop in the target Dutta  _et al._ (2020). The estimate of those photons
involves non-trivial nuclear effects, so we employ the GEANT4 package for a
systematic estimate of the photon flux, simply considering the target geometry
and the beam specifications of a DUNE-like experiment with a 1.2-MW, 120-GeV
proton beam impinging on a graphite target. We assumed a proton-on-target
(POT) rate of $1.1\times 10^{21}$ POT$\cdot$year-1. We find that most of the
photons are moving in the forward direction, hence one can expect that a large
fraction of ALP flux is directed toward the near detector complex.

ALP event rate at the detector. ALPs can be detected by inverse Primakoff
scattering or decays to $\gamma\gamma$ at the detector.

For the decay to a photon pair, we compute the probability $P_{\rm decay}$
that the ALP decays within the detector volume. This can be calculated by
integrating the decay probability density between the front and the back of
the detector, $(\ell,\ell+\Delta\ell)$;

$P_{\rm decay}=e^{-\ell/(\tau v_{a})}\left[1-e^{-\Delta\ell/(\tau
v_{a})}\right]\,,$ (8)

where $\tau$ is the ALP lifetime in the laboratory frame, based on the decay
width

$\Gamma(a\to\gamma\gamma)=\dfrac{g_{a\gamma}^{2}m_{a}^{3}}{64\pi},$ (9)

and velocity in the laboratory frame $v_{a}=p_{a}/E_{a}$. For a given exposure
time $\mathcal{E}$, the total event rate from ALP decays $N_{\rm decay}$ is
then given by a convolution of the decay probability with the ALP flux derived
in Eq. (7):

$N_{\rm decay}=\mathcal{E}\int\frac{d\Phi_{a}}{dE_{a}}P_{\rm decay}dE_{a}\,.$
(10)

Similarly, ALPs can scatter, yielding a single $\gamma$ signal, through the
inverse Primakoff scattering process. This process has a cross-section
$\sigma_{IP}$, identical to $\sigma_{P}$ but larger by a factor of 2 to
account for the sum over polarization states.

$N_{\rm
scatter}=N_{T}\mathcal{E}\int\sigma_{IP}(E_{a})\frac{d\Phi_{a}}{dE_{a}}P_{\rm
surv}dE_{a}\,,$ (11)

where $N_{T}$ is the number of argon targets. $P_{\rm surv}$ is the survival
probability that an ALP reaches the detector without decaying, that is,
$P_{\rm surv}=e^{-\ell/(\tau v_{a})}$.

Background consideration. In Fig. 3, one can infer two qualitatively different
exclusion regions corresponding to scatterings and decays of ALPs. The limits
will be discussed in more detail in the next section and here we note that
they were obtained under a zero-background hypothesis. In what follows, we
justify such a choice for the case of decaying ALPs in a DUNE-like GAr
detector and also discuss the situation in a DUNE-like LAr detector for
scattering ALPs. For a DUNE-like GAr detector, we make use of its volume and
consider ALP decays into two photons. The main SM background is from
NC-$\pi^{0}$ production where the final-state pion decays also to two photons.

For the 7-year exposure, the DUNE Collaboration predicts $\sim 3\times 10^{6}$
NC events in the detector Carneiro (2020), summing up the contributions from
positive and negative horn polarity modes. Let us conservatively assume that
all of those would contain a single pion. Among those, the events where there
is some hadronic activity in the interaction vertex could be vetoed as ALPs
decay would not give rise to additional energetic charged particles. The study
in Berryman  _et al._ (2020) estimated that the cut on the hadronic activity
removes $\sim 80\%$ of the background.

For a further background suppression, one can employ kinematics: while the
ALPs and hence its decay products (photons) will appear in the detector at a
very small opening angle (see also the left panel of Fig. 2) and are
practically collimated with the beamline, the $\pi^{0}$s in the detector (and
hence photons from their decay) would be distributed with a greater preference
than our signal toward larger angles. We have checked that explicitly, using
Monte Carlo generator GENIE Andreopoulos  _et al._ (2015). We found that by
imposing a cut on the angle $\epsilon=4$ mrad, which matches resolution of the
detector Berryman  _et al._ (2020), background is reduced by the factor of
$10^{3}$.

For the coherent single photon production (not yet implemented in GENIE), we
expect negligible contributions, considering its 3-4 orders of magnitude
smaller cross-sections than those of NC-$\pi^{0}$ events Rein and Sehgal
(1981); Wang  _et al._ (2014). For completeness, let us emphasize that the
background from misidentified electrons is suppressed; the $\nu_{e}$ flux
producing electrons in CC processes is suppressed with respect to the
$\nu_{\mu}$ flux. Further, the dirt events which could be present in both LAr
and GAr analyses are subdominant with respect to the in-detector $\pi^{0}$
background; this is mainly because to generate dirt events, neutrinos would
need to interact very close to the upstream end of the LAr detector and this
dramatically reduces effective volume of the dirt. The few events that may
sneak into the detector can be removed using topological and spatial cuts. See
the related discussion for MiniBooNE Aguilar-Arevalo  _et al._ (2021).

Figure 2: Left: Opening angle $\Delta\theta_{\gamma\gamma}$ from
$a\to\gamma\gamma$ decays for representative ALP masses in the decay-dominated
signal regime for a fixed lifetime. Middle: Decaying ALP fluxes as a function
of the angle with respect to the beam axis. The black dotted line represents
the geometric acceptance for the GAr-detector solid angle. Right: Normalized
energy spectra for the decay-dominated ($a\to\gamma\gamma$) region of
parameter space for a fixed ALP lifetime.

Finally, the likelihood for a photon conversion in the gas is only 12%, which
means that there is 23% chance that at least one of the photons would convert.
This reduces the background counts by a further factor of 4, but we note that
our signal also suffers from this. Putting all the numbers together leaves us
with $\mathcal{O}(10)$ background events. We note that this number could be
further reduced by employing kinematic variables (which depend on $m_{a}$).
For example, it is possible to use the opening angle,
$\Delta\theta_{\gamma\gamma}$, between the photons (the left panel of Fig. 2)
and the photon energy spectra (the right panel of Fig. 2) emerging from the
ALP decays for various $m_{a}$ to suppress the backgrounds. Still, note that
by even taking into account background of this level would basically not yield
an observable change of the presented limits.

As far as the limit from ALP scattering in LAr is concerned, here the
situation is vastly different. First, the fiducial mass of LAr detector is 50
times larger than that of GAr, hence giving $\mathcal{O}(10^{8})$ $\pi^{0}$
events before cuts. We anticipate that most of those would typically yield two
observable showers and hence would be rejected since, in contrast, there would
be only one final-state photon from the ALP scattering. For eliminating the
remaining $\pi^{0}$ events where one of the photons would be missed, one can
employ kinematical distributions. Regarding those, unlike in the decaying
case, one cannot assume that all the photons from ALPs will be within an
$\epsilon$ angle because photons emerging from ALP scattering in the detector
would have distribution similar to those arising from decaying $\pi^{0}$. This
makes the background consideration more involved than in the case of decaying
ALPs. We reserve a more quantitative study in this regard for a follow-up work
and present a background-free limit in what follows.

Results. We are now in the position to discuss the expected sensitivity of a
DUNE-like detector to the ALP scattering and decay signals in the presence of
SM backgrounds. We consider the LAr detector of a 50-ton fiducial mass for the
scattering limit and the GAr detector of a cylindrical fiducial volume (5.2-m
diameter and 5-m length) for the decay limit. Our sensitivity reach for decays
are then evaluated at 90% C.L. and reported in Fig. 3 by the red lines, while
the reach using inverse Primakoff scattering is shown in blue. Various
existing astrophysical and laboratory-based limits are also shown by tan and
gray regions, respectively. For lower ALP masses, notice that the scattering
limit becomes much better than the decay limit near $m_{a}\sim 100$ keV.

This study suggests that a DUNE-like detector, especially, the gas-phase
detector due to its significantly low background contamination, will be
greatly capable of probing a wide range of the $(m_{a},g_{a\gamma})$ parameter
space. A few remarkable potentials are worth pointing out. Above all, a DUNE-
like detector can probe the regions beyond the current beam-dump limits, i.e.,
up to $m_{a}\sim 3-4$ GeV and down to $g_{a\gamma}\sim 10^{-8}~{}{\rm
GeV}^{-1}$, which have never been explored before by existing direct or
indirect searches. Second, the cosmological triangle surrounded by beam-dump,
HB stars, and supernova limits can be (completely) covered. Although this
region would be constrained by standard cosmological considerations, they are
highly model-dependent, hence can be evaded by non-standard cosmology Carenza
_et al._ (2020); Depta  _et al._ (2020). Third, a large portion of
astrophysical limits (HB stars and SN1987a) can be constrained (see also
Lucente  _et al._ (2020) for a recent development in the supernova limit
calculation and Bar  _et al._ (2020) for its criticism). Like the cosmological
limits, the astrophysical limits highly depend on the underlying ALP model
details (see, for example, Refs. Jaeckel  _et al._ (2007); Khoury and Weltman
(2004); Masso and Redondo (2005, 2006); Dupays  _et al._ (2007); Mohapatra and
Nasri (2007); Brax  _et al._ (2007); DeRocco  _et al._ (2020)), so a DUNE-like
detector is expected to constrain those regions in a model-independent
fashion. A number of ongoing/upcoming experiments, e.g., Belle II and neutrino
experiments having smaller beam energy (e.g., CCM, SBND, T2K/T2HK) would also
be able to investigate the cosmological triangle using the similar production
process as described here but DUNE is going to have much better reach for
smaller $g_{a\gamma}$/larger $m_{a}$ values.

Figure 3: The 90% C.L. sensitivity reaches in the $(m_{a},g_{a\gamma})$ plane
expected for a DUNE-like GAr detector with 1 and 7 year exposures
($a\to\gamma\gamma$ decays) and at a LAr detector for 1 and 7 year exposures
(inverse Primakoff scattering, $a+A\to\gamma+A$). Existing laboratory limits
Blumlein  _et al._ (1992, 1991); Jaeckel and Spannowsky (2016) are shown in
gray, while astrophysical limits Raffelt (1986); Raffelt and Dearborn (1987);
Raffelt (2008); Payez  _et al._ (2015); Jaeckel  _et al._ (2018) are shown in
tan.

Outlook. We plan to extend our analysis in several directions in future work.
For example, DUNE-PRISM may offer advantages for further background
discrimination given its ability to move several degrees off-axis. We have
also not considered a detailed background estimation for the scattering regime
relevant for ALPs below 0.1 MeV. Lastly, we expect other neutrino facilities
to play a complementary role in the future hunt for ALPs.

###### Acknowledgements.

We would like to thank Alan Bross, André de Gouvêa, Laura Fields, Vladimir
Ivantchenko, Soon Yung Jun and Shirley Li for useful discussions. We also
thank Patrick Huber to be part of this work in the initial phase. JY also
acknowledges discussions with several other members of the DUNE Collaboration.
BD and AT acknowledge support from the U.S. Department of Energy (DOE) Grant
DE-SC0010813. The work of DK is supported by DOE under Grant No. DE-
FG02-13ER41976/DE-SC0009913/DE-SC0010813. The work of IMS is supported by DOE
under the award number DE-SC0020250. The work of ZT is supported by DOE under
the award numbers DE-SC0020250 and DE-SC0020262. WJ and JY acknowledge the
support from DOE under Grant No. DE-SC0011686. Fermilab is operated by the
Fermi Research Alliance, LLC under contract No. DE-AC02-07CH11359 with the
United States Department of Energy. Portions of this research were conducted
with the advanced computing resources provided by Texas A&M High Performance
Research Computing.

## References

  * Peccei and Quinn (1977) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Preskill  _et al._ (1983) J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Abbott and Sikivie (1983) L. Abbott and P. Sikivie, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler (1983) M. Dine and W. Fischler, Phys. Lett. B 120, 137 (1983).
  * Duffy and van Bibber (2009) L. D. Duffy and K. van Bibber, New J. Phys. 11, 105008 (2009), arXiv:0904.3346 [hep-ph] .
  * Marsh (2016) D. J. E. Marsh, Phys. Rept. 643, 1 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Battaglieri  _et al._ (2017) M. Battaglieri _et al._ , in  _U.S. Cosmic Visions: New Ideas in Dark Matter College Park, MD, USA, March 23-25, 2017_ (2017) arXiv:1707.04591 [hep-ph] .
  * Zioutas  _et al._ (1999) K. Zioutas _et al._ , Nucl. Instrum. Meth. A425, 480 (1999), arXiv:astro-ph/9801176 [astro-ph] .
  * Anastassopoulos  _et al._ (2017) V. Anastassopoulos _et al._ (CAST), Nature Phys. 13, 584 (2017), arXiv:1705.02290 [hep-ex] .
  * Irastorza  _et al._ (2013) I. Irastorza _et al._ (IAXO), (2013).
  * Kahn  _et al._ (2016) Y. Kahn, B. R. Safdi, and J. Thaler, Phys. Rev. Lett. 117, 141801 (2016), arXiv:1602.01086 [hep-ph] .
  * Salemi (2019) C. P. Salemi (ABRACADABRA), in  _54th Rencontres de Moriond on Electroweak Interactions and Unified Theories (Moriond EW 2019) La Thuile, Italy, March 16-23, 2019_ (2019) arXiv:1905.06882 [hep-ex] .
  * Asztalos  _et al._ (2001) S. J. Asztalos _et al._ (ADMX), Phys. Rev. D 64, 092003 (2001).
  * Du  _et al._ (2018) N. Du _et al._ (ADMX), Phys. Rev. Lett. 120, 151301 (2018), arXiv:1804.05750 [hep-ex] .
  * Jackson Kimball  _et al._ (2017) D. F. Jackson Kimball _et al._ , (2017), arXiv:1711.08999 [physics.ins-det] .
  * Brubaker  _et al._ (2017) B. M. Brubaker _et al._ , Phys. Rev. Lett. 118, 061302 (2017), arXiv:1610.02580 [astro-ph.CO] .
  * Droster and van Bibber (2019) A. Droster and K. van Bibber (HAYSTAC), in  _13th Conference on the Intersections of Particle and Nuclear Physics (CIPANP 2018) Palm Springs, California, USA, May 29-June 3, 2018_ (2019) arXiv:1901.01668 [physics.ins-det] .
  * Spector (2019) A. Spector (ALPS), in  _14th Patras Workshop on Axions, WIMPs and WISPs (AXION-WIMP 2018) (PATRAS 2018) Hamburg, Germany, June 18-22, 2018_ (2019) arXiv:1906.09011 [physics.ins-det] .
  * Melissinos (2009) A. C. Melissinos, Phys. Rev. Lett. 102, 202001 (2009), arXiv:0807.1092 [hep-ph] .
  * DeRocco and Hook (2018) W. DeRocco and A. Hook, Phys. Rev. D 98, 035021 (2018), arXiv:1802.07273 [hep-ph] .
  * Liu  _et al._ (2019) H. Liu, B. D. Elwood, M. Evans, and J. Thaler, Phys. Rev. D 100, 023548 (2019), arXiv:1809.01656 [hep-ph] .
  * Obata  _et al._ (2018) I. Obata, T. Fujita, and Y. Michimura, Phys. Rev. Lett. 121, 161301 (2018), arXiv:1805.11753 [astro-ph.CO] .
  * Feng  _et al._ (2018) J. L. Feng, I. Galon, F. Kling, and S. Trojanowski, Phys. Rev. D 98, 055021 (2018), arXiv:1806.02348 [hep-ph] .
  * Berlin  _et al._ (2019) A. Berlin, N. Blinov, G. Krnjaic, P. Schuster, and N. Toro, Phys. Rev. D 99, 075001 (2019), arXiv:1807.01730 [hep-ph] .
  * Åkesson  _et al._ (2018) T. Åkesson _et al._ (LDMX), (2018), arXiv:1808.05219 [hep-ex] .
  * Volpe (2019) R. Volpe, in  _Meeting of the Division of Particles and Fields of the American Physical Society (DPF2019) Boston, Massachusetts, July 29-August 2, 2019_ (2019) arXiv:1910.10429 [hep-ex] .
  * Dusaev  _et al._ (2020) R. Dusaev, D. Kirpichnikov, and M. Kirsanov, Phys. Rev. D 102, 055018 (2020), arXiv:2004.04469 [hep-ph] .
  * Banerjee  _et al._ (2020) D. Banerjee _et al._ (NA64), Phys. Rev. Lett. 125, 081801 (2020), arXiv:2005.02710 [hep-ex] .
  * Berlin  _et al._ (2018) A. Berlin, S. Gori, P. Schuster, and N. Toro, Phys. Rev. D 98, 035011 (2018), arXiv:1804.00661 [hep-ph] .
  * Alekhin  _et al._ (2016) S. Alekhin _et al._ , Rept. Prog. Phys. 79, 124201 (2016), arXiv:1504.04855 [hep-ph] .
  * Bonivento  _et al._ (2020) W. M. Bonivento, D. Kim, and K. Sinha, Eur. Phys. J. C 80, 164 (2020), arXiv:1909.03071 [hep-ph] .
  * Dent  _et al._ (2020a) J. B. Dent, B. Dutta, D. Kim, S. Liao, R. Mahapatra, K. Sinha, and A. Thompson, Phys. Rev. Lett. 124, 211804 (2020a), arXiv:1912.05733 [hep-ph] .
  * Aristizabal Sierra  _et al._ (2020) D. Aristizabal Sierra, V. De Romeri, L. Flores, and D. Papoulias, (2020), arXiv:2010.15712 [hep-ph] .
  * Aprile  _et al._ (2020) E. Aprile _et al._ (XENON), (2020), arXiv:2006.09721 [hep-ex] .
  * Dent  _et al._ (2020b) J. B. Dent, B. Dutta, J. L. Newstead, and A. Thompson, (2020b), arXiv:2006.15118 [hep-ph] .
  * Aralis  _et al._ (2020) T. Aralis _et al._ (SuperCDMS Collaboration), Phys. Rev. D 101, 052008 (2020).
  * Fu  _et al._ (2017) C. Fu, X. Zhou, X. Chen, Y. Chen, X. Cui, D. Fang, K. Giboni, F. Giuliani, K. Han, X. Huang, and et al., Physical Review Letters 119 (2017), 10.1103/physrevlett.119.181806.
  * Kelly  _et al._ (2020) K. J. Kelly, S. Kumar, and Z. Liu, (2020), arXiv:2011.05995 [hep-ph] .
  * Agostinelli  _et al._ (2003) S. Agostinelli _et al._ (GEANT4), Nucl. Instrum. Meth. A506, 250 (2003).
  * Dutta  _et al._ (2020) B. Dutta, D. Kim, S. Liao, J.-C. Park, S. Shin, L. E. Strigari, and A. Thompson, (2020), arXiv:2006.09386 [hep-ph] .
  * Carneiro (2020) M. F. Carneiro, “Neutrino Interaction Physics & the DUNE Near Detector,” https://indico.cern.ch/event/868940/contributions/3816980/attachments/2082771/3498525/XSec_and_the_DUNE_ND_-_Mateus_F._Carneiro_2.pdf (2020).
  * Berryman  _et al._ (2020) J. M. Berryman, A. de Gouvea, P. J. Fox, B. J. Kayser, K. J. Kelly, and J. L. Raaf, JHEP 02, 174 (2020), arXiv:1912.07622 [hep-ph] .
  * Andreopoulos  _et al._ (2015) C. Andreopoulos, C. Barry, S. Dytman, H. Gallagher, T. Golan, R. Hatcher, G. Perdue, and J. Yarba, “The genie neutrino monte carlo generator: Physics and user manual,” (2015), arXiv:1510.05494 [hep-ph] .
  * Rein and Sehgal (1981) D. Rein and L. M. Sehgal, Phys. Lett. B 104, 394 (1981), [Erratum: Phys.Lett.B 106, 513 (1981)].
  * Wang  _et al._ (2014) E. Wang, L. Alvarez-Ruso, and J. Nieves, Phys. Rev. C 89, 015503 (2014), arXiv:1311.2151 [nucl-th] .
  * Aguilar-Arevalo  _et al._ (2021) A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. D 103, 052002 (2021), arXiv:2006.16883 [hep-ex] .
  * Carenza  _et al._ (2020) P. Carenza, O. Straniero, B. Döbrich, M. Giannotti, G. Lucente, and A. Mirizzi, Phys. Lett. B 809, 135709 (2020), arXiv:2004.08399 [hep-ph] .
  * Depta  _et al._ (2020) P. F. Depta, M. Hufnagel, and K. Schmidt-Hoberg, JCAP 05, 009 (2020), arXiv:2002.08370 [hep-ph] .
  * Lucente  _et al._ (2020) G. Lucente, P. Carenza, T. Fischer, M. Giannotti, and A. Mirizzi, (2020), arXiv:2008.04918 [hep-ph] .
  * Bar  _et al._ (2020) N. Bar, K. Blum, and G. D’Amico, Phys. Rev. D 101, 123025 (2020), arXiv:1907.05020 [hep-ph] .
  * Jaeckel  _et al._ (2007) J. Jaeckel, E. Masso, J. Redondo, A. Ringwald, and F. Takahashi, Phys. Rev. D 75, 013004 (2007), arXiv:hep-ph/0610203 [hep-ph] .
  * Khoury and Weltman (2004) J. Khoury and A. Weltman, Phys. Rev. Lett. 93, 171104 (2004), arXiv:astro-ph/0309300 [astro-ph] .
  * Masso and Redondo (2005) E. Masso and J. Redondo, JCAP 0509, 015 (2005), arXiv:hep-ph/0504202 [hep-ph] .
  * Masso and Redondo (2006) E. Masso and J. Redondo, Phys. Rev. Lett. 97, 151802 (2006), arXiv:hep-ph/0606163 [hep-ph] .
  * Dupays  _et al._ (2007) A. Dupays, E. Masso, J. Redondo, and C. Rizzo, Phys. Rev. Lett. 98, 131802 (2007), arXiv:hep-ph/0610286 [hep-ph] .
  * Mohapatra and Nasri (2007) R. N. Mohapatra and S. Nasri, Phys. Rev. Lett. 98, 050402 (2007), arXiv:hep-ph/0610068 [hep-ph] .
  * Brax  _et al._ (2007) P. Brax, C. van de Bruck, and A.-C. Davis, Phys. Rev. Lett. 99, 121103 (2007), arXiv:hep-ph/0703243 [HEP-PH] .
  * DeRocco  _et al._ (2020) W. DeRocco, P. W. Graham, and S. Rajendran, (2020), arXiv:2006.15112 [hep-ph] .
  * Blumlein  _et al._ (1992) J. Blumlein _et al._ , Int. J. Mod. Phys. A 7, 3835 (1992).
  * Blumlein  _et al._ (1991) J. Blumlein _et al._ , Z. Phys. C 51, 341 (1991).
  * Jaeckel and Spannowsky (2016) J. Jaeckel and M. Spannowsky, Phys. Lett. B 753, 482 (2016), arXiv:1509.00476 [hep-ph] .
  * Raffelt (1986) G. G. Raffelt, Phys. Rev. D 33, 897 (1986).
  * Raffelt and Dearborn (1987) G. G. Raffelt and D. S. Dearborn, Phys. Rev. D 36, 2211 (1987).
  * Raffelt (2008) G. G. Raffelt, Lect. Notes Phys. 741, 51 (2008), arXiv:hep-ph/0611350 .
  * Payez  _et al._ (2015) A. Payez, C. Evoli, T. Fischer, M. Giannotti, A. Mirizzi, and A. Ringwald, JCAP 02, 006 (2015), arXiv:1410.3747 [astro-ph.HE] .
  * Jaeckel  _et al._ (2018) J. Jaeckel, P. Malta, and J. Redondo, Phys. Rev. D 98, 055032 (2018), arXiv:1702.02964 [hep-ph] .


#  Raman Study of Cooper Pairing Instabilities in (Li1-xFex)OHFeSe

G. He<EMAIL_ADDRESS>Walther Meissner Institut, Bayerische Akademie der
Wissenschaften, 85748 Garching, Germany  D. Li  Beijing National Laboratory
for Condensed Matter Physics, Institute of Physics, Chinese Academy of
Sciences, Beijing 100190, China.  School of Physical Sciences, University of
Chinese Academy of Sciences, Beijing 100049, China.  D. Jost  Walther Meissner
Institut, Bayerische Akademie der Wissenschaften, 85748 Garching, Germany
Fakultät für Physik E23, Technische Universität München, 85748 Garching,
Germany  A. Baum  Walther Meissner Institut, Bayerische Akademie der
Wissenschaften, 85748 Garching, Germany  P. P. Shen  Beijing National
Laboratory for Condensed Matter Physics, Institute of Physics, Chinese Academy
of Sciences, Beijing 100190, China.  School of Physical Sciences, University
of Chinese Academy of Sciences, Beijing 100049, China.  X. L. Dong  Beijing
National Laboratory for Condensed Matter Physics, Institute of Physics,
Chinese Academy of Sciences, Beijing 100190, China.  School of Physical
Sciences, University of Chinese Academy of Sciences, Beijing 100049, China.
Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China  Z. X.
Zhao  Beijing National Laboratory for Condensed Matter Physics, Institute of
Physics, Chinese Academy of Sciences, Beijing 100190, China.  School of
Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049,
China.  Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China
R. Hackl<EMAIL_ADDRESS>Walther Meissner Institut, Bayerische Akademie der
Wissenschaften, 85748 Garching, Germany  Fakultät für Physik E23, Technische
Universität München, 85748 Garching, Germany

(August 0 d , 2024)

###### Abstract

We studied the electronic Raman spectra of (Li1-xFex)OHFeSe as a function of
light polarization and temperature. In the ${B_{1g}}$ spectra alone we observe
the redistribution of spectral weight expected for a superconductor and two
well-resolved peaks below $T_{c}$. The nearly resolution-limited peak at 110
cm-1 (13.6 meV) is identified as a collective mode. The peak at 190 cm-1 (23.6
meV) is presumably another collective mode since the line is symmetric and its
energy is significantly below the gap energy observed by single-particle
spectroscopies. Given the experimental band structure of (Li1-xFex)OHFeSe, the
most plausible explanations include conventional spin-fluctuation pairing
between the electron bands and the incipient hole band and pairing between the
hybridized electron bands. The absence of gap features in ${A_{1g}}$ and
${B_{2g}}$ symmetry favors the second case. Thus, in spite of various
differences between the pnictides and chalcogenides, this Letter demonstrates
the proximity of pairing states and the importance of band structure effects
in the Fe-based compounds.

The mechanism of Cooper pairing in Fe-based superconductors (FeSCs) or the
copper-oxygen compounds is among the most vexing problems in condensed matter
physics. The at least partial understanding of these unconventional
superconductors would pave the way toward new materials. In either case
superconductivity occurs close to magnetic order Paglione and Greene (2010).
Consequently, spin fluctuations are among the candidates for supporting
electron pairing Mazin  _et al._ (2008); Scalapino (2012). Alternatively,
charge Onari and Kontani (2009) or orbital fluctuations Kontani and Onari
(2010) between the Fe $3d$ orbitals, spin-orbit coupling Borisenko  _et al._
(2016) and/or nematic fluctuations Lederer  _et al._ (2015) may support Cooper
pairing. In all cases the Fermi surface topology strongly influences the
pairing tendencies and qualitative differences between the pnictides and
chalcogenides may be expected and were scrutinized in doped BaFe2As2 (122) and
FeSe-based (11) compounds.

Intercalated FeSe superconductors show $T_{c}$ values higher than 40 K (refs.
Lu  _et al._ (2015); Burrard-Lucas  _et al._ (2013)) but a Fermi surface
topology different from the pnictides. In (Li1-xFex)OHFeSe, as shown in Fig. 1
(a2), the holelike Fermi surface encircling the $\Gamma$ point in 122
compounds and bulk FeSe cannot be resolved in angle-resolved photoemission
spectroscopy (ARPES) any further Zhao  _et al._ (2016) while it is still
present in density functional theoryNekrasov and Sadovskii (2015) marking one
of the similarities between intercalated and monolayer FeSe Zhao  _et al._
(2016); Shi  _et al._ (2017). This similarity triggers the question as to the
pairing interactions and the directly related gap structure. Neither the
recent ARPES nor the tunneling experiments yielded clear answers here but show
only that there are essentially two rather different gap energies Du  _et al._
(2016, 2018); Chen  _et al._ (2019) on the presumably hybridized concentric
electronlike Fermi surfaces Zhao  _et al._ (2016).

Here electronic Raman scattering can contribute useful information Devereaux
and Hackl (2007); Scalapino and Devereaux (2009); Maiti  _et al._ (2016);
Kretzschmar  _et al._ (2013); Böhm  _et al._ (2014). In addition to the gap
formation and the pair-breaking peaks at approximately twice the gap energy
Abrikosov and Fal’kovskii (1961); Klein and Dierker (1984), collective
excitations appear in the Raman response which are related to details of the
pairing potential $V_{{\bf k}{\bf k}^{\prime}}$. Collective excitations in
superconductors were first discussed by Bardasis and Schrieffer (BS) Bardasis
and Schrieffer (1961) and by Leggett Leggett (1966). The BS mode stems from a
subleading pairing interaction that is orthogonal to the ground state. The
Leggett mode is best thought of as interband Josephson-like number-phase
fluctuation, the absolute energy of which corresponds to the relative coupling
strength between the bands in comparison to the intraband coupling Klein
(2010); Blumberg  _et al._ (2007). Thus, the careful study of putative
collective modes offers an opportunity to clarify the competing
superconducting instabilities and the related pairing glue.

Figure 1: (a1) 1 Fe (dashes) and 2 Fe (solid line) unit cells. (a2) The
crystal structure of (Li1-xFex)OHFeSe. (b1)-(b3) Polarization configurations
of ${A_{1g}}$, ${B_{1g}}$ and ${B_{2g}}$ symmetries. The polarization
configurations are indicated on the FeSe-layer. (c1)-(c3) First- and (c4)-(c6)
second-order Raman vertices in the first Brillouin Zone (BZ, 1 Fe unit cell)
of ${A_{1g}}$, ${B_{1g}}$ and ${B_{2g}}$ symmetry, respectively, for the
$D_{4h}$ point group. The unfolded electron pockets are shown as red half
ellipses.

In this Letter, we present polarization-dependent Raman spectra for
temperatures between 7.2 K and 48 K in high-quality single-crystalline
(Li1-xFex)OHFeSe thin films with $T_{c}=42$ K. At 7.2 K we observe two well-
defined features at 110 and 190 cm-1 in ${B_{1g}}$ but not in ${A_{1g}}$ and
${B_{2g}}$ symmetry. We conclude that at least the resolution-limited line at
110 cm-1 is a collective mode being either related to a subleading pairing
interaction or a number-phase oscillation between the electron bands. The
superconducting ground state in (Li1-xFex)OHFeSe may result from either spin
fluctuations between the electron bands and the incipient hole band or the
interaction between the hybridized electron bands.

(Li1-xFex)OHFeSe ($x\sim 0.18$) thin films were grown epitaxially on
$(00l)$-oriented LaAlO3 substrates as reported previously Huang  _et al._
(2017a, b). The thin films have a typical thickness of 100 nm and were
characterized by x-ray diffraction and magnetization measurements showing a
high crystallinity and a superconducting transition at $T_{c}=42\pm 1$ K (see
Supplemental Material A). The Raman experiments were carried out with a
standard light scattering equipment Devereaux and Hackl (2007). For excitation
we used a solid-state and an Ar+ laser emitting at 577 and 457 nm,
respectively. All spectra were measured with an absorbed laser power of
$P_{\rm abs}=2$ mW limiting the heating in the spot to below 1.5 K/mW (see
Supplemental Material B). The polarizations of the incoming and scattered
photons will be defined with respect to the 1 Fe unit cell as shown in Figs.
1(b1)- 1(b3), which are more appropriate for electronic excitations. We show
Raman susceptibilities
$R\chi^{\prime\prime}(T,\Omega)=S(\Omega,T)\\{1+n(T,\Omega)\\}^{-1}$ , where
$R$ is an experimental constant, $S(\Omega,T)$ is the dynamical structure
factor that is proportional to the rate of scattered photons, and
$n(T,\Omega)$ is the Bose-Einstein distribution function. The first- and
second-order crystal harmonics of each symmetry which Raman vertices are
proportional to and the position of the Fermi pockets of (Li1-xFex)OHFeSe are
shown in Figs. 1(c1)- 1(c6) to illustrate the relation between electronic
Raman response in different symmetries and the Fermi surface topology.

Figure 2 shows the polarization-dependent Raman response of (Li1-xFex)OHFeSe
above (red) and below (blue) $T_{c}$. The peaks observed at approximately 167
cm-1 in ${A_{1g}}$ symmetry [Fig. 2 (e)] and 205 cm-1 in ${B_{2g}}$ symmetry
[Fig. 2 (c)] correspond to Se in-phase and Fe out-of-phase vibrations along
the $c$ axis, respectively Zhang  _et al._ (2019). The peak at 165 cm-1 in the
${B_{1g}}$ spectra [Fig. 2 (a)] can be identified as a leakage from the
${A_{1g}}$ phonon. The small line width and high intensity of the phonon lines
underpin the excellent crystalline quality of the sample. Because of surface
contamination, broad peaks at 240 cm-1 [Fig. 2 (a)], 224 and 256 cm-1 [Fig. 2
(c)], and 252 cm-1 [Fig. 2 (e)] in ${B_{1g}}$, ${B_{2g}}$ and ${A_{1g}}$
symmetry, respectively, appear in the spectra. These peaks probably originate
from Fe oxide phonons since they disappear after cleaving and are discussed in
more detail in Supplemental Material C.

Figure 2: (a),(c),(e) Raman spectra of (Li1-xFex)OHFeSe in ${B_{1g}}$,
${B_{2g}}$, and ${A_{1g}}$ symmetry. (b),(d),(f) Difference spectra between
superconducting and normal state.

If the normal state spectra are subtracted from the superconducting spectra
all phonons and extra lines disappear since they do not change appreciably
upon crossing $T_{c}$ as shown in Figs. 2 (b), 2(d), and 2(f). However, due to
the high spectral resolution even small changes of the phonon lines can be
identified.

The small changes in the phonon lines are the only detectable effects of
superconductivity in ${A_{1g}}$ and ${B_{2g}}$ symmetry. The increase toward
zero energy is an artifact resulting from insufficient rejection of the laser
line. The spectral changes in ${B_{1g}}$ symmetry are resolved clearly since
the extra peaks at approximately 110 cm-1 and 190 cm-1 have an intensity
comparable to that of the phonons. In addition to the peaks, the continuum is
suppressed below 90 cm-1 and is nearly energy independent. None of the
excitations display appreciable resonance behavior, and the data can be
reproduced in different regions of the sample, as shown in Supplemental
Material D and E. The suppression and the additional peaks indicate a relation
to superconductivity.

The most important observations include the following:
(i) There is no intensity redistribution below $T_{c}$ in ${A_{1g}}$ and
${B_{2g}}$ symmetry, as observed earlierKretzschmar  _et al._ (2013); Muschler
_et al._ (2009). The phenomenon can be understood qualitatively in terms of
the related polarization-dependent Raman form factors, the Fermi surface
topology of (Li1-xFex)OHFeSe as shown in Figs. 1(c1) and 1(c3), and screening
effects.
(ii) The${B_{1g}}$ spectrum is suppressed below 90 cm-1. This indicates a
nearly isotropic superconducting gap as already observed for ${\rm
Rb_{0.8}Fe_{1.6}Se_{2}}$ Kretzschmar  _et al._ (2013). The residual intensity
of approximately 2-3 counts(s mW)-1 is not entirely clear but may either
originate from the substrate, surface layers, or luminescence. In agreement
with the scanning tunneling spectroscopy (STS) results Du  _et al._ (2016),
there is no reason to assume that there are states inside the gap. Similar
residual Raman intensities are also observed in single crystals of pnictides
and other chalcogenides Kretzschmar  _et al._ (2013); Böhm  _et al._ (2014);
Jost  _et al._ (2018) and can safely be assumed to be extrinsic.
(iii) There are two superconductivity-induced features separated by some 80
cm-1. Gap features at a similar separation but slightly higher energies were
observed by STS and ARPES as summarized in Table 1. The significant
differences in the derived gap energies presumably have a real physical reason
such as a substantial energy difference between the gap and collective modes
Böhm  _et al._ (2018). The line at 190 cm-1 is nearly symmetric and
10-12$\,\mathrm{cm}^{-1}$ wide while that at 110 cm-1 is rather sharp around
the maximum but asymmetric.
(iv) The peak at 110$\,\mathrm{cm}^{-1}$ is nearly resolution-limited. For its
width, it cannot result from pair breaking alone. Rather the symmetric line at
110 cm-1 having a FWHM of less than 5 cm-1 is a collective mode while the
shoulder on the high-energy side originates from pair breaking. We will
explore this possibility later by a phenomenological analysis.

Figure 3: (a) Temperature dependence of the difference spectra in ${B_{1g}}$ symmetry. For clarity data other than those measured at the lowest temperature are shifted vertically. The peak positions are marked by the black vertical lines. Two dashed lines indicate the peak positions at 7.2 K. The increase of the 13 K spectrum toward higher energies originates from a local surface contamination. (b) Normalized temperature dependence of peak positions. They are extracted by a fit with a Lorentzian function. For normalization of the energy the $\Omega_{0}$ = 110 and 190 cm-1 for peak 1 and peak 2, respectively. Table 1: Comparison between the peak positions in the ${B_{1g}}$ Raman spectra and the gap energies obtained by STS and ARPES.  | Raman (cm-1) | Raman (meV) |  | STS (meV) Du  _et al._ (2016); Chen  _et al._ (2019) | ARPES (meV) Zhao  _et al._ (2016)
---|---|---|---|---|---
$\Omega_{\rm peak}^{(1)}$ | 110$\pm$0.5 | 13.75$\pm$0.06 | $2\Delta_{1}$ | 17.2$\pm$2.0 | –
$\Omega_{\rm peak}^{(2)}$ | 190$\pm$0.5 | 23.75$\pm$0.06 | $2\Delta_{2}$ | 28.4$\pm$4.0 | 26.0$\pm$4.0
$\Omega_{\rm peak}^{(1)}/k_{B}T_{c}$ | 3.80$\pm$0.02 | $2\Delta_{1}/k_{B}T_{c}$ | 4.75$\pm$0.55 | –
$\Omega_{\rm peak}^{(2)}/k_{B}T_{c}$ | 6.56$\pm$0.02 | $2\Delta_{2}/k_{B}T_{c}$ | 7.84$\pm$1.10 | 7.18$\pm$1.10

In Fig. 3, we show the variation with temperature of the ${B_{1g}}$ difference
spectra. With increasing temperature, the two peaks shift to lower energy and
cannot be resolved any further above 38 K. Both peaks depend more weakly on
temperature than expected from the BCS theory. As opposed to the results in
Ba1-xKxFe2As2 Böhm  _et al._ (2014) this temperature dependence does not allow
us to clearly identify the origin of the superconducting structures. Whereas
the pair-breaking peaks do not necessarily follow the BCS prediction, except
in weakly interacting systems Devereaux and Einzel (1995), at least BS
collective modes are expected to follow the related single-particle gap
$\Delta(T)$ Monien and Zawadowski (1990). The temperature dependence of
Leggett modes has not been analyzed yet, but is presumably more complicated
since the coupling of at least two gaps has to be considered Suhl  _et al._
(1959); Khodas and Chubukov (2012). Thus, the variation with temperature is
not an identification criterion.

Figure 4: Interactions and related phenomenological spectra. (a),(c)
Illustration of the superconducting ground state and possible subdominant
interactions in the first BZ. The blue solid lines indicate the first quadrant
of the 1 Fe BZ boundary in (a) and (c). The 2 Fe BZ boundary in (c) is
illustrated by the red dashed line. The Fermi pockets are shown by the black
solid curves. The red double arrows indicate the pairing interaction in (a)
and (c). The green double arrow in (a) and the orange and dark blue ones in
(c) indicate the possible subdominant pairing interactions. (b), (d) The
fitting results correspond to the situations in (a) and (c). The red curves
are the fitting results. The gray solid dots are the experimental data. The BS
modes or Leggett mode and pair-breaking (PB) peaks are indicated in (b) and
(d).

What are the possible explanations and the implications thereof for
superconductivity in (Li1-xFex)OHFeSe? Given the STS and ARPES results a
scenario with two isotropic $s$-wave gaps appears to be natural and compatible
with the Raman spectra. However, the gap energies observed by Raman scattering
are significantly too small (Table 1), and the phenomenology returns poor
agreement with the data (see Supplemental Material G) . Rather, the shape of
the line at 110 cm-1 is strongly indicative of a collective mode. An undamped
quadrupolar excitation inside the gap Thorsmølle  _et al._ (2016) or a nematic
resonance Gallais  _et al._ (2016); Kang and Fernandes (2016) are not very
likely since the related fluctuations above $T_{c}$ could not be observed (see
Supplemental Material F). The distinction between a BS and a Leggett mode is
less obvious. Since the mode is apparently below the edge of the smaller gap,
the damping is small in either case.

A BS mode indicates orthogonal pairing channels Bardasis and Schrieffer
(1961). Given the possible interactions between the two electron bands at the
$X$ and $Y$ points, as shown in Fig. 4 (a), only a $d_{x^{2}-y^{2}}$ gap of
lowest (first) order is possible if spin fluctuations are relevant. In this
case the gaps on the two bands have opposite sign and, since there are no
nodes on the Fermi surface, the gap is clean in agreement with the STS and
Raman results. Then the BS mode would result from a $d_{x^{2}-y^{2}}$
interaction of second or higher order [for second-order ${B_{1g}}$ symmetry,
cf. Fig. 1(c5)], having a smaller coupling strength than the ground state. The
description of the asymmetric maximum at 110 cm-1 is reasonable, as shown in
Fig. 4 (b), and yields a coupling strength of approximately 0.16 for the
subleading channel (see Supplemental Material G).

This interpretation favors pair breaking as a possible explanation for the
high-energy peak at 190 cm-1. However, the symmetric shape is untypical for a
pair-breaking feature and may only be explained by a broader gap distribution
as suggested by the STS data for instance (see Table 1). Yet, the energy of
190 cm-1 (23.6 meV) is significantly below the single-particle gap. More
importantly, it is difficult to explain why there should be two rather
distinct gaps for this scenario of equivalent bands. It is, in fact, more
likely that there is a $(\pi,\pi)$ reconstruction of the Fermi surface and a
hybridization between the electron bands, as suggested by Khodas et al. Khodas
and Chubukov (2012); Khodas  _et al._ (2014). Then one expects two concentric
Fermi surfaces with distinctly different gaps as seen here in the Raman data
and also in the STS dataDu  _et al._ (2016); Chen  _et al._ (2019). ARPES Zhao
_et al._ (2016) and STS Du  _et al._ (2016) tell us that the outer Fermi
surface has the larger gap. Superconductivity would then arise from the
comparably strong interaction between these hybridized bands, inducing a
repulsion of the gap energies Du  _et al._ (2018) and either a collective
Leggett mode in the ${A_{1g}}$ channel Cea and Benfatto (2016) or a double-
peak structure well below the gap for certain parameter ranges in the case of
a sign change of the gap between the bands ($s_{\pm}$ gap) Khodas  _et al._
(2014).

Since only the ${B_{1g}}$ channel displays a distinct redistribution of
spectral weight below $T_{c}$, scenarios which include other channels are less
likely to explain the results. Thus, in addition to the strongly coupled
ground state resulting from the interaction between the concentric bands there
must be a weaker $(\pi,\pi)$ interaction leading to collective modes in
${B_{1g}}$ symmetry, as shown in Fig. 4(c). In this scenario, the sharp mode
is either a Leggett mode from a weak $(\pi,\pi)$ coupling between $X$ and $Y$
on top of the strongly coupled ground state resulting from the strong coupling
of the concentric bands or a BS mode having a similar origin. Huang et al.
Huang  _et al._ (2018) indeed argue that the distinction between Leggett and
BS modes becomes obsolete here.

The only remaining issue concerns the positions of the Raman peaks that appear
at smaller energies than in the single-particle spectroscopies. Whereas the
energy of the Raman maximum at 110 cm-1 is naturally explained in terms of a
BS mode appearing below the related single-particle gap at
136$\,\mathrm{cm}^{-1}$ (17 meV) and manifesting itself as a shoulder in the
Raman spectrum, the position of the mode at 190 cm-1 is less obvious since
there is no additional pair-breaking feature in the spectra. In principle it
could be another BS mode pulled down by 10% from the gap edge at approximately
220$\,\mathrm{cm}^{-1}$ (27 meV) and suppressing the pair breaking almost
entirely. Also in this case the ground state would be induced by a strong
hybridization of the two electron bands. The resulting description of the
experimental data is in fact much better in this case [see Fig. 4 (d)].

Finally, since the hole band is rather close to the Fermi surface, spin-
fluctuation pairing between the incipient hole band and the electron bands can
still be rather strong Linscheid  _et al._ (2016); Mishra  _et al._ (2016).
Then the ground state is $s_{\pm}$ with all the electron bands having the same
sign, and the ${B_{1g}}$ modes are collective $d$-wave modes from the
hybridized electron bands. Whether or not the magnitudes of the gaps observed
on the electron bands and the related subleading pairing strengths are
compatible with these considerations needs to be worked out theoretically. The
resulting spectra could be very similar to those in the previous case.
However, due to the pairing-induced renormalization of the central hole band
below $T_{c}$ one would not expect the ${A_{1g}}$ spectra to be entirely
insensitive to superconductivity. Thus the scenario of an incipient band is
less supported by the present experiment.

In conclusion, we studied the polarization- and temperature-dependent Raman
spectra in (Li1-xFex)OHFeSe. Superconductivity affects only the ${B_{1g}}$
spectra. One of the observed modes is resolution limited arguing strongly for
its collective character. For the surprisingly successful description of the
data in terms of BS modes there are essentially two scenarios: (i) dominant
pairing between the hybridized electron bands Khodas and Chubukov (2012) and
subleading ($\pi$, $\pi$) interactions between the electron bands. Then
Leggett and BS modes cannot be distinguished Huang  _et al._ (2018). (ii) If
the ground state originates from the interaction between the incipient hole
band and the electron bands a similar collective mode may be expected, but
completely inert ${A_{1g}}$ spectra are unlikely in this case making scenario
(i) more likely. Yet, the distinction between the two scenarios requires
quantitative theoretical studies.

###### Acknowledgements.

We thank L. Benfatto, P. Hirschfeld, T. Maier and L. Zhao for fruitful
discussions. This work is supported by the Deutsche Forschungsgemeinschaft
(DFG) through the coordinated programme TRR80 (Projekt-ID 107745057) and
project HA2071/12-1. G. H. would like to thank the Alexander von Humboldt
Foundation for support from a research fellowship. The work at China was
supported by National Natural Science Foundation of China (No. 11834016), and
the National Key Research and Development Program of China (Grant No.
2017YFA0303003) and Key Research Program of Frontier Sciences of the Chinese
Academy of Sciences (Grant No. QYZDY-SSW-SLH001). ).

## Supplementary information

A. Sample characterization

X-ray diffraction (XRD) measurements were carried out on a 9 kW Rigaku
SmartLab x-ray diffractometer. Magnetization measurements were done using a
Quantum Design MPMS-XL1 system with a measuring field of 1 Oe to characterize
the superconducting state. Figure 5 (a) and (b) show the XRD characterizations
of the (Li1-xFex)OHFeSe films. The XRD pattern exhibits a single preferred
orientation of (00$l$) for (Li1-xFex)OHFeSe, and the peaks of LaAlO3
substrates are marked with LAO. The corresponding full width at half maximum
(FWHM) of the rocking curve for the (006) reflection is 0.22∘, indicative of
the high crystalline quality. Figure 5 (c) presents the typical temperature
dependent magnetization of a small piece of a (Li1-xFex)OHFeSe film. Sharp
transition was reproduced for several small pieces (approximately 2$\times$2
mm2) of thin films of the same batch and indicates that the sample is
homogeneous. The $T_{c}$, obtained from the onset temperature of diamagnetism,
of optimal (Li1-xFex)OHFeSe films is 42 K. The transition width (10 to 90 % of
the full signal) is less than 2 K. The film thickness is on the order of 100
nm as reported in Ref. Huang  _et al._ (2017a).

Figure 5:  Characterizations of the (Li1-xFex)OHFeSe film on LaAlO3 substrate.
(a) $\theta-2\theta$ scan. The scan shows only (00$l$) peaks. The LAO peaks
are indicated. (b) Rocking curve of (006) reflection. The FWHM is as small as
0.22∘. (c) Temperature dependence of field-cooling (FC) magnetization with the
magnetic field of 1 Oe. The diamagnetic transition occurs at $\sim$42$\pm 1$
K.

B. Determination of the spot temperature

The real temperature in the illuminated spot is higher than the holder
temperature due to the laser heating effect. In order to figure out the
heating, we measure the spectra in ${B_{1g}}$ symmetry at 30 K and 33 K with
different laser power. As seen in Fig. 6, the peak positions measured at 30 K
with an absorbed laser power of 4 mW and at 33 K measured with 2 mW are almost
the same. Therefore, we can estimate that the laser induced heating is
approximately 1.5 K/mW.

Figure 6: Determination of the spot temperature. The spectra were obtained in
${B_{1g}}$ symmetry.

C. Surface quality

To clarify the origin of the additional peaks appearing in the range of
220-260 cm-1 in ${A_{1g}}$, ${B_{1g}}$ and ${B_{2g}}$ symmetry, we compared
the spectra in cleaved and as-grown regions as shown in Fig. 7. It can be
clearly seen that the additional peaks disappear in the cleaved region in both
${A_{1g}}$ and ${B_{2g}}$ symmetry. It means that they stem from surface
contamination or degradation, and originate probably from the ${A_{1g}}$,
${E_{g}}$ (1), ${E_{g}}$ (2) and ${E_{g}}$ (3) phonons of Fe oxides Shim and
Duffy (2002). We can rule out a luminescence effect since these peaks stay
pinned and only change intensity for different laser lines (see SI.D) as shown
in Fig. 9. Because the film is too thin, the ${A_{1g}}$ and ${E_{g}}$ phonons
of the substrate (LaAlO3) Sathe and Dubey (2007) become stronger after
cleavage.

Figure 7:  Raman spectra in cleaved and as-grown regions of the
(Li1-xFex)OHFeSe sample. The ${A_{1g}}$ and ${E_{g}}$ phonons of the substrate
(LaAlO3) can be observed in the cleaved region. The ${A_{1g}}$ and ${B_{1g}}$
phonons of (Li1-xFex)OHFeSe sample are identified in both cleaved and as-grown
region. The blue arrows indicate the additional peaks induced by surface
contamination.

D. Resonances

Figure 8 shows the Raman spectra in ${B_{1g}}$ symmetry for the excitation
lines at 457 and 577 nm. The peak at 250 cm-1 is enhanced for blue laser
excitation. The peak at 50 cm-1 is a plasma line in the blue. The peak at 35
cm-1 with yellow excitation, is of unknown origin. However, neither of them is
related to superconductivity since they also appear above the superconducting
transition temperature. The superconductivity-induced peaks are located at the
same positions for both yellow and blue laser (see dashed lines in Fig. 9),
which means that resonance effects play a minor or no role.

Figure 8:  Comparison between blue laser excitation and yellow laser
excitation in ${B_{1g}}$ symmetry at 7.2 K and 2 mW.  Figure 9:
Superconducting and normal state Raman spectra of (Li1-xFex)OHFeSe thin film
in ${B_{1g}}$ channel taken from the second spot. (a) Raman spectra at 7.2 K
and 46 K. (b) The difference spectra between 7.2 K and 46 K, i.e.
$R\chi^{{}^{\prime\prime}}(\Omega,7.2K)-R\chi^{{}^{\prime\prime}}(\Omega,46K)$.
Figure 10: Raman spectra in ${B_{1g}}$ symmetry at various temperatures.
Partially due to the limited rejection of stray light fluctuations cannot be
detected in (Li1-xFex)OHFeSe.

E. Data reproducibility

Figure 9 shows the Raman spectra in ${B_{1g}}$ symmetry taken from the second
spot on the same sample. The two superconductivity-induced peaks can be
observed at 110 and 190 cm-1, which are consistent with the results of the
first spot (Fig. 2 in the main text). The peaks at 165 cm-1 and 240 cm-1 are
ascribed to a leakage of the ${A_{1g}}$ phonon and phonons of Fe oxides
respectively, as described in the main text and SI. C. All these features
indicate that our data are well reproducible.

F. Fluctuations

Gallais et al. proposed that nematic fluctuations may enhance the pair-
breaking peak resonantly Gallais  _et al._ (2016). In fact, if there are
fluctuations in the system, an extra contribution may be superimposed on the
electronic continuum in the Raman spectra Kretzschmar  _et al._ (2016). This
contribution cannot be observed in the ${B_{1g}}$ Raman spectra of
(Li1-xFex)OHFeSe above $T_{c}$ as shown in Fig. 10. One could argue that the
low-energy part is not as well resolved in (Li1-xFex)OHFeSe as in
Ba(Fe1-xCox)2As2 as a consequence of the poor stray-light rejection. Yet, it
is fair to say that there is no clear evidence for temperature dependence
comparable to that in Ba(Fe1-xCox)2As2 (see, e.g. our recent paper in Phil.
Mag. Lederer  _et al._ (2020)). In addition, the spectra in (Li1-xFex)OHFeSe
and optimally doped Ba(Fe1-xCox)2As2 look very different, and, although more
work on (Li1-xFex)OHFeSe samples having better surfaces is necessary, we tend
to conclude that fluctuations cannot be observed in (Li1-xFex)OHFeSe above
$T_{c}$. We note that the line shape of a nematic resonance cannot a priory be
distinguished from a BS mode on the basis of just one doping level.

Table 2: Fitting parameters of difference spectra in ${B_{1g}}$ symmetry. Parameters | 2$\Delta_{0}^{1}$(cm-1) | 2$\Delta_{0}^{2}$(cm-1) | $\lambda_{d}^{1}$ | $\lambda_{d}^{2}$ | $\lambda_{s}$ | $\delta$ | $b_{1}$ | $b_{2}$ |  $\sigma$(cm-1)
---|---|---|---|---|---|---|---|---|---
BS or Leggett mode + PB peak | 112 | 186 | 0.16 | 0 | 1 | 0.1 | 0.53 | 0.30 | 3
Two BS modes | 112 | 224 | 0.16 | 0.40 | 1 | 0.1 | 0.53 | 0.12 | 3
Two PB peaks | 108 | 186 | 0 | 0 | 1 | 0.1 | 0.82 | 0.20 | 3
Figure 11:  Plots of the Raman response for different coupling strength
$\lambda_{d}$. Here the parameters are set as follows, $\delta$ = 0.1,
$\lambda_{s}$ = 1, $\sigma$ = 3 cm-1. Figure 12: Difference Raman spectrum
between 7.2 K and 46 K in ${B_{1g}}$ symmetry (grey points) and fitting with
two istropic s-wave gaps (red curve).

G. Phenomenology of collective modes

For describing the data we adopted the phenomenology proposed by Scalapino and
Devereaux Scalapino and Devereaux (2009). In this theory, the $s$-wave state
pairing strength exceeds that of the $d$-wave state. Therefore, except for the
$s$-wave ground state, there could be a $d$-wave (${B_{1g}}$) collective
excitonic mode. The Raman scattering is given by

$\displaystyle
R\chi^{{}^{\prime\prime}}(\omega)=\frac{4\pi}{\omega}\left\langle\frac{\Delta(\theta)^{2}}{\sqrt{\omega^{2}-(2\Delta(\theta))^{2}}}\right\rangle+Im\left\\{\frac{2\left\langle\frac{2\Delta(\theta)}{\omega}\overline{P}(\omega,\theta)\right\rangle^{2}}{\left(\frac{1}{\lambda_{d}}-\frac{1}{\widetilde{\lambda_{s}}}\right)-\langle\overline{P}(\omega,\theta)\rangle-i\delta}\right\\}.$
(1)

Here, $\Delta(\theta)$ is the $s$-wave ground-state gap. $\lambda_{s}$ and
$\lambda_{d}$ are the $s$-wave and $d$-wave pairing strength respectively. A
small damping term $\delta$ has been added to control the integrals. The
function $\overline{P}(\omega,\theta)$ is defined by

$\displaystyle\overline{P}(\omega,\theta)=\left\\{\begin{aligned}
&\frac{\omega/2\Delta(\theta)}{\sqrt{1-[\omega/2\Delta(\theta)]^{2}}}\sin^{-1}\left[\frac{\omega}{2\Delta(\theta)}\right],\left|\frac{\omega}{2\Delta(\theta)}\right|<1\\\
&\frac{\omega/2\Delta(\theta)}{\sqrt{[\omega/2\Delta(\theta)]^{2}-1}}\left[\ln\left\\{\left|\frac{\omega}{2\Delta(\theta)}\right|-\sqrt{[\omega/2\Delta(\theta)]^{2}-1}\right\\}+i\frac{\pi}{2}\right],\left|\frac{\omega}{2\Delta(\theta)}\right|>1\\\
\end{aligned}\right.$ (2)

Since the line at 110 cm-1 is very narrow the resolution of the spectrometer
cannot be neglected and a convolution with a Gaussian is necessary for
properly describing the Raman scattering intensity

$I(\omega_{0})=\frac{R}{\sqrt{2\pi}\sigma}\int\chi^{\prime\prime}{(\omega)}e^{\frac{(\omega-\omega_{0})^{2}}{2\sigma^{2}}}d\omega,$
(3)

where $\sigma^{2}$ is the variance of the Gaussian.

Figure 11 shows plots of the Raman response for different coupling strength
$\lambda_{d}$. With the increase of the $\lambda_{d}$, the pair breaking peak
gets suppressed. In our case, there are two distinct superconducting energy
gaps. To simplify, we treat them as two independent isotropic gaps
$\Delta_{i}$ on the different bands as suggested by the tunneling experiments,
that contribute to Raman scattering separately. The total Raman intensity is
expressed as

$I_{\rm total}(\omega)=\sum_{i=1,2}b_{i}I_{i}(\omega),$ (4)

where $b_{i}$ is a dimensionless constant.

Based on this phenomenological model, we fit the difference Raman spectrum in
${B_{1g}}$ symmetry by considering three situations, i.e. two pair-breaking
peaks (see Fig. 11), one BS or Leggett mode plus one pairing-breaking mode
(see Fig.4(b) in the main text) and two BS modes (see Fig.4(d) in the main
text). The fitting parameters are listed in Table 2. It is obvious that the
two peaks in ${B_{1g}}$ spectrum can not be explained by only considering two
isotropic s-wave gaps.

## References

  * Paglione and Greene (2010) J. Paglione and R. L. Greene, Nature Phys. 6, 645 (2010).
  * Mazin  _et al._ (2008) I. I. Mazin, D. J. Singh, M. D. Johannes, and M. H. Du, Phys. Rev. Lett. 101, 057003 (2008).
  * Scalapino (2012) D. J. Scalapino, Rev. Mod. Phys. 84, 1383 (2012).
  * Onari and Kontani (2009) S. Onari and H. Kontani, Phys. Rev. Lett. 103, 177001 (2009).
  * Kontani and Onari (2010) H. Kontani and S. Onari, Phys. Rev. Lett. 104, 157001 (2010).
  * Borisenko  _et al._ (2016) S. V. Borisenko, D. V. Evtushinsky, Z. H. Liu, I. Morozov, R. Kappenberger, S. Wurmehl, B. Buchner, A. N. Yaresko, T. K. Kim, M. Hoesch, T. Wolf, and N. D. Zhigadlo, Nature Phys. 12, 311 (2016).
  * Lederer  _et al._ (2015) S. Lederer, Y. Schattner, E. Berg, and S. A. Kivelson, Phys. Rev. Lett. 114, 097001 (2015).
  * Lu  _et al._ (2015) X. F. Lu, N. Z. Wang, H. Wu, Y. P. Wu, D. Zhao, X. Z. Zeng, X. G. Luo, T. Wu, W. Bao, G. H. Zhang, F. Q. Huang, Q. Z. Huang, and X. H. Chen, Nature Mater. 14, 325 (2015).
  * Burrard-Lucas  _et al._ (2013) M. Burrard-Lucas, D. G. Free, S. J. Sedlmaier, J. D. Wright, S. J. Cassidy, Y. Hara, A. J. Corkett, T. Lancaster, P. J. Baker, S. J. Blundell, and S. J. Clarke, Nature Mater. 12, 15 (2013).
  * Zhao  _et al._ (2016) L. Zhao, A. J. Liang, D. N. Yuan, Y. Hu, D. F. Liu, J. W. Huang, S. L. He, B. Shen, Y. Xu, X. Liu, L. Yu, G. D. Liu, H. X. Zhou, Y. L. Huang, X. L. Dong, F. Zhou, K. Liu, Z. Y. Lu, Z. X. Zhao, C. T. Chen, Z. Y. Xu, and X. J. Zhou, Nature Communi. 7, 10608 (2016).
  * Nekrasov and Sadovskii (2015) I. A. Nekrasov and M. V. Sadovskii, JETP Lett. 101, 47 (2015).
  * Shi  _et al._ (2017) X. Shi, Z.-Q. Han, X.-L. Peng, P. Richard, T. Qian, X.-X. Wu, M.-W. Qiu, S. C. Wang, J. P. Hu, Y.-J. Sun, and H. Ding, Nature Commun. 8, 14988 (2017).
  * Du  _et al._ (2016) Z. Y. Du, X. Yang, H. Lin, D. L. Fang, G. Du, J. Xing, H. Yang, X. Y. Zhu, and H. H. Wen, Nature Communi. 7, 10565 (2016).
  * Du  _et al._ (2018) Z. Y. Du, X. Yang, D. Altenfeld, Q. Q. Gu, H. Yang, I. Eremin, P. J. Hirschfeld, I. I. Mazin, H. Lin, X. Y. Zhu, and H. H. Wen, Nature Phys. 14, 134 (2018).
  * Chen  _et al._ (2019) C. Chen, Q. Liu, T. Z. Zhang, D. Li, P. P. Shen, X. L. Dong, Z. X. Zhao, T. Zhang, and D. L. Feng, Chin. Rev. Lett. 36, 057403 (2019).
  * Devereaux and Hackl (2007) T. P. Devereaux and R. Hackl, Rev. Mod. Phys. 79, 175 (2007).
  * Scalapino and Devereaux (2009) D. J. Scalapino and T. P. Devereaux, Phys. Rev. B 80, 140512(R) (2009).
  * Maiti  _et al._ (2016) S. Maiti, T. A. Maier, T. Böhm, R. Hackl, and P. J. Hirschfeld, Phys. Rev. Lett. 117, 257001 (2016).
  * Kretzschmar  _et al._ (2013) F. Kretzschmar, B. Muschler, T. Böhm, A. Baum, R. Hackl, H. H. Wen, V. Tsurkan, J. Deisenhofer, and A. Loidl, Phys. Rev. Lett. 110, 187002 (2013).
  * Böhm  _et al._ (2014) T. Böhm, A. F. Kemper, B. Moritz, F. Kretzschmar, B. Muschler, H. M. Eiter, R. Hackl, T. P. Devereaux, D. J. Scalapino, and H. H. Wen, Phys. Rev. X 4, 041046 (2014).
  * Abrikosov and Fal’kovskii (1961) A. A. Abrikosov and L. A. Fal’kovskii, Zh. Eksp. Teor. Fiz. 40, 262 (1961), [Sov. Phys. JETP 13, 179 (1961)].
  * Klein and Dierker (1984) M. V. Klein and S. B. Dierker, Phys. Rev. B 29, 4976 (1984).
  * Bardasis and Schrieffer (1961) A. Bardasis and J. R. Schrieffer, Phys. Rev. 121, 1050 (1961).
  * Leggett (1966) A. J. Leggett, Prog.Theor. Phys. 36, 901 (1966).
  * Klein (2010) M. V. Klein, Phys. Rev. B 82, 014507 (2010).
  * Blumberg  _et al._ (2007) G. Blumberg, A. Mialitsin, B. S. Dennis, M. V. Klein, N. D. Zhigadlo, and J. Karpinski, Phys. Rev. Lett. 99, 227002 (2007).
  * Huang  _et al._ (2017a) Y. L. Huang, Z. P. Feng, S. L. Ni, J. Li, W. Hu, S. B. Liu, Y. Y. Mao, H. X. Zhou, F. Zhou, K. Jin, H. B. Wang, J. Yuan, X. L. Dong, and Z. X. Zhao, Chin. Rev. Lett. 34, 077404 (2017a).
  * Huang  _et al._ (2017b) Y. L. Huang, Z. P. Feng, J. Yuan, W. Hu, J. Li, S. L. Ni, S. B. Liu, Y. Y. Mao, H. X. Zhou, H. B. Wang, F. Zhou, G. M. Zhang, K. Jin, X. L. Dong, and Z. X. Zhao, arXiv: 1711.02920 (2017b).
  * Zhang  _et al._ (2019) A. Zhang, X. Ma, Y. Wang, S. Sun, B. Lei, H. Lei, X. Chen, X. Wang, C. Chen, and Q. Zhang, Phys. Rev. B 100, 060504(R) (2019).
  * Muschler  _et al._ (2009) B. Muschler, W. Prestel, R. Hackl, T. P. Devereaux, J. G. Analytis, J. H. Chu, and I. R. Fisher, Phys. Rev. B 80, 180510(R) (2009).
  * Jost  _et al._ (2018) D. Jost, J.-R. Scholz, U. Zweck, W. R. Meier, A. E. Böhmer, P. C. Canfield, N. Lazarević, and R. Hackl, Phys. Rev. B 98, 020504(R) (2018).
  * Böhm  _et al._ (2018) T. Böhm, F. Kretzschmar, A. Baum, M. Rehm, D. Jost, R. H. Ahangharnejhad, R. Thomale, C. Platt, T. A. Maier, W. Hanke, B. Moritz, T. P. Devereaux, D. J. Scalapino, S. Maiti, P. J. Hirschfeld, P. Adelmann, T. Wolf, H. H. Wen, and R. Hackl, Npj Quantum Mater. 3, 48 (2018).
  * Devereaux and Einzel (1995) T. P. Devereaux and D. Einzel, Phys. Rev. B 51, 16336 (1995).
  * Monien and Zawadowski (1990) H. Monien and A. Zawadowski, Phys. Rev. B 41, 8798 (1990).
  * Suhl  _et al._ (1959) H. Suhl, B. T. Matthias, and L. R. Walker, Phys. Rev. Lett. 3, 552 (1959).
  * Khodas and Chubukov (2012) M. Khodas and A. V. Chubukov, Phys. Rev. Lett. 108, 247003 (2012).
  * Thorsmølle  _et al._ (2016) V. K. Thorsmølle, M. Khodas, Z. P. Yin, C. L. Zhang, S. V. Carr, P. C. Dai, and G. Blumberg, Phys. Rev. B 93, 054515 (2016).
  * Gallais  _et al._ (2016) Y. Gallais, I. Paul, L. Chauvière, and J. Schmalian, Phys. Rev. Lett. 116, 017001 (2016).
  * Kang and Fernandes (2016) J. Kang and R. M. Fernandes, Phys. Rev. Lett. 117, 217003 (2016).
  * Khodas  _et al._ (2014) M. Khodas, A. V. Chubukov, and G. Blumberg, Phys. Rev. B 89, 245134 (2014).
  * Cea and Benfatto (2016) T. Cea and L. Benfatto, Phys. Rev. B 94, 064512 (2016).
  * Huang  _et al._ (2018) W. Huang, M. Sigrist, and Z. Y. Weng, Phys. Rev. B 97, 144507 (2018).
  * Linscheid  _et al._ (2016) A. Linscheid, S. Maiti, Y. Wang, S. Johnston, and P. J. Hirschfeld, Phys. Rev. Lett. 117, 077003 (2016).
  * Mishra  _et al._ (2016) V. Mishra, D. J. Scalapino, and T. A. Maier, Sci. Rep. 6, 32078 (2016).
  * Shim and Duffy (2002) S. H. Shim and T. S. Duffy, American Mineralogist 87, 318 (2002).
  * Sathe and Dubey (2007) V. G. Sathe and A. Dubey, J. Phys.: Condens. Matter 19, 382201 (2007).
  * Kretzschmar  _et al._ (2016) F. Kretzschmar, T. Böhm, U. Karahasanović, B. Muschler, A. Baum, D. Jost, J. Schmalian, S. Caprara, M. Grilli, C. Di Castro, J. G. Analytis, J. H. Chu, I. R. Fisher, and R. Hackl, Nat. Phys. 12, 560 (2016).
  * Lederer  _et al._ (2020) S. Lederer, D. Jost, T. Böhm, R. Hackl, E. Berg, and S. A. Kivelson, Phil. Mag. 100, 2477 (2020).


# Study of Mach reflection in inviscid flows

###### Abstract

In this paper, we study the Mach reflection phenomenon in inviscid flows using
a higher order discontinuous Galerkin method and overset grids. We use the
shock capturing procedure proposed in [1] using overset grids to capture the
discontinuities occurring in the supersonic flow over a wedge accurately. In
this procedure, we obtain a coarse grid solution first and using the troubled
cell data, we construct an overset grid which is approximately aligned to all
the discontinuities. We rerun the solver with the coarse grid solution as the
initial condition while using the troubled cell indicator and the limiter only
on the overset grid. This allows us to capture the discontinuities accurately.
Using this procedure, we have obtained the solution for Mach $3.0$ and $4.0$
flow over a wedge for various wedge angles and determined the detachment
criterion and the Von Neumann condition accurately. We have also determined
the Mach stem height for various wedge angles for these Mach numbers. We have
also demonstrated the hysteresis that occurs in the transition from regular
reflection to Mach reflection.

Keywords: Mach reflection, discontinuous Galerkin method, overset grids,
hysteresis

##  1 Introduction

In this paper, we study the transition between regular reflection (RR) and
Mach reflection (MR) of steady shock waves in inviscid flows using
discontinuous Galerkin method (DGM) along with overset grids. This has been
studied quite extensively in literature experimentally [2], numerically [3],
and analytically [4]. We use the shock capturing procedure proposed in [1]
using overset grids and a higher order method (DGM) to capture the
discontinuities occurring in the supersonic flow over a wedge accurately. In
this procedure, we obtain a coarse grid solution first and using the solution
and the troubled cell data, we construct an overset grid which is
approximately aligned to all the discontinuities. We rerun the solver with the
coarse grid solution as the initial condition while using the troubled cell
indicator and the limiter only on the overset grid. This allows us to capture
the discontinuities accurately. We solve the Euler equations in the
computational domain shown in Figure 1 using DGM and overset grids for Mach
3.0 and 4.0 flow over a wedge for different wedge angles to determine the
transition between RR and MR. In this way, we determine the detachment
criterion and the Von Neumann condition accurately to demonstrate the
hysteresis that occurs in such a flow.

$H_{m}$$w$$H$SupersonicOutflowWallInflowSymmetryPlane$y$$x$$\theta_{w}$

Figure 1: Computational domain for Mach reflection showing the shock structure
and the Mach stem height ($H_{m}$)

The paper is organized as follows. We describe the formulation of the
discontinuous Galerkin method used for all our results in Section 2, the
procedure used for shock capturing using overset grids is described in Section
3, the results are described in Section 4 and we conclude the paper in Section
5.

##  2 Description of discontinuous Galerkin method

Consider the Euler equations in conservative form as given by

$\frac{\partial\textbf{Q}}{\partial t}+\frac{\partial\textbf{F(Q)}}{\partial
x}+\frac{\partial\textbf{G(Q)}}{\partial y}=0\quad\text{in the
domain}\quad\Omega$ (1)

where $\textbf{Q}=(\rho,\rho u,\rho v,E)^{T}$,
$\textbf{F(Q)}=u\textbf{Q}+(0,p,0,pu)^{T}$ and
$\textbf{G(Q)}=v\textbf{Q}+(0,0,p,pv)^{T}$ with
$p=(\gamma-1)(E-\frac{1}{2}\rho(u^{2}+v^{2}))$ and $\gamma=1.4$. Here, $\rho$
is the density, $(u,v)$ is the velocity, $E$ is the total energy and $p$ is
the pressure. We approximate the domain $\Omega$ by $K$ non overlapping
elements given by $\Omega_{k}$.

We look at solving (1) using the discontinuous Galerkin method. We approximate
the local solution in an element $\Omega_{k}$, where $k$ is the element
number, as a polynomial of order $N$ which is given by:

$Q_{h}^{k}(r,s)=\sum_{i=0}^{N_{p}-1}Q_{i}^{k}\psi_{i}(r,s)$ (2)

where $N_{p}=(N+1)(N+1)$ and $r$ and $s$ are the local coordinates. Here, the
subscript $i$ represents the particular degree of freedom, $h$ represents the
grid size, and the superscript $k$ is the element number. The polynomial basis
used ($\psi_{i}(r,s)$) is the tensor product orthonormalized Legendre
polynomials of degree $N$. The number of degrees of freedom are given by
$N_{p}=(N+1)(N+1)$. Now, using $\psi_{j}(r,s)$ as the test function, the weak
form of the equation (1) is obtained as

$\sum_{i=0}^{N_{p}-1}\frac{\partial Q_{i}^{k}}{\partial
t}\int_{\Omega_{k}}\psi_{i}\psi_{j}d\Omega+\int_{\partial\Omega_{k}}\hat{F}\psi_{j}ds-\int_{\Omega_{k}}\vec{F}\cdot\nabla\psi_{j}d\Omega=0\quad
j=0,\ldots,N_{p}-1$ (3)

where $\partial\Omega_{k}$ is the boundary of $\Omega_{k}$,
$\vec{F}=(\textbf{F(Q)},\textbf{G(Q)})$ and $\hat{F}=\bar{F^{*}}\cdot\hat{n}$
where $\bar{F^{*}}$ is the monotone numerical flux at the interface which is
calculated using an exact or approximate Riemann solver and $\hat{n}$ is the
unit outward normal. This is termed to be $\mathbf{P}^{N}$ based discontinuous
Galerkin method.

Equation (3) is integrated using an appropriate Gauss Legendre quadrature and
is discretized in time by using the fifth order Runge-Kutta time
discretization given in [5] unless otherwise specified. To control spurious
oscillations which occur near discontinuities, a limiter is used with a
troubled cell indicator. We have used the KXRCF troubled cell indicator [6]
and the compact subcell WENO (CSWENO) limiter proposed in [7] for all our
calculations.

##  3 Overset grids and shock capturing

Overset grids consist of multiple grids which overlap each other as shown in
Figure 2. When using DGM on overset grids, there are two possible approaches
to handle data communication between the grids. One is a face based
communication approach developed in [8], where solutions at an overset
interface are obtained from the donor element, and then the boundary condition
is applied weakly by imposing a numerical flux at the flux interpolation
points. The other is an element based communication approach developed in [9],
where the internal degrees of freedom of cells near the overset interface are
obtained from the donor element. We use the new element based communication
approach developed in [10] for the data communication between the grids.

Figure 2: Two overlapping grids (Grid 1 in black and Grid 2 in Red) Figure 3:
Coarse grid solution for Mach 3.0 flow over a $24^{\circ}$ wedge using
$\mathbf{P}^{1}$ based discontinuous Galerkin method Figure 4: Troubled cells
obtained using the KXRCF troubled cell indicator [6] for Mach 3.0 flow over a
$24^{\circ}$ wedge

For capturing the shocks accurately, we use the procedure developed in [1]
using overset grids. A brief explanation of the procedure is given below:

Step 1: Run the solver on a coarse grid with a given troubled cell indicator
and limiter to steady state and obtain the solution. As an example, we show
the coarse grid solution obtained for Mach 3.0 flow over a $24^{\circ}$ wedge
in Figure 3

Step 2: Look at the troubled cells to locate the discontinuities (shocks) that
occur in the solution. The troubled cell profile obtained for the Mach 3.0
flow over a $24^{\circ}$ wedge using the KXRCF troubled cell indicator [6] is
shown in Figure 4. From this figure, we can see that the troubled cells give
us a good idea of the location of the shocks, the contact discontinuity and
the initial expansion that forms near the wedge.

Step 3: Construct an overset grid conforming to the computational domain which
is refined in a direction perpendicular to the discontinuities such that they
are approximately parallel to a grid line. This overset grid also encompasses
all the troubled cells. An example overset grid constructed in such a fashion
for the Mach 3.0 flow over a $24^{\circ}$ wedge is shown in Figure 5.

Figure 5: Example overset grid for Mach 3.0 flow over a $24^{\circ}$ wedge
constructed such that the discontinuities are approximately parallel to a grid
line

Step 4: Using this overset grid, we rerun the solver with the coarse grid
solution as the initial condition. While running the solver, we use the
troubled cell indicator and the limiter only on the overset grid. We also use
a high resolution numerical flux on the overset grid to capture the shock
accurately. We have used the SLAU2 [11] numerical flux in the overset grid and
the less expensive Lax-Friedrichs flux elsewhere. Using this procedure, we
obtain a more accurate solution with the discontinuities approximately aligned
to a grid line. The final solution obtained in this fashion for Mach 3.0 flow
over a $24^{\circ}$ wedge is shown in Figure 6.

Figure 6: Final solution obtained for Mach 3.0 flow over a $24^{\circ}$ wedge
using an overset grid and $\mathbf{P}^{4}$ based discontinuous Galerkin method
with discontinuities captured accurately

Regarding the computational cost of the scheme proposed for the current
problem, we note that the troubled cell data shown in Figure 4 are obtained
from a coarse grid of 20000 elements using $\mathbf{P}^{1}$ based DGM after
converging to a residue of about $1\times 10^{-6}$ which happens in about 8000
iterations. This coarse grid solution is not a fully converged solution but is
good enough to be used as an initial condition for the overset grid higher
order ($\mathbf{P}^{4}$ based DGM) solution. The time taken per iteration per
degree of freedom using $\mathbf{P}^{1}$ based DGM for this problem is about
$1.084\times 10^{-6}$s. This tells us that the coarse grid troubled cell data
is obtained in about $693.76$s which is a little less than 12 minutes. Most of
the work involved is in constructing the overset grid which is done manually.
After the overset grid is constructed and the solver is run on the new grid
with the coarse grid solution as the initial condition, the residue converges
to $1\times 10^{-16}$ in about 32000 iterations. For $\mathbf{P}^{4}$ based
DGM, the time taken per iteration per degree of freedom is about $1.345\times
10^{-6}$s for this problem. This tells us that the final solution (which is
fifth order accurate) is obtained in about $21520$s which is about six hours.
All these calculations are done on a 3.60 GHz, Intel(R) Core(TM) i7-7700 CPU
with a single thread. The calculations for remaining angles are quite similar.
We also note that we have a parallelised code and we obtain the solution much
faster based on the number of threads used.

##  4 Results

1) Mach number 3.0: We solve the two-dimensional Euler equations given by (1)
using the discontinuous Galerkin method for various wedge angles between
$\theta_{w}=19.5^{\circ}$ and $\theta_{w}=24^{\circ}$ near the transition
criterion (transition from Regular reflection to Mach reflection) for Mach
number 3.0 and $w/H=1.0$ in the computational domain shown in Figure 1. We
consider two cases to demonstrate the hysteresis phenomenon. We solve the
equations using an impulsive start as the initial condition for the first
case, and the converged solution for $\theta_{w}=24^{\circ}$ as the second
case. The first case is a numerical model to obtain the detachment criterion
and the second case demonstrates the Von Neumann condition. Using the
procedure outlined in Section 3, we obtain a very accurate solution for each
of the wedge angles and obtain the transition criterion clearly. Note that the
initial conditions mentioned above are used to obtain the coarse grid solution
first and the coarse grid solution is used as an initial condition to get the
fine grid solution using $\mathbf{P}^{4}$ based DGM. To show that the
transition criterion has been captured accurately, we show the solution
obtained for Mach 3.0 flow over a wedge with wedge angles
$\theta_{w}=21.45^{\circ}$ and $\theta_{w}=21.46^{\circ}$ with the case 1
initial conditions in Figures 7 and 8 respectively.

Figure 7: Final solution obtained for Mach 3.0 flow over a $21.45^{\circ}$
wedge with an overset grid and $\mathbf{P}^{4}$ based discontinuous Galerkin
method with discontinuities captured accurately using an impulsive start
Figure 8: Final solution obtained for Mach 3.0 flow over a $21.46^{\circ}$
wedge with an overset grid and $\mathbf{P}^{4}$ based discontinuous Galerkin
method with discontinuities captured accurately using an impulsive start

We also show the table containing the Mach stem height obtained using this
procedure for various wedge angles for both cases of initial conditions in
Table 1. From Table 1, we can see that the transition criterion occurs between
wedge angles $21.45^{\circ}$ and $21.46^{\circ}$ and the Von Neumann
transition criterion occurs between wedge angles $19.6^{\circ}$ and
$19.7^{\circ}$. From the three-shock theory, the transition criterion is at
$\theta_{w}=21.456^{\circ}$ and the Von Neumann condition occurs at
$\theta_{w}=19.656^{\circ}$. This closely agrees with what we have obtained
and validates our procedure.

Table 1: Non-dimensionalised height of the Mach stem ($H_{m}/H$) as a function of wedge angle ($\theta_{w}$) for $M$=3.0 and $w/H=1.0$ using two sets of initial conditions:- Case 1: Impulsive start, Case 2: Converged solution for $\theta_{w}=24^{\circ}$ $\theta_{w}$ |  $H_{m}/H$ for case 1 |  $H_{m}/H$ for case 2
---|---|---
(impulsive start) |  (converged solution for $\theta_{w}=24^{\circ}$)
$24^{\circ}$ |  0.274 |  -
$23.5^{\circ}$ |  0.233 |  0.233
$23^{\circ}$ |  0.188 |  0.188
$22.5^{\circ}$ |  0.153 |  0.153
$22^{\circ}$ |  0.122 |  0.122
$21.46^{\circ}$ |  0.078 |  0.078
$21.45^{\circ}$ |  RR |  0.066
$21^{\circ}$ |  RR |  0.041
$20.5^{\circ}$ |  RR |  0.028
$20^{\circ}$ |  RR |  0.009
$19.7^{\circ}$ |  RR |  0.002
$19.6^{\circ}$ |  RR |  RR
$19.5^{\circ}$ |  RR |  RR

2) Mach number 4.0: We repeat the same procedure for various wedge angles
between $\theta_{w}=20.5^{\circ}$ and $\theta_{w}=27^{\circ}$ near the
transition criterion (transition from Regular reflection to Mach reflection)
for Mach number 4.0 and $w/H=1.0$ in the computational domain shown in Figure
1. Again, we consider two cases to demonstrate the hysteresis phenomenon. We
solve the equations using an impulsive start as the initial condition for the
first case, and the converged solution for $\theta_{w}=27^{\circ}$ as the
second case. We show the table containing the Mach stem height obtained using
this procedure for various wedge angles for both cases of initial conditions
in Table 2. From Table 2, we can see that the transition criterion occurs
between wedge angles $25.6^{\circ}$ and $25.7^{\circ}$ and the Von Neumann
transition criterion occurs between wedge angles $20.9^{\circ}$ and
$20.8^{\circ}$. From the three-shock theory, the transition criterion for Mach
number $4.0$ is at $\theta_{w}=25.61^{\circ}$ and the Von Neumann condition
occurs at $\theta_{w}=20.86^{\circ}$. This closely agrees with what we have
obtained and validates our procedure.

Table 2: Non-dimensionalised height of the Mach stem ($H_{m}/H$) as a function of wedge angle ($\theta_{w}$) for $M$=4.0 and $w/H=1.0$ using two sets of initial conditions:- Case 1: Impulsive start, Case 2: Converged solution for $\theta_{w}=27^{\circ}$ $\theta_{w}$ |  $H_{m}/H$ for case 1 |  $H_{m}/H$ for case 2
---|---|---
(impulsive start) |  (converged solution for $\theta_{w}=27^{\circ}$)
$27^{\circ}$ |  0.344 |  -
$26.5^{\circ}$ |  0.313 |  0.313
$26^{\circ}$ |  0.285 |  0.285
$25.7^{\circ}$ |  0.245 |  0.245
$25.6^{\circ}$ |  RR |  0.224
$25.5^{\circ}$ |  RR |  0.203
$25^{\circ}$ |  RR |  0.183
$24.5^{\circ}$ |  RR |  0.144
$24^{\circ}$ |  RR |  0.115
$23^{\circ}$ |  RR |  0.093
$22^{\circ}$ |  RR |  0.045
$21^{\circ}$ |  RR |  0.005
$20.9^{\circ}$ |  RR |  0.002
$20.8^{\circ}$ |  RR |  RR
$20.5^{\circ}$ |  RR |  RR

##  5 Conclusion

We have demonstrated the transition between regular reflection (RR) and Mach
reflection (MR) of steady shock waves in inviscid flows using discontinuous
Galerkin method (DGM) along with overset grids using an accurate shock
capturing procedure [1]. We have identified the detachment criterion and the
Von Neumann condition accurately and demonstrated the hysteresis which occurs
in the transition. We have also calculated the Mach stem height for various
wedge angles for Mach numbers 3.0 and 4.0. As future work, this shock
capturing procedure can be used to further study the different flow phenomena
that occur when the inflow Mach number is small ($<2.0$).

## References

  * [1] Siva Prasad Kochi S.R. and Ramakrishna M., “Shock capturing with discontinuous Galerkin Method using Overset grids for two-dimensional Euler equations,” 2020.  Paper under review in J. Comput. Phys. Preprint at https://arxiv.org/abs/2003.01378.
  * [2] Chpoun A., Passerel D., Li H. and Ben-Dor G., “Reconsideration of the state-of-the-art of oblique shock wave reflection in steady flows. Part I: Experimental investigation,” J. Fluid. Mech., vol. 301, pp. 19–35, 1995\.
  * [3] Vuillon J., Zeitoun D. and Ben-Dor G., “Reconsideration of the state-of-the-art of oblique shock wave reflection in steady flows. Part 2: Numerical investigation,” J. Fluid. Mech., vol. 301, pp. 37–50, 1995\.
  * [4] Ben-Dor G., Shock wave reflection phenomena.  New York: Springer, 2007.
  * [5] Butcher J.C., Numerical Methods for Ordinary Differential Equations.  3rd edition: John Wiley and Sons, 2016.
  * [6] Krivodonova L., Xin J., Remacle J.-F., Chevaugeon N. and Flaherty J., “Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws,” Appl. Numer. Math., vol. 48, pp. 323–338, 2004.
  * [7] Siva Prasad Kochi S.R. and Ramakrishna M., “A compact subcell WENO limiting strategy using immediate neighbours for Runge-Kutta discontinuous Galerkin methods,” Int. J. Comput. Math., vol. 98(3), pp. 608–626, 2021.
  * [8] Galbraith M.C., Benek J.A., Orkwis P.D. and Turner M.G., “A Discontinuous Galerkin Chimera scheme,” Computers and Fluids, vol. 98, pp. 27–53, 2014.
  * [9] Nastase C.R., Mavriplis D.J. and Sitaraman J., “An Overset Unstructured Mesh Discontinuous Galerkin Approach for Aerodynamic Problems,” AIAA 2011-195, 2011.
  * [10] Siva Prasad Kochi S.R. and Ramakrishna M., “A Discontinuous Galerkin Overset Scheme Using WENO Reconstruction and Subcells for Two-Dimensional Problems,” J. Sci. Comput., vol. 93:35, 2022.
  * [11] Kitamura K. and Shima E., “Towards shock-stable and accurate hypersonic heating computations: A new pressure flux for AUSM-family schemes,” J. Comput. Phys., vol. 245, pp. 62–83, 2013.