#  Diversity in Fermi/GBM Gamma Ray Bursts: New insights from Machine Learning

Dimple<EMAIL_ADDRESS>Chennai Mathematical Institute, Siruseri,
603103 Tamilnadu, India.  Aryabhatta Research Institute of Observational
Sciences, Manora Peak, Nainital-263002, India.  K. Misra Aryabhatta Research
Institute of Observational Sciences, Manora Peak, Nainital-263002, India.  K.
G. Arun Chennai Mathematical Institute, Siruseri, 603103 Tamilnadu, India.

###### Abstract

Classification of gamma-ray bursts (GRBs) has been a long-standing puzzle in
high-energy astrophysics. Recent observations challenge the traditional short
vs. long viewpoint, where long GRBs are thought to originate from the collapse
of massive stars and short GRBs from compact binary mergers. Machine learning
(ML) algorithms have been instrumental in addressing this problem, revealing
five distinct GRB groups within the Swift/BAT light curve data, two of which
are associated with kilonovae (KNe). We corroborate these five classes by
extending this analysis to the Fermi/GBM data using unsupervised ML
techniques. These five clusters are well separated in fluence-duration plane,
hinting at a potential link between fluence, duration and complexities (or
structures) in the light curves of GRBs. Further, we confirm two distinct
classes of KN-associated GRBs. The presence of GRB 170817A in one of the two
KNe-associated clusters lends evidence to the hypothesis that this class of
GRBs could potentially be produced by binary neutron star (BNS) mergers. The
second KN-associated GRB cluster could potentially originate from NS-BH
mergers. Future multimessenger observations of compact binaries in
gravitational waves (GWs) and electromagnetic waves can be paramount in
understanding these clusters better.

GRBs — Progenitors — Kilonova – Machine Learning — PCA-UMAP

##  1 Introduction

Over five decades after their discovery, the origins of gamma-ray bursts
(GRBs) remain a widely debated topic within the astronomy community. The
bimodal distribution of GRB duration within the BATSE data hinted at two
classes of GRBs (Kouveliotou et al., 1993). Based on the duration, $T_{\rm
90}$ (the duration over which 90% of the energy in prompt emission is
released), GRBs are broadly divided into two classes: short and long GRBs,
with a boundary at $T_{\rm 90}$ of 2 sec (Fishman & Meegan, 1995). Over the
last decades, a broad consensus has emerged that long-duration GRBs are
associated with the core-collapse of massive stars (Bloom et al., 2002; Hjorth
et al., 2003; Woosley & Bloom, 2006; Cano et al., 2017), while short-duration
GRBs likely arise from a distinct phenomenon, attributed to the mergers of
compact objects (Paczynski, 1986; Berger et al., 2013; Tanvir et al., 2013;
Abbott et al., 2017; Troja et al., 2019).

Despite the initial success of the classification based on $T_{\rm 90}$, with
data from various GRB detectors covering different energy bands, the
limitations of this classification are widely acknowledged (Fynbo et al.,
2006; Zhang et al., 2009; Qin & Chen, 2013). An overlap between two
populations is also observed, where it is hard to determine if a burst with
mixed characteristics belongs to the short or long class (Fynbo et al., 2006;
Dimple et al., 2022a; Troja et al., 2022; Becerra et al., 2023). Recent
observations suggest further diversity within these classes. For example, some
short GRBs exhibit properties that deviate from the merger model and vice
versa (Bromberg et al., 2013; Dimple et al., 2022b, 2023b; Petrosian &
Dainotti, 2024). Long GRBs 211211A (Rastinejad et al., 2022; Troja et al.,
2022; Gompertz et al., 2022; Yang et al., 2022) and 230307A (Dichiara et al.,
2023; Sun et al., 2023; Levan et al., 2024) have observational evidence of a
kilonova (KN) association in their optical light curves. Similarly, the short
GRB 200826A was found to be associated with a supernova (SN) hinting at its
origin from a collapsar (Ahumada et al., 2021; Rastinejad et al., 2022). These
cases hint at a more complex picture, suggesting that the progenitors
responsible for GRBs might not be as straightforward as previously believed.

Recent advances in machine learning (ML) and artificial intelligence (AI) are
fostering a more detailed understanding of GRBs (Chattopadhyay & Maitra, 2017;
Acuner & Ryde, 2018; Jespersen et al., 2020a; Salmon et al., 2022; Steinhardt
et al., 2023; Garcia-Cifuentes et al., 2023; Mehta & Iyyani, 2024; Negro et
al., 2024). When applied to the analysis of prompt emission light curves,
these techniques have revealed hidden patterns within the GRB population,
suggesting the existence of subclasses beyond the traditional short–long GRB
classification (Jespersen et al., 2020b; Steinhardt et al., 2023; Garcia-
Cifuentes et al., 2023; Dimple et al., 2023a).

An ML-based analysis of the prompt emission light curves using Swift/BAT data
has revealed five distinct clusters within the GRB population, which might
represent distinct progenitors (Dimple et al., 2023a). Among these clusters,
the GRBs associated with KNe were found to lie in two distinct clusters,
suggesting that they may have been produced by different progenitors. The
clusters could correspond to binary neutron star (BNS) and neutron star–black
hole (NS-BH) mergers as GRB progenitors or subclasses of either of them. These
findings, if confirmed, can have profound implications for our understanding
of GRBs and compact binaries. One way to validate these findings is to wait
for a sizeable number of multimessenger GRBs such as GW170817. Given the
current sensitivities of gravitational wave (GW) detectors, this may take a
few years. Another way forward is to validate these findings using data from
other datasets, which constitutes the theme of this paper.

One such dataset is the Fermi Gamma-ray Burst Monitor (GBM) catalog (von
Kienlin et al., 2020), which compiles data of GRBs detected by the Fermi/GBM.
Fermi/GBM data covers a wide energy spectrum, ranging from a few keVs to
several MeVs, and includes nearly twice as many GRBs compared to the Swift/BAT
catalog (Lien et al., 2016). Further, Fermi/GBM observed the multimessenger
event GW170817 (Goldstein et al., 2017), originating from a BNS merger (Abbott
et al., 2017). Studying the GRBs in this catalog can provide a more
comprehensive understanding of GRB progenitors.

This paper presents a comprehensive study of the GRB population using ML
techniques on Fermi/GBM prompt emission light curves. Previous studies by
Acuner & Ryde (2018) and Mehta & Iyyani (2024) employed ML techniques on the
Fermi/GBM data to classify GRBs using spectral parameters. However, our
approach differs from these works. We utilize light curves across different
energy bands and apply unsupervised ML algorithms to cluster GRBs based on
similarities and dissimilarities in their light curves. The paper is
structured as follows: Section 2 describes the methodology employed for the
data analysis and clustering of GRBs. Section 3 illustrates the outcomes of
our ML analysis and discusses their implications. Lastly, Section 4 summarizes
the study, highlighting the main findings and their prospective influence on
future research.

##  2 Methodology

We compiled the prompt emission data files for the GRBs observed from July
2008 to August 2023 from the Fermi/GBM catalog, inclusive of the time-tagged
event (TTE) and response files for the detector configuration sourced from the
Fermi/GBM
archive111https://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermigbrst.html. This
section outlines the methodology used to analyze the GRB data and create the
clustering through ML algorithms.

Figure 1: This figure visualizes the embeddings obtained using PCA-UMAP on the
prompt emission light curves (in three energy bins) of 3516 GRBs. Left:
Connectivity map where each point represents a GRB and edges represent how it
is connected to other bursts. The bursts are color-coded with the clusters
identified by AutoGMM. Right: The respective density map, the deeper the shade
of blue, the higher the density of GRBs, meaning more GRBs are found in that
region. The density map is consistent with the groups identified by AutoGMM.

###  2.1 Data Reduction Pipeline

We designed a systematic pipeline implemented in Python to extract light
curves in various energy bands and prepare the data for the ML algorithms. The
pipeline uses the GbmDataTools package (Goldstein et al., 2022) to access and
analyze GBM data. This pipeline automates several key steps, as listed below:

  1. a)

Detector selection: The process begins by estimating the angles between each
detector and the direction of the burst. It then identifies the Sodium Iodide
(NaI) detector with the smallest angle to the burst and selects the
corresponding Bismuth Germanium Oxide (BGO) detector. The data from these
detectors are used for further analysis.

  2. b)

Light curve extraction: Light curves are extracted at different energy bands:
8–50 keV and 50–300 keV for the NaI detector and 300–1000 keV for the BGO
detector. The pipeline uses a temporal binning of 16 milliseconds from TTE
data files using the TTE function within GbmDataTools.

  3. c)

Background subtraction: The background signal is fitted and subtracted from
the corresponding light curves. The BackgroundFitter function in GbmDataTools
is used for background fitting. The resulting light curves are used for
subsequent processing.

  4. d)

Data denoising: Before using the light curves for ML analysis, unwanted noise,
which could obscure the true variations in the data, is removed. Wavelet
analysis is employed for denoising. The optimal wavelet and decomposition
level are chosen to minimize the mean squared error between the original and
denoised data. This process effectively reduces noise in the time series data.

  5. e)

Dataset preparation for ML algorithms: The background-subtracted light curves
are normalized with the fluence to mitigate the potential algorithmic
distractions caused by the typical disparity in fluence between long and short
GRBs. Following normalization, a standardization procedure is applied to
ensure uniformity in the length of the light curves. All light curves are
shifted to a common starting point and padded with zeros. Subsequently, the
normalized and standardized light curves are concatenated, forming a cohesive
time-series dataset with a shared axis and consistent lengths across different
GRBs. Finally, a discrete-time Fourier transform (DTFT) is employed on the
concatenated data set to preserve the time delay information between various
light curves.

After standardization, we could recover 3516 GRBs from the Fermi/GBM catalog.
For the ML algorithm, we used the data in three energy bands, two energy bands
of NaI (8–50 keVand 50–300 keV) and one energy band of BGO (300–100 keV) as
NaI is more sensitive in the lower energy band, and BGO is more sensitive in
the high energy bands. For both datasets, we used binning of 16ms.

###  2.2 Clustering using unsupervised learning

We employed Uniform Manifold Approximation and Projection (UMAP; McInnes et
al. 2018) with Principal Component Analysis (PCA) for dimensionality
reduction. PCA-UMAP (discussed in detail in Dimple et al. 2023a) allows the
analysis of high-dimensional data in a lower-dimensional space while
preserving its underlying structure. For that, we first extracted the first
500 principal components using PCA, capturing 99% of the data variance. Then
UMAP is used to reduce the dimensionality of the data further using topology.
UMAP builds a network of connections between data points based on similarity
and dissimilarity. This network reflects the inherent structure of the data,
capturing local and global relationships. Then, it projects the data to lower-
dimensional space while maintaining these relationships. The algorithm relies
on two key parameters: n_neighbors, which limits the size of the local
neighborhood considered to build the topological structure, and min_dist,
which controls how tightly packed nearby points are. By adjusting these
parameters, we obtained a more precise and insightful representation of our
data in lower dimensions. Specifically, we set the value for n_neighbors at 25
and min_dist at 0.01.

We generated the connectivity and density maps to further analyze the
resulting low-dimensional representation. Connectivity maps visualize the
probability of a connection between two data points, essentially depicting the
aforementioned network in a two-dimensional format. On the other hand, density
maps, along with the 2-dimensional embeddings, offer additional insights into
the data structure. These maps help to understand the distribution and
concentration of data points in a low-dimensional space. Finally, we
identified clusters within the two-dimensional embedding using the AutoGMM
module (Athey et al., 2019), which applies a Gaussian Mixture Model (GMM) to
automatically determine the optimal number of clusters present in the data.

Figure 2: GRB locations on fluence–$T_{\rm 90}$ plane: GRBs are color coded
based on distinct clusters identified by AutoGMM. Notably, these clusters
exhibit clear separation in the fluence-duration plane, implying a potential
correlation between light curve features and GRB fluence and duration. Figure
3: The connectivity map color coded as traditional long and short duration
GRBs based on their $T_{\rm 90}$ duration. KN and SN-associated GRBs are shown
with red diamonds and blue squares, respectively, on the embeddings. KN-
associated GRBs are labelled with their names and respective colors as per
their traditional classification. Short and long-duration KN-associated GRBs
are well separated into two distinct groups. SN-associated GRBs are scattered
across the map with an unclear clustering pattern. Notably, several SN and KN-
associated GRBs lie at the top of cluster 2. The inset offers a closer view,
revealing that GRBs 211211A and 230307A are adjacent, suggesting similarity in
their light curves. Additionally, some bright SNe are also situated nearby
labeled with their names in grey color.

##  3 Results and discussion

We obtained the two-dimensional embeddings using the method described in
Section 2. The left panel of Figure 1 shows these embeddings with connectivity
maps color-coded with the clusters identified by AutoGMM. The right panel of
the figure displays the respective density maps, providing additional
information about the density of GRBs in the map.

###  3.1 Distinct groups within GRBs

Our method identifies five distinct groups within the embedding, evident from
the density plots (see Figure 1). This is remarkably similar to the findings
of Dimple et al. (2023a), where five distinct clusters were identified using
the same method on the Swift/BAT data. It is also interesting to note that,
based on spectral parameters that govern the emission mechanism, Acuner & Ryde
(2018) showed the existence of five clusters in the Fermi/GBM data. However,
we do not find an obvious overlap of the clusters we identify here with those
of Acuner & Ryde (2018). Besides the Swift/BAT and Fermi/GBM datasets, the
BATSE data also revealed five clusters in a study by Chattopadhyay & Maitra
(2017). Given the differences in the detectors and their energy bands, among
other things, the match in the number of clusters is very intriguing.

The consistency of five distinct clusters across multiple datasets and
methodologies suggests that the GRB population exhibits diversity beyond just
two classes. To further investigate these clusters, we examined their prompt
emission properties. We observed that these clusters exhibit minimal
differences in their spectral properties, such as $E_{\rm p}$ and $\alpha$
(See Figure A.1 in Appendix A). To explore any correlation between fluence and
duration for these clusters, we plotted the GRBs from the five clusters on the
fluence–$T_{\rm 90}$ plane (Figure 2). Interestingly, despite using fluence-
normalized light curves for classification, they are well separated in this
plane. This suggests a strong correlation between the brightness and
complexity of the GRB light curves, a prediction that aligns with the findings
of Stern et al. (1997). This also allows the identification of various
classes: short-duration bursts with low fluence (cluster 0), intermediate-
duration bursts with low fluence (clusters 1 and 3), long-duration bursts with
low fluence (cluster 4), and intermediate to long-duration bursts with high
fluence (cluster 2). Though the correlation between fluence and duration is
well known (Horváth et al., 2001), our results hint at a correlation between
light curve shapes and brightness, which may have a deeper connection with the
physics of the central engines of GRBs. A detailed study of the implications
of this correlation for the central engines of GRBs may be the subject of
future work.

Given that our clustering is based solely on the prompt emission light curves,
these distinct classes may represent different progenitors, central engines,
emission mechanisms, or a combination of these. To date, evidence from optical
bumps (SNe/KNe) suggests that the collapse of massive stars (collapsars) or
the merging of compact objects (NS/BH/WD) are two possible scenarios for the
origin of GRBs. We locate all the candidates for KNe and SNe on the
connectivity maps. SN-associated GRBs are spread across the map, as evident
from Figure 3. Interestingly, some of these candidates are adjacent to KN-
associated GRBs. Li et al. (2023) also reported similarity in the temporal and
spectral properties of some of the KN- and SN-associated GRBs. This indicates
that despite differences in the progenitors, they might have similar central
engines and/or emission mechanisms, resulting in similar light curves. It is
possible that collapsars and compact binary mergers may result in the
formation of a central engine, such as a hyperaccreting black hole or a
magnetar, which could be responsible for launching a relativistic jet (Duncan
& Thompson, 1992; Usov, 1992; Popham et al., 1999; Narayan et al., 2001;
Metzger, 2019).

###  3.2 Two distinct populations of KN-associated GRBs

Next, we look for the locations of the KNe-associated GRBs in the Fermi/GBM
map. This is motivated by the findings of Dimple et al. (2023a), who argued
for two distinct populations of GRBs associated with KNe in the Swift/BAT
data. Our KN sample comprises four confirmed cases of KN-associated GRBs:
160821B (Kasliwal et al., 2017; Lamb et al., 2019; Troja et al., 2019),
170817A (Abbott et al., 2017; Goldstein et al., 2017; Valenti et al., 2017),
211211A (Rastinejad et al., 2022; Troja et al., 2022), and 230307A (Dichiara
et al., 2023; Sun et al., 2023; Levan et al., 2024).

The results are presented in Figure 3, which interestingly confirms the
existence of two clusters of KNe-associated GRBs based on Swift/BAT. Unlike
Swift/BAT, Fermi/GBM has observed GRB 170817A, known to be a BNS merger from
GW observations (Abbott et al., 2017) with a confirmed detection of a KN.
Therefore, we find evidence favouring the conjecture of Dimple et al. (2023a)
that the low fluence, small duration cluster may correspond to BNS mergers.
Similarly, GRBs 170817A and 160821B are closely located in the same cluster,
implying they could possibly originate from the same kind of progenitor
systems. Dimple et al. (2023a) argued that the second KNe-associated cluster
could potentially be arising from NS-BH mergers. As GRB 230307A is not
detected by Swift/BAT, the proximity of GRBs 211211A and 230307A in the
Fermi/GBM map suggests the similarity in their light curves and the
progenitors. They are lying in a close vicinity in the low-dimensional
embeddings obtained using deep learning on Fermi/GBM data (Negro et al.,
2024). An independent study of the temporal and spectral properties of these
GRBs also suggests that they might have the same origin (Peng et al., 2024).
However, multimessenger observation of a NS-BH system is necessary to confirm
the prediction that GRBs in this cluster are products of NS-BH mergers. For
the typical NS masses, the BH must be of mass $\sim 3M_{\odot}$ for the NS-BH
to produce a GRB. Recent observation of an NS-BH merger GW230529 (Abac et al.,
2024) with a BH between $2.5-4.5M_{\odot}$ indicates there may be a population
of this low mass BH in the universe. With all these, the first multimessenger
NS-BH merger will revolutionize our understanding of compact objects and GRBs.

##  4 summary and future aspects

ML algorithms are slated to play an important role in classifying GRBs, which
is an open issue. Our current study looked into the classification of
Fermi/GBM GRBs using the latest catalog. Listed below are the three key
takeaways from this work.

  1. 1.

An ML-based classification of the GRB light curves of Fermi/GBM suggests five
distinct clusters of GRBs. Independent analyses in the case of Swift/BAT and
BATSE also find the number of clusters to be five, potentially hinting there
is something unique and fundamental about this number.

  2. 2.

Of these five clusters, two distinct clusters contain KNe-associated GRBs. One
cluster contains GRB 170817A, which is associated with a BNS merger. This
strengthens the case for these two clusters to be associated with BNS and NS-
BH mergers as proposed by Dimple et al. (2023a). However, there are a couple
of SNe-associated GRBs in these clusters which may possibly indicate the
universality of the central engines of GRBs regardless of their progenitors.

  3. 3.

Despite the use of fluence-normalized light curves, the ML-based clusters are
well-separated in the fluence-$T_{\rm 90}$ plane, hinting that fluence played
a pivotal role in clustering and may have a strong correlation with the
duration and complexities or structures in the light curve of the GRBs.

We next discuss the implications of these findings. Even though the Fermi and
Swift datasets contain mostly distinct events, with only a small overlap,
their clustering patterns are remarkably similar. This consistency highlights
the robustness of the ML-based classification. Further, KNe-associated GRBs
are tightly clustered at two locations in both datasets. These clustering
patterns can be used for early identification of such GRBs. This may be
achieved by reanalyzing the existing data with the light curve of a new GRB in
question and asking whether it lies in one of the two clusters. Such an
algorithm can be developed and applied to every GRB, assigning a probability
of a KNe association. This should enable dedicated deep searches for KNe
counterparts associated with these GRBs, even without a GW observation. This
will require developing a statistical framework (and a related analysis
pipeline) that will carry out this exercise in a self-consistent manner.
Moreover, these predictions will save the telescope time as it is impossible
to perform a targeted search of all the GRBs, given the faint and rapid decay
nature of KNe. The recent discovery of KNe associated with long GRBs further
complicates the search, as previously thought to be exclusive to short GRBs.

Secondly, assuming the two clusters correspond to BNS and NS-BH mergers, the
rate of GRBs in these clusters could be translated into rates of BNS and NS-BH
mergers under reasonable assumptions about the detector’s sky coverage,
beaming angle range, etc. This may be compared against the rates of BNS and
NS-BH mergers from GW observations to see if these two estimates are
consistent (Kapadia et al., 2024). As the next-generation GW detectors such as
Voyager (Adhikari et al., 2020), Cosmic Explorer (Evans et al., 2021), and the
Einstein Telescope (Sathyaprakash et al., 2012; Kalogera et al., 2021) that
are capable of detecting BNS and NS-BH mergers even at higher redshifts are
operational, this method can give us valuable insights about the relation
between GRBs and compact binary mergers.

Finally, as a caveat, one should note that we are not using redshift-corrected
light curves, as the number of GRBs with redshift estimation is limited. In
the future, with a large dataset with known redshifts, we can explore how time
dilation affects this clustering. Further, to advance this paradigm, it may be
beneficial to employ the ML approach on the afterglows of these GRBs and study
what we can learn from combining the multiband data.

## acknowledgments

We thank the Fermi/GBM team for providing extensive data through public
archives. Dimple thanks P. Sreekumar, Ranjeev Mishra, Siddharth Pritam, Ben
Gompertz for fruitful discussions. K.G.A. and Dimple acknowledge the
Swarnajayanti grant DST/SJF/PSA-01/2017-18 of the Department of Science and
Technology and SERB, India, and support from Infosys Foundation. KM and K.G.A.
acknowledge the support from BRICS grant DST/ICD/BRICS/Call-5/CoNMuTraMO/2023
(G) funded by the Department of Science and Technology (DST), India.

## References

  * Abac et al. (2024) Abac, A. G., et al. 2024.  https://arxiv.org/abs/2404.04248
  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 848, L13, doi: 10.3847/2041-8213/aa920c
  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, The Astrophysical Journal, 848, L12, doi: 10.3847/2041-8213/aa91c9
  * Acuner & Ryde (2018) Acuner, Z., & Ryde, F. 2018, MNRAS, 475, 1708, doi: 10.1093/mnras/stx3106
  * Adhikari et al. (2020) Adhikari, R. X., Arai, K., Brooks, A. F., et al. 2020, Classical and Quantum Gravity, 37, 165003, doi: 10.1088/1361-6382/ab9143
  * Ahumada et al. (2021) Ahumada, T., Singer, L. P., Anand, S., et al. 2021, Nature Astronomy, 5, 917, doi: 10.1038/s41550-021-01428-7
  * Athey et al. (2019) Athey, T. L., Liu, T., Pedigo, B. D., & Vogelstein, J. T. 2019, arXiv e-prints, arXiv:1909.02688.  https://arxiv.org/abs/1909.02688
  * Becerra et al. (2023) Becerra, R. L., Troja, E., Watson, A. M., et al. 2023, MNRAS, 522, 5204, doi: 10.1093/mnras/stad1372
  * Berger et al. (2013) Berger, E., Fong, W., & Chornock, R. 2013, ApJL, 774, L23, doi: 10.1088/2041-8205/774/2/L23
  * Bloom et al. (2002) Bloom, J. S., Kulkarni, S. R., & Djorgovski, S. G. 2002, Astrophysical Journal, 123, 1111, doi: 10.1086/338893
  * Bromberg et al. (2013) Bromberg, O., Nakar, E., Piran, T., & Sari, R. 2013, Astrophysical Journal, 764, doi: 10.1088/0004-637X/764/2/179
  * Cano et al. (2017) Cano, Z., Izzo, L., de Ugarte Postigo, A., et al. 2017, Astronomy & Astrophysics, 605, A107, doi: 10.1051/0004-6361/201731005
  * Chattopadhyay & Maitra (2017) Chattopadhyay, S., & Maitra, R. 2017, MNRAS, 469, 3374, doi: 10.1093/mnras/stx1024
  * Dichiara et al. (2023) Dichiara, S., Tsang, D., Troja, E., et al. 2023, ApJ, 954, L29, doi: 10.3847/2041-8213/acf21d
  * Dimple et al. (2023a) Dimple, Misra, K., & Arun, K. G. 2023a, ApJ, 949, L22, doi: 10.3847/2041-8213/acd4c4
  * Dimple et al. (2023b) Dimple, Misra, K., & Yadav, L. 2023b, arXiv e-prints, arXiv:2309.12788, doi: 10.48550/arXiv.2309.12788
  * Dimple et al. (2022a) Dimple, Misra, K., Ghosh, A., et al. 2022a, Journal of Astrophysics and Astronomy, 43, 39, doi: 10.1007/s12036-022-09822-x
  * Dimple et al. (2022b) Dimple, Misra, K., Kann, D. A., et al. 2022b, MNRAS, 516, 1, doi: 10.1093/mnras/stac2162
  * Duncan & Thompson (1992) Duncan, R. C., & Thompson, C. 1992, ApJ, 392, L9, doi: 10.1086/186413
  * Evans et al. (2021) Evans, M., Adhikari, R. X., Afle, C., et al. 2021, A Horizon Study for Cosmic Explorer: Science, Observatories, and Community.  https://arxiv.org/abs/2109.09882
  * Fishman & Meegan (1995) Fishman, G. J., & Meegan, C. A. 1995, araa, 33, 415, doi: 10.1146/annurev.aa.33.090195.002215
  * Fynbo et al. (2006) Fynbo, J. P. U., Watson, D., Thöne, C. C., et al. 2006, Nature, 444, 1047, doi: 10.1038/nature05375
  * Garcia-Cifuentes et al. (2023) Garcia-Cifuentes, K., Becerra, R. L., De Colle, F., Cabrera, J. I., & del Burgo, C. 2023, arXiv e-prints, arXiv:2304.08666, doi: 10.48550/arXiv.2304.08666
  * Goldstein et al. (2022) Goldstein, A., Cleveland, W. H., & Kocevski, D. 2022, Fermi GBM Data Tools: v1.1.1.  https://fermi.gsfc.nasa.gov/ssc/data/analysis/gbm
  * Goldstein et al. (2017) Goldstein, A., Veres, P., Burns, E., et al. 2017, The Astrophysical Journal, 848, L14, doi: 10.3847/2041-8213/aa8f41
  * Gompertz et al. (2022) Gompertz, B. P., Ravasio, M. E., Nicholl, M., et al. 2022, Nature Astronomy, doi: 10.1038/s41550-022-01819-4
  * Hjorth et al. (2003) Hjorth, J., Sollerman, J., Møller, P., et al. 2003, Nature, 423, 847, doi: 10.1038/nature01750
  * Horváth et al. (2001) Horváth, I., Balázs, L. G., Mészáros, P., Bagoly, Z., & Mészáros, A. 2001, in Gamma-ray Bursts in the Afterglow Era, ed. E. Costa, F. Frontera, & J. Hjorth, 66, doi: 10.1007/10853853_17
  * Jespersen et al. (2020a) Jespersen, C. K., Severin, J. B., Steinhardt, C. L., et al. 2020a, ApJ, 896, L20, doi: 10.3847/2041-8213/ab964d
  * Jespersen et al. (2020b) —. 2020b, ApJL, 896, L20, doi: 10.3847/2041-8213/ab964d
  * Kalogera et al. (2021) Kalogera, V., Sathyaprakash, B. S., Bailes, M., et al. 2021, The Next Generation Global Gravitational Wave Observatory: The Science Book.  https://arxiv.org/abs/2111.06990
  * Kapadia et al. (2024) Kapadia, S., Dimple, Jain, D., et al. 2024, In Preparation
  * Kasliwal et al. (2017) Kasliwal, M. M., Korobkin, O., Lau, R. M., Wollaeger, R., & Fryer, C. L. 2017, ApJL, 843, L34, doi: 10.3847/2041-8213/aa799d
  * Kouveliotou et al. (1993) Kouveliotou, C., Meegan, C. A., Fishman, G. J., et al. 1993, ApJ, 413, L101, doi: 10.1086/186969
  * Lamb et al. (2019) Lamb, G. P., Tanvir, N. R., Levan, A. J., et al. 2019, The Astrophysical Journal, 883, 48, doi: 10.3847/1538-4357/ab38bb
  * Levan et al. (2024) Levan, A. J., Gompertz, B. P., Salafia, O. S., et al. 2024, Nature, 626, 737, doi: 10.1038/s41586-023-06759-1
  * Li et al. (2023) Li, Q. M., Zhang, Z. B., Han, X. L., et al. 2023, MNRAS, 524, 1096, doi: 10.1093/mnras/stad1648
  * Lien et al. (2016) Lien, A., Sakamoto, T., Barthelmy, S. D., et al. 2016, ApJ, 829, 7, doi: 10.3847/0004-637X/829/1/7
  * McInnes et al. (2018) McInnes, L., Healy, J., & Melville, J. 2018, arXiv e-prints, arXiv:1802.03426.  https://arxiv.org/abs/1802.03426
  * Mehta & Iyyani (2024) Mehta, N., & Iyyani, S. 2024, arXiv e-prints, arXiv:2402.15260, doi: 10.48550/arXiv.2402.15260
  * Metzger (2019) Metzger, B. D. 2019, Living Reviews in Relativity, 23, 1, doi: 10.1007/s41114-019-0024-0
  * Narayan et al. (2001) Narayan, R., Piran, T., & Kumar, P. 2001, ApJ, 557, 949, doi: 10.1086/322267
  * Negro et al. (2024) Negro, M., Cibrario, N., Burns, E., et al. 2024, Prompt GRB recognition through waterfalls and deep learning.  https://arxiv.org/abs/2406.03643
  * Paczynski (1986) Paczynski, B. 1986, The Astrophysical Journal, 308, L43, doi: 10.1086/184740
  * Peng et al. (2024) Peng, Z.-Y., Chen, J.-M., & Mao, J. 2024, arXiv e-prints, arXiv:2404.17913, doi: 10.48550/arXiv.2404.17913
  * Petrosian & Dainotti (2024) Petrosian, V., & Dainotti, M. G. 2024, ApJ, 963, L12, doi: 10.3847/2041-8213/ad2763
  * Popham et al. (1999) Popham, R., Woosley, S. E., & Fryer, C. 1999, ApJ, 518, 356, doi: 10.1086/307259
  * Qin & Chen (2013) Qin, Y.-P., & Chen, Z.-F. 2013, Monthly Notices of the Royal Astronomical Society, 430, 163, doi: 10.1093/mnras/sts547
  * Rastinejad et al. (2022) Rastinejad, J. C., Gompertz, B. P., Levan, A. J., et al. 2022, arXiv e-prints, arXiv:2204.10864.  https://arxiv.org/abs/2204.10864
  * Salmon et al. (2022) Salmon, L., Hanlon, L., & Martin-Carrillo, A. 2022, Galaxies, 10, 77, doi: 10.3390/galaxies10040077
  * Sathyaprakash et al. (2012) Sathyaprakash, B., Abernathy, M., Acernese, F., et al. 2012, Classical and Quantum Gravity, 29, 124013, doi: 10.1088/0264-9381/29/12/124013
  * Steinhardt et al. (2023) Steinhardt, C. L., Mann, W. J., Rusakov, V., & Jespersen, C. K. 2023, arXiv e-prints, arXiv:2301.00820.  https://arxiv.org/abs/2301.00820
  * Stern et al. (1997) Stern, B., Poutanen, J., & Svensson, R. 1997, ApJ, 489, L41, doi: 10.1086/310947
  * Sun et al. (2023) Sun, H., Wang, C. W., Yang, J., et al. 2023, arXiv e-prints, arXiv:2307.05689, doi: 10.48550/arXiv.2307.05689
  * Tanvir et al. (2013) Tanvir, N. R., Levan, A. J., Fruchter, A. S., et al. 2013, Nature, 500, 547, doi: 10.1038/nature12505
  * Troja et al. (2019) Troja, E., Castro-Tirado, A. J., Becerra González, J., et al. 2019, MNRAS, 489, 2104, doi: 10.1093/mnras/stz2255
  * Troja et al. (2022) Troja, E., Fryer, C. L., O’Connor, B., et al. 2022, Nature, 612, 228, doi: 10.1038/s41586-022-05327-3
  * Usov (1992) Usov, V. V. 1992, Nature, 357, 472, doi: 10.1038/357472a0
  * Valenti et al. (2017) Valenti, S., Sand, D. J., Yang, S., et al. 2017, ApJL, 848, L24, doi: 10.3847/2041-8213/aa8edf
  * von Kienlin et al. (2020) von Kienlin, A., Meegan, C. A., Paciesas, W. S., et al. 2020, ApJ, 893, 46, doi: 10.3847/1538-4357/ab7a18
  * Woosley & Bloom (2006) Woosley, S. E., & Bloom, J. S. 2006, ARAA, 44, 507, doi: 10.1146/annurev.astro.43.072103.150558
  * Yang et al. (2022) Yang, J., Ai, S., Zhang, B.-B., et al. 2022, Nature, 612, 232, doi: 10.1038/s41586-022-05403-8
  * Zhang et al. (2009) Zhang, B., Zhang, B.-B., Virgili, F. J., et al. 2009, ApJ, 703, 1696, doi: 10.1088/0004-637X/703/2/1696

##  Appendix A Characteristics of the clusters

We probe the known properties of GRBs to investigate what these clusters might
represent. We took the values of the parameters $T_{\rm 90}$, fluence, $E_{\rm
p}$ and $\alpha$ from the Fermi/GBM
catalog222https://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermigbrst.html. Figure
A.1 shows the violin plots for $T_{90}$ (upper left panel), fluence (upper
right panel), and $E_{p}$ (lower left panel), $\alpha$ (lower right panel) for
each cluster.

Figure A.1: Violin plots for the known GRB properties corresponding to the
five groups of GRBs identified by ML algorithms.


# Transport in p-wave interacting Fermi gases

Jeff Maki  Department of Physics and HKU-UCAS Joint Institute for Theoretical
and Computational Physics, The University of Hong Kong, Hong Kong, China
Pitaevskii BEC Center, CNR-INO and Dipartimento di Fisica, Università di
Trento, I-38123 Trento, Italy  Tilman Enss  Institut für Theoretische Physik,
Universität Heidelberg, 69120 Heidelberg, Germany

(August 27, 2024)

###### Abstract

The scattering properties of spin-polarized Fermi gases are dominated by
p-wave interactions. Besides their inherent angular dependence, these
interactions differ from their s-wave counterparts as they also require the
presence of a finite effective range in order to understand the low-energy
properties of the system. In this article we examine how the shear viscosity
and thermal conductivity of a three-dimensional spin-polarized Fermi gas in
the normal phase depend on the effective range and the scattering volume in
both the weakly and strongly interacting limits. We show that although the
shear viscosity and thermal conductivity both explicitly depend on the
effective range near resonance, the Prandtl number which parametrizes the
ratio of momentum to thermal diffusivity does not have an explicit interaction
dependence both at resonance and for weak interactions in the low-energy
limit. In contrast to s-wave systems, p-wave scattering exhibits an additional
resonance at weak attraction from a quasi-bound state at positive energies,
which leads to a pronounced dip in the shear viscosity at specific
temperatures.

##  I Introduction

Spin-polarized Fermi gases have become an excellent platform for studying
quantum systems with higher partial-wave interactions. Due to the Pauli
exclusion principle, the leading short-ranged interactions are p-wave in
nature ($\ell=1$, where $\ell$ is the angular momentum quantum number). For
this reason, there have been numerous theoretical investigations into the
physics of p-wave Fermi gases, such as p-wave superconductivity and
topological physics in two-dimensions [1, 2, 3, 5, 6, 4], the three-body loss
rate [7, 8, 9, 10], and its application to one-dimensional physics [11, 12,
13, 14, 16, 15, 17]. Such physics is not purely theoretical, but has become
the focus of recent experimental investigations both in 40K and in 6Li, as
there is a p-wave Feshbach resonance [18, 19, 20, 22, 23, 25, 26, 27, 28, 29,
21, 24].

The main difference between p-wave ($\ell=1)$ and s-wave ($\ell=0$) scattering
in 3D is the presence of the centrifugal barrier. The centrifugal barrier
limits the overlap between two-body bound states and scattering states to a
region at short inter-particle distances, making the resonance inherently
narrow. As a corollary, the two-body bound state is quite long-lived near
resonance, even for positive energies when the two-body bound state is inside
the scattering continuum, i.e., the bound state becomes a long lived quasi-
bound state. In other words, the wavefunction for the two-body bound state
remains localized even near resonance, in contrast to s-wave systems where the
two-body bound state size approaches infinity as one approaches resonance.

The presence of the centrifugal barrier and the narrowness of the resonance
are related to the relevancy of the effective range in the scattering
amplitude. For p-wave scattering in 3D, the inverse p-wave scattering
amplitude is:

$f_{\ell=1}^{-1}=-i-\frac{1}{p^{3}v}-\frac{R}{p}+O(p)$ (1)

where $p$ is the magnitude of the relative momentum, related to the relative
scattering energy, $E=p^{2}/m$ with $m$ the single-particle mass. We also
define $v$ as the 3D p-wave scattering volume with units of volume, and $R$ as
the p-wave effective range parameter with units of momentum. The first term
represents the unitary scattering, while the second and third terms define the
scattering parameters. As one can see from the effective range expansion, both
the scattering volume and effective range terms in the scattering amplitude
are parametrically more important than the unitary term at low energies.
Furthermore, one can show that the low-energy limit, $E\ll R^{2}/m$, and the
zero-range limit $R\to 0$, can not be taken simultaneously, in
contradistinction to s-wave interactions. Thus the low-energy scattering
physics depends on both the scattering volume and the effective range.

Since the effective range is a relevant quantity to understand the low-energy
scattering, the energetics and dynamics will also crucially depend on the
effective range. Previous studies have examined this in the context of the
necessity of two thermodynamic contacts in describing the energetics in the
normal phase [30, 31], the Landau liquid parameters [32], and the three-body
recombination rate [8, 10]. In these studies, the effective range was found to
be important describing the leading behaviour. This ought to be compared to
s-wave physics where the effective range merely adds a perturbative correction
[33, 34, 35, 36].

In terms of transport, a previous study examined the bulk viscosity [37] for
the p-wave Fermi gas. There it was shown that the effective range produces a
finite bulk viscosity at resonance (i.e., when $v^{-1}=0$) proportional to
$\zeta\propto T^{5/2}/R^{2}$. In the weakly interacting limit, the bulk
viscosity is proportional to $v^{2}T^{9/2}$. If one tries to take the zero-
range limit, the bulk viscosity at resonance diverges, which is a hallmark of
the relevancy of the effective range. This is in contrast to the s-wave case
where the bulk viscosity only depends on the s-wave scattering length, $a$, in
the strongly interacting limit, and vanishes at resonance, $a^{-1}=0$ [40, 39,
41, 38, 42, 43].

An interesting open question is how the remaining transport coefficients, the
shear viscosity and thermal conductivity, depend on the scattering parameters
in both the weakly and strongly interacting limit. For s-wave interactions, it
is known that these two quantities become divergent near the non-interacting
point ($a=0$) as $a^{-2}$, while near resonance the shear viscosity and
thermal conductivity depend only on the equation of state, i.e., the density
and the temperature [44, 45, 46, 47, 48, 49, 50, 51, 41, 52]. It is unclear
how this picture is modified for p-wave Fermi gases in 3D, even in the
experimentally applicable limit of a small but finite effective range.

In this work we consider this issue and evaluate the shear viscosity, $\eta$,
and thermal conductivity, $\kappa$, using the kinetic theory approach [53,
54]. We find that the shear viscosity and thermal conductivity scattering
times explicitly depend on the scattering parameters, even at resonance.
Although the shear viscosity and thermal conductivity explicitly depend on the
interaction parameters, we show that the Prandtl number, which describes the
ratio of momentum and thermal diffusion, approaches a universal constant in
these two limits that does not explicitly depend on the interaction, similar
to the case of s-wave physics.

The remainder of this article is organized as follows. In Sec. II we present
the two-body scattering properties for spin-polarized Fermi gases and obtain
the two-body T-matrix. From there we give a brief overview of the kinetic
theory approach and how it applies to the shear viscosity and thermal
conductivity in Sec. III. We then present the results for the shear and
thermal scattering times in Sec. IV for arbitrary values of the scattering
volume. For negative scattering volumes, when the bound state becomes a long-
lived quasi-bound state we find non-monotonic behaviour for the transport
properties which is further discussed in Sec. V. From there we discuss the
Prandtl number in Sec. VI, and finally conclude our discussions in Sec. VII.

##  II Two-Body p-Wave Scattering

In this article we consider a single-channel model for a spin-polarized Fermi
gas with p-wave interactions:

$\displaystyle H$ $\displaystyle=\int d^{3}{\bf x}\frac{1}{2}\nabla_{\bf
x}\psi^{\dagger}({\bf x})\nabla_{\bf x}\psi({\bf x})$ $\displaystyle+\int
d^{3}{\bf x}\frac{g}{4}\psi^{\dagger}({\bf x})\overleftrightarrow{\nabla}_{\bf
x}\psi^{\dagger}({\bf x})\psi({\bf x})\overleftrightarrow{\nabla}_{\bf
x}\psi({\bf x})$ (2)

where $\psi^{(\dagger)}({\bf x})$ is the annihilation (creation) operator for
spin-polarized Fermions, the bidirectional gradient is
$\overleftrightarrow{\nabla}_{\bf x}=(\overleftarrow{\nabla}_{\bf
x}-\overrightarrow{\nabla}_{\bf x})/2$, $g$ is the p-wave coupling constant,
and we have set $\hbar$ and the atomic mass, $m$, to unity.

As is custom, we renormalize this theory by examining the two-body scattering.
Consider the two-body T-matrix between states with center-of-mass momentum
${\bf Q}$, relative momenta ${\bf p}$ and ${\bf q}$, $|{\bf Q}/2\pm{\bf p\
(q)}\rangle$, and a total complex frequency $E$. In the presence of the many-
body background the T-matrix has the form [31, 37]:

$\displaystyle\langle{\bf\frac{Q}{2}\pm p}|$
$\displaystyle\mathcal{T}|{\bf\frac{Q}{2}\pm q}\rangle={\bf p\cdot q}\,T({\bf
Q},E)$ $\displaystyle T^{-1}({\bf Q},E)$
$\displaystyle=\frac{1}{24\pi}\left[\frac{1}{v}+\left(E-\frac{Q^{2}}{4}+2\mu\right)R\right.$
$\displaystyle\left.+\left(-E+\frac{Q^{2}}{4}-2\mu\right)^{3/2}+A_{\text{mb}}({\bf
Q},E)\right]$ (3)

where $\mu$ is the chemical potential, and $A_{\text{mb}}({\bf Q},E)$ is the
contribution due to the many-body background which is reported in Appendix A.
Such a term captures the effects of Fermi blocking in the intermediate
scattering states on an equal footing with the Fermi factors in the collision
integral (see below, Eq. (LABEL:eq:colision_integral)).

The T-matrix in the first line of Eq. (3) splits into two pieces [55]. The
first piece is the form factor of the p-wave interaction potential, ${\bf
p\cdot q}$, see also Eq. (2). The second piece, $T({\bf Q},E)$, describes the
dependence of the scattering on the center-of-mass momentum, ${\bf Q}$, and
total energy, $E$. Due to the presence of the many-body background, the
scattering is no longer Galilean invariant, which is described by the
nontrivial dependence on the center-of-mass momentum in $A_{\text{mb}}({\bf
Q},E)$.

Equation (3) is already renormalized, and both the scattering volume, $v$, and
effective range, $R$, are defined in terms of the coupling constant $g$ and
the ultraviolet cutoff in the theory, $\Lambda$:

$\displaystyle\frac{1}{24\pi v}$
$\displaystyle=\frac{1}{g}-\frac{\Lambda^{3}}{18\pi^{2}},$
$\displaystyle\frac{R}{24\pi}$ $\displaystyle=\frac{\Lambda}{6\pi^{2}}.$ (4)

This renormalization ensures that the p-wave scattering amplitude in the
absence of the many-body background has the form shown in Eq. (1) since:

$f_{\ell=1}=-\frac{p^{3}T({\bf Q},E_{o.s})}{24\pi},$ (5)

where now $E_{o.s}=\frac{Q^{2}}{4}+p^{2}-2\mu+i\delta$ is the on-shell energy
for two-particle scattering.

In this theory there are two two-body bound states at energies defined as the
poles of Eq. (3). The first is a shallow dimer with $E_{b}=-1/(vR)$. For
positive values of the scattering volume, this is a true two-body bound state.
Therefore we call the regime where $E_{b}<0$ the BEC side. For negative values
of the scattering volume, $E_{b}>0$, the dimer is a long-lived quasi-bound
state. We label this side as the BCS side. The second bound state is a deep
dimer with energy $-R^{2}$. Both the shallow and deep dimer bound states are
three-fold degenerate for $\ell=1$. The state with energy $-R^{2}$ is actually
unphysical as it possesses a negative norm [56, 57]. To avoid this issue, we
will work in the low-energy limit where all energy scales $E$ satisfy: $E\ll
R^{2}$, so that the presence of such an unphysical state is unimportant.

We note that often a single-channel model is inadequate for describing p-wave
scattering in terms of the microscopic parameters that parametrize the two-
body interaction potential, for e.g. the van der Waals length. In this study
we are only interested in the general properties of the transport in terms of
the low-energy scattering parameters, $v$ and $R$. For our purposes, the
single-channel model suffices after the proper renormalization of the T-matrix
is taken into account, see Eq. (8). However, for more in-depth knowledge of
the individual p-wave scattering parameters, a two-channel model is required
[31, 37].

##  III Kinetic Theory

In order to calculate the shear viscosity and thermal conductivity we employ
the standard kinetic theory approach using the Boltzmann equation [53, 54]

$\frac{\partial n_{\bf p}}{\partial t}+{\bf p}\cdot\nabla_{\bf x}n_{\bf
p}=I\left[n_{\bf p}\right],$ (6)

where $n_{\bf p}=n_{\bf p}({\bf x},t)$ is the local quasiparticle distribution
function. The collision integral, $I\left[n_{\bf p}\right]$, is defined as

$\displaystyle I$ $\displaystyle\left[n_{\bf p}\right]=\int\frac{d^{3}{\bf
q}}{(2\pi)^{3}}\int\frac{d^{3}{\bf q^{\prime}}}{(2\pi)^{3}}\int\frac{d^{3}{\bf
p^{\prime}}}{(2\pi)^{3}}W({\bf p,q,p^{\prime},q^{\prime}})$
$\displaystyle\times\left[(1-n_{\bf p})(1-n_{\bf q})n_{\bf p^{\prime}}n_{\bf
q^{\prime}}-(1-n_{\bf p^{\prime}})(1-n_{\bf p^{\prime}})n_{\bf p}n_{\bf
q}\right].$

Equation (LABEL:eq:colision_integral) depends on the transition rate $W({\bf
p,q,p^{\prime},q^{\prime}})$ between particles with incoming momenta $\bf p,q$
and outgoing momenta $\bf p^{\prime},q^{\prime}$, defined in terms of the on-
shell T-matrix, $T({\bf Q,p})=T({\bf Q},E_{o.s.})$:

$\displaystyle W$ $\displaystyle({\bf
p,q,p^{\prime},q^{\prime}})=(2\pi)^{4}\delta({\bf p+q-p^{\prime}-q^{\prime}})$
$\displaystyle\times\delta\left(\frac{p^{2}}{2}+\frac{q^{2}}{2}-\frac{p^{\prime
2}}{2}-\frac{q^{\prime 2}}{2}\right)\left(\frac{\bf p-q}{2}\cdot\frac{\bf
p^{\prime}-q^{\prime}}{2}\right)^{2}$ $\displaystyle\times\left|T\left({\bf
p+q},\frac{\bf p-q}{2}\right)\right|^{2}.$ (8)

Following kinetic theory we linearize the Boltzmann equation by writing
$n_{\bf p}=n_{\bf p}^{0}+\delta n_{\bf p}$, where $n_{\bf p}^{0}$ is the local
equilibrium distribution function,

$n_{\bf p}^{0}=\left[e^{\beta\left(({\bf
p-v})^{2}/2-\mu\right)}+1\right]^{-1},$ (9)

that depends on the local inverse temperature, $\beta({\bf x},t)=1/T({\bf
x},t)$, velocity, ${\bf v}({\bf x},t)$ and chemical potential, $\mu({\bf
x},t)$. We have muted the spatial and temporal coordinates of the
thermodynamic variables for simplicity. The correction to the distribution
function in response to an external perturbation is denoted by $\delta n_{\bf
p}$. It is subject to the constraints of conserved number, momentum and
energy,

$\displaystyle 0=\int\frac{d^{3}{\bf p}}{(2\pi)^{3}}\\{1,{\bf
p},\frac{p^{2}}{2}\\}\ \delta n_{\bf p}.$ (10)

It can be written in the form $\delta n_{\bf p}=\beta n_{\bf p}^{0}(1-n_{\bf
p}^{0})\phi({\bf p})$ with

$\displaystyle\phi({\bf p})$ $\displaystyle=-\left[\phi_{ij}({\bf
p-v})\frac{V_{i,j}}{2}-{\boldsymbol{\phi}}({\bf p-v})\cdot\nabla_{\bf
x}\ln\beta\right]$ (11)

and $V_{i,j}=\partial_{i}v_{j}+\partial_{j}v_{i}-2/3\delta_{i,j}\nabla_{\bf
x}\cdot{\bf v}$.

Due to rotational invariance, we only need to consider perturbations of the
form $\phi_{xy}({\bf p})$ (shear) and $\phi_{x}({\bf p})$ (heat current).
Following the kinetic theory approach [53, 42] the linearized Boltzmann
equation leads to equations that determine $\phi_{xy}({\bf p})$ and
$\phi_{x}({\bf p})$:

$\displaystyle p_{x}p_{y}=\int\frac{d^{3}{\bf
q}}{(2\pi)^{3}}\int\frac{d^{3}{\bf p^{\prime}}}{(2\pi)^{3}}\int\frac{d^{3}{\bf
q^{\prime}}}{(2\pi)^{3}}n_{\bf q}^{0}(1-n_{\bf q}^{0})$ $\displaystyle W({\bf
p,q,p^{\prime},q^{\prime}})\left(\phi_{xy}({\bf p})+\phi_{xy}({\bf
q})-\phi_{xy}({\bf p^{\prime}})-\phi_{xy}({\bf q^{\prime}})\right),$
$\displaystyle p_{x}\left(\frac{p^{2}}{2}-w\right)=\int\frac{d^{3}{\bf
q}}{(2\pi)^{3}}\int\frac{d^{3}{\bf p^{\prime}}}{(2\pi)^{3}}\int\frac{d^{3}{\bf
q^{\prime}}}{(2\pi)^{3}}n_{\bf q}^{0}(1-n_{\bf q}^{0})$ $\displaystyle W({\bf
p,q,p^{\prime},q^{\prime}})\left(\phi_{x}({\bf p})+\phi_{x}({\bf
q})-\phi_{x}({\bf p^{\prime}})-\phi_{x}({\bf q^{\prime}})\right).$

In Eq. (LABEL:eq:phix), $w$ is the enthalpy per particle,
$w=(\varepsilon+p)/n$, with $\varepsilon$ as the energy density, $p$ the
pressure, and $n$ the density. The functions $\phi_{xy}({\bf p})$ and
$\phi_{x}({\bf p})$ which solve Eqs. (LABEL:eq:phixy-LABEL:eq:phix) are then
related to the shear viscosity, $\eta$, and thermal conductivity, $\kappa$,
respectively:

$\displaystyle\eta$ $\displaystyle=\beta\int\frac{d^{3}{\bf
p}}{(2\pi)^{3}}n_{\bf p}^{0}(1-n_{\bf p}^{0})p_{x}p_{y}\phi_{xy}({\bf p}),$
(14) $\displaystyle\frac{\kappa}{\beta}$
$\displaystyle=\beta\int\frac{d^{3}{\bf p}}{(2\pi)^{3}}n_{\bf p}^{0}(1-n_{\bf
p}^{0})p_{x}\left(\frac{p^{2}}{2}-w\right)\phi_{x}({\bf p}).$ (15)

To solve Eqs. (LABEL:eq:phixy-15), we expand the functions $\phi_{xy}({\bf
p})$ and $\phi_{x}({\bf p})$ in terms of a set of orthogonal basis functions:

$\displaystyle\phi_{xy}({\bf p})$
$\displaystyle=\beta\sum_{j}c_{j}^{\eta}U_{j}^{\eta}({\bf p}),$ (16)
$\displaystyle\phi_{x}({\bf p})$
$\displaystyle=\beta\sum_{j}c_{j}^{\kappa}U_{j}^{\kappa}({\bf p}),$ (17)

where $c_{j}^{(\eta,\kappa)}$ are expansion coefficients, and the first basis
functions are $U_{1}^{\eta}=p_{x}p_{y}$ and $U_{1}^{\kappa}=p_{x}(p^{2}/2-w)$.
In terms of these modes, Eqs. (LABEL:eq:phixy-15) are more conveniently
expressed as

$\displaystyle\delta_{i,1}$
$\displaystyle=\sum_{j}\mathcal{A}^{\left(\eta,\kappa\right)}_{i,j}c_{j}^{(\eta,\kappa)},$
$\displaystyle\mathcal{A}^{\left(\eta,\kappa\right)}_{i,j}=\frac{\left(U_{i}^{\left(\eta,\kappa\right)},\mathcal{L}U_{j}^{\left(\eta,\kappa\right)}\right)}{\left(U_{1}^{\left(\eta,\kappa\right)},U_{1}^{\left(\eta,\kappa\right)}\right)}.$
(18)

In Eq. (18) we have introduced the inner product

$\displaystyle\left(A,B\right)=\int\frac{d^{3}{\bf p}}{(2\pi)^{3}}n_{\bf
p}^{0}(1-n_{\bf p}^{0})A({\bf p})B({\bf p})$ (19)

as well as the linearized collision integral operator

$\displaystyle\mathcal{L}B$ $\displaystyle=\frac{\beta}{1-n_{\bf
p}^{0}}\int\frac{d^{3}{\bf q}}{(2\pi)^{3}}\int\frac{d^{3}{\bf
p^{\prime}}}{(2\pi)^{3}}\int\frac{d^{3}{\bf q^{\prime}}}{(2\pi)^{3}}n_{\bf
q}^{0}(1-n_{\bf q}^{0})$ $\displaystyle W({\bf
p,q,p^{\prime},q^{\prime}})\left(B({\bf p})+B({\bf q})-B({\bf
p^{\prime}})-B({\bf q^{\prime}})\right).$ (20)

In this notation the shear viscosity and thermal conductivity have a simple
form:

$\displaystyle\eta$
$\displaystyle=\beta^{2}\left(\mathcal{A}^{\eta}\right)^{-1}_{1,1}\left(U_{1}^{\eta},U_{1}^{\eta}\right),$
$\displaystyle\frac{\kappa}{\beta}$
$\displaystyle=\beta^{2}\left(\mathcal{A}^{\kappa}\right)^{-1}_{1,1}\left(U_{1}^{\kappa},U_{1}^{\kappa}\right).$
(21)

The quantity $\beta\left(\mathcal{A}^{\eta,\kappa}\right)^{-1}_{1,1}$ has
units of time and defines the shear scattering time, $\tau_{\eta}$, and the
thermal scattering time, $\tau_{\kappa}$ respectively:

$\tau_{\eta,\kappa}=\beta\left(\mathcal{A}^{\eta,\kappa}\right)^{-1}_{1,1}.$
(22)

In principle, as the number of basis modes increases, the accuracy of the
calculation improves. The case of a single mode is special, and is equivalent
to the relaxation time approximation (RTA) [53, 54]. For our purposes we will
work with two basis modes as that provides a drastic improvement in the
calculation. This is discussed further below. We label this as beyond the
relaxation time approximation (BRTA).

##  IV Shear and Thermal Scattering Times

The linearized Boltzmann equation (18) is amenable to numerical calculation of
the shear and thermal scattering times, $\tau_{\eta}$ and $\tau_{\kappa}$ via
Eq. (22). For our purposes we will consider two basis modes both on the BCS
and on the BEC side. The details of the calculation are shown in Appendix B,
while we only discuss the results here.

The shear and thermal scattering times are presented in Fig. 1 as functions of
fugacity $z$ for various values of the binding energy $\beta
E_{b}=-\beta/(vR)$. The dashed straight lines correspond to the high-
temperature limit where $\beta/\tau_{\eta,\kappa}\propto z$. As the
temperature is lowered (larger $z$), the deviations from the high-temperature
limit become more pronounced, and we find that the shear and thermal
scattering rates become smaller than the predictions from the high-temperature
theory. This reduction of scattering is a consequence of Fermi blocking.

At resonance, $\beta E_{b}=0$, the shear and thermal scattering times are
proportional to $\beta R^{2}$ in the low-energy limit, up to a correction of
order $O(1)$, as can be seen from Appendix B. This is a consequence of
defining the low-energy physics as: $\\{T,E,\mu,...\\}\ll R^{2}$. This leads
us to an important result: the presence of the effective range is necessary
for understanding the shear viscosity and the thermal conductivity in the
strongly interacting limit. Such a conclusion is not necessarily obvious, as
one might expect that for a small effective range the transport would
predominantly depend on the equation of state since scale symmetry is only
broken slightly. However, this is not the case as the low-energy constraint
renders the effective range a relevant quantity. This should be contrasted to
the spin-1/2 s-wave Fermi gas, where the resonant scattering times are only
functions of the temperature and density because of the scale symmetry.

In the weakly interacting limit, a similar analysis shows that the shear and
thermal scattering times are proportional to $v^{-2}$. This is the standard
result for weakly interacting systems; the scattering times become infinite as
one approaches the non-interacting limit.

Figure 1: Shear and thermal scattering rates as functions of fugacity, $z$. In
all calculations we set the effective range $\beta^{1/2}R=100$. The lines from
bottom to top indicate the direction towards resonance, $|\beta E_{b}|\to 0$.
The solid lines are the full numerical solution to the Boltzmann equation,
while the dashed lines represent the corresponding high-temperature
approximation which is linear in $z$.

Although the shear and thermal scattering times have explicit interaction
dependencies in the strongly and weakly interacting limits, the two scattering
times actually have the same interaction dependencies to leading order. This
is readily seen in Appendix B where we provide explicit formulas in the high-
temperature limit. This leads us to our second main result: since both
scattering times have the same leading-order dependence on the scattering
parameters both in the strongly and the weakly interacting limits, the ratio
of the two scattering times in the strongly and weakly interacting limits is
only a function of the equation of state:
$\tau_{\eta}/\tau_{\kappa}=\mathcal{F}\left(\beta\mu\right)+O(1/(\beta
R^{2}))$ with a dimensionless function $\mathcal{F}(x)$. Such behaviour also
occurs for spin-1/2 Fermi gases with s-wave interactions, but this is because
scale symmetry requires that the scattering times themselves are only
functions of the equation of state. In the p-wave case, scale symmetry is
still broken and thus the individual scattering times depend on the
interactions, but their ratio will not depend on the scattering volume or
effective range at leading order.

Let us first consider the results for the ratio of the shear to thermal
scattering times in the high-temperature limit shown in Fig. 2 as a function
of the bound state energy, $E_{b}$. In the RTA, the ratio of the scattering
times is independent of the fugacity and the scattering parameters for
arbitrary interaction strength, $\tau_{\eta}/\tau_{\kappa}=2/3$, see Appendix
B. In the BRTA, the ratio of the scattering times becomes interaction
dependent. In the strongly and weakly interacting limits, the ratios can be
evaluated analytically to give:

Figure 2: BEC-BCS Crossover for the ratio of scattering times,
$\tau_{\eta}/\tau_{\kappa}$ in the high-temperature limit with
$\beta^{1/2}R=100$. In this limit the ratio of scattering times coincides with
the Prandtl number. The left hand side corresponds to the BCS side, while the
right hand side is the BEC side. The solid, dashed, and dashed-dotted lines
correspond to the RTA prediction, the BRTA prediction at resonance, and the
BRTA prediction for weak interactions, respectively. The vertical dashed lines
correspond to the points where the Prandtl number saturates the RTA value,
$\beta E_{b}\approx 0,7/2$.

$\displaystyle\left.\frac{\tau_{\eta}}{\tau_{\kappa}}\right|_{\textrm{res.}}$
$\displaystyle=\frac{2}{3}-\frac{14}{2727}\approx 0.662$
$\displaystyle\left.\frac{\tau_{\eta}}{\tau_{\kappa}}\right|_{\textrm{weak
int.}}$ $\displaystyle=\frac{2}{3}-\frac{250}{3729}\approx 0.600$ (23)

In these limits the ratio of the scattering times does not explicitly depend
on the interaction parameters, $v$ and $R$, the fugacity, $z$, and the inverse
temperature, $\beta$. The individual values of $\tau_{\eta}$ and
$\tau_{\kappa}$ are shown in Appendix B.

On the BEC side (right) the behaviour of the ratio of scattering times is
monotonic, while on the BCS side (left) there are two point where it saturates
the RTA value. This can be understood as a consequence of the quasi-long-lived
bound state, see Appendix B. Specifically these resonances occur when the off-
diagonal matrix elements of $\mathcal{A}^{(\eta,\kappa)}$ vanish. To leading
order this occurs when the two-body bound state energy is approximately:
$\beta E_{b}\approx 0,7/2$. Since this resonance requires the two-body bound
state energy to be positive, this secondary resonance where the ratio of
scattering time saturates the RTA value can not occur for s-wave interacting
systems.

As one lowers the temperature however, we do expect deviations to occur from
the Fig. 2. In particular, we expect that the ratio of the scattering times to
leading order in $1/(\beta R^{2})$ will be a dimensionless function of the
equation of state, but not of the interaction strength, in the strongly and
weakly interacting limit. In Fig. 3 we examine the ratio of the two scattering
times both a) as a function of fugacity and b) as a function of binding energy
$E_{b}$ on the BEC side. The universal limits in Eq. (23) are shown as the
dashed and dash-dotted lines, respectively. The black solid line corresponds
to the RTA result of $2/3$. The two major trends are an increase of the ratio
of scattering times as the temperature is lowered, and a decrease in the ratio
as one goes to the weakly interacting limit. On the BCS side, $v<0$, we find
that the many-body corrections are not as important as the scattering physics
is highly dominated by the quasi-long-lived bound state. Hence the physics is
accurately captured by the high-temperature physics. Thus even in the presence
of the many-body background, there will still be a non-monotonic behaviour of
the ratio of the scattering times, with the BRTA saturating to the RTA value
at around $\beta E_{b}\approx 7/2$.

Figure 3: Ratio of shear to thermal scattering times as a function of the
fugacity $z$ (top), and the binding energy $E_{b}$ (bottom). The black solid
line, dashed line, and dashed-dotted line correspond to the high-temperature
results in the RTA, the resonant limit in the BRTA, and the weakly interacting
limit in the BRTA, respectively. The lines from bottom to top indicate the
direction towards: top) resonance, $|\beta E_{b}|\to 0$, and bottom) the low-
temperature limit.

##  V Minimum in the Shear Viscosity and Thermal Conductivity on the BCS Side

The presence of the quasi-bound state and the non-monotonic behaviour of the
ratio of the scattering times on the BCS side have important consequences for
the shear viscosity and thermal conductivity. In particular, if the average
energy of the atoms, which is proportional to $T$, is of the order of the
quasi-bound state energy, the scattering will become maximal and equivalently
the transport coefficient will exhibit a minimum. Thus it is instructive to
look at the shear viscosity at fixed density and variable temperature, as we
expect a dip to occur on the BCS side when the temperature is comparable to
the quasi-bound state energy.

We investigated this issue by calculating the shear viscosity at fixed density
as a function of temperature within the RTA. Corrections from the BRTA are
smaller than the Fermi blocking effect included in the RTA at lower
temperatures and do not produce qualitative differences, see Fig. 7 of
Appendix B. In order to address the thermodynamics, we have assumed a
noninteracting equation of state of spin-polarized Fermions. Such an
approximation is reasonably valid in the BCS limit where no molecules exist,
in contrast to the BEC limit discussed below where a Bose-Fermi model is used
[31]. The results of this calculation are shown in Fig. 4 for $E_{b}=20E_{F}$,
$E_{b}=10E_{F}$ and $E_{b}=0$ (resonance).

As one can see from Fig. 4, there is a dip in the shear viscosity due to the
resonant scattering at the quasi-bound state energy. The minimum in the shear
viscosity occurs roughly for $T/T_{F}\approx 0.2E_{b}/E_{F}$. At resonance,
there is no quasi-bound state and hence there is no minimum of the shear
viscosity in the normal state. It is interesting to note that for p-wave
scattering this minimum in the shear viscosity appears for $T>T_{F}$, while
for s-wave Fermi gases it occurs at lower temperatures [48]. Similarly, when
the shear viscosity is expressed in units of entropy density $s$, it will also
have a minimum. For example when $E_{b}/E_{F}=5$, $\eta/s$ reaches a minimum
value of approximately $0.5\hbar/k_{B}$ at a temperature of $T\approx T_{F}$,
comparable to the value found for the s-wave unitary Fermi gas [48] and
slightly larger than the Kovtun-Son-Starinets bound [58].

Figure 4: Shear viscosity as a function of temperature for fixed density and
various values of the binding energy, $E_{b}$, on the BCS side ($E_{b}>0$). In
this figure, $T_{F}$ and $E_{F}$ are the Fermi temperature and energy, while
the solid and dashed lines correspond to the RTA calculation and the high-
temperature limit respectively. The dip in the shear viscosity is associated
with an increased scattering at the quasi-bound state energy. This dip in the
shear viscosity does not appear at resonance, $E_{b}=0$, as zero-energy
scattering is strongly suppressed for p-wave interactions. Similar physics
will occur for the thermal conductivity.

##  VI Prandtl Number

Given the ratio of the shear to thermal scattering times, it is
straightforward to calculate the Prandtl number. The Prandtl number is the
ratio of the shear to thermal diffusivity, defined as

$\Pr=\frac{\eta}{\kappa}\frac{C_{p}}{n},$ (24)

where $C_{p}$ is the specific heat at constant pressure per unit volume, and
$n$ is the density. In this work we calculate the shear viscosity and thermal
conductivity according to Eq. (21), while we calculate the specific heat at
constant pressure and density using the Bose-Fermi model, see below.

At high temperatures, $C_{p}/n=5/2$ while
$\eta/\kappa=2\tau_{\eta}/5\tau_{\kappa}$, so that
$\Pr=\tau_{\eta}/\tau_{\kappa}$. This was presented in Fig. 2 for the whole
BEC-BCS crossover. As one can see, the Prandtl number is a nonmonotonic
function of $E_{b}$ on the BCS side, while it is monotonic on the BEC side.

We evaluate the thermodynamics using the Nozieres-Schmitt-Rink (NSR)
approximation [59]. This scheme involves calculating the free energy using the
in-medium T-matrix. In this way we calculate the free energy and the
scattering times consistently at the same level of approximation. Such a
calculation was done previously in the context of p-wave Fermi gases using a
two-channel model [31], but we have verified that their predictions are
equivalent to the single-channel model after renormalization of the many-body
T-matrix, see Eq. (4).

The free-energy calculated in the NSR scheme is equivalent to a model
describing a non-interacting mixture of bosons and fermions (see Appendix D).
The Bose-Fermi model is an accurate description of the thermodynamics of a
normal state p-wave Fermi gas in 3D for the entire BEC-BCS crossover [31].
Here we restrict ourselves to the high-temperature limit which is defined as:
$E_{b}\ll T\ll R^{2}$. The results of our calculations are shown in Fig. 5. As
one can see for the BEC side, the Prandtl number is non-monotonic, but
ultimately decreases as temperature decreases or density increases to about
$25\%$ of the resonant value. At fixed fugacity, the Prandtl number also
decreases slightly as the interaction parameters are decreased.

Figure 5: Prandtl number on the BEC side as a function of fugacity, $z$, and
density through the Fermi temperature $T_{F}$, for various values of the
binding energy, $\beta E_{b}$. For the top (bottom) figures, the lines from
bottom to top denote towards (away from) resonance. In general, decreasing the
interaction decreases the Prandtl number, while lowering the temperature
(increasing the fugacity or increasing the density) tends to an ultimate
decrease of the Prandtl number by about $25\%$.

##  VII Conclusions

In this article we have examined the shear viscosity and thermal conductivity
of a 3D p-wave spin-polarized Fermi gas. We found that the scattering times
are proportional to $R^{2}$ near resonance and $v^{-2}$ for weakly interacting
systems. This means that the transport properties explicitly depend on the
scattering parameters for arbitrary interaction strengths. However, the
Prandtl number in the weakly and strongly interacting limits is indeed a
function of only the equation of state, like the Prandtl number for spin-1/2
s-wave Fermi gases. Unlike the spin-1/2 s-wave Fermi gas, there is no scale
symmetry restricting the Prandtl number to be this way.

Our analysis is valid at high temperatures in the normal phase. It quantifies
the role of two-body correlations in the transport properties of spin-
polarized Fermi gases in 3D. As one goes to lower temperatures, it is
important to also account for three-body correlations and losses. Such losses
have been previously examined [7, 8, 9, 10, 18, 19, 21, 22, 23, 24], and were
found to become quite strong at the p-wave Feshbach resonance. These three-
body losses are suppressed in the high-temperature limit as they are of order
$z^{3}$. Obviously these losses will have an important contribution to the
transport at lower temperatures, but for temperatures $T\gtrsim T_{F}$ the
transport properties are determined by the two-body scattering processes
discussed in this work. We predict a pronounced dip in the viscosity at
temperatures above $T_{F}$ that could be observed in experiment.

##  VIII Acknowledgements

This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation), project-ID 273811115 (SFB1225 ISOQUANT) and under
Germany’s Excellence Strategy EXC2181/1-390900948 (the Heidelberg STRUCTURES
Excellence Cluster), as well as HK GRF 17304820, 17304719, 17305218 and CRF
C6009-20G and C7012-21G. The authors would like to thank Shizhong Zhang for
useful discussions.

##  Appendix A Many-Body Contribution to the Two-Body T-matrix

The many-body contribution to the two-body T-matrix has the form:

$\displaystyle A_{\text{mb}}(Q,$ $\displaystyle E)=24\pi$
$\displaystyle\times\int\frac{d^{3}{\bf
k}}{(2\pi)^{3}}k^{2}\frac{1-u^{2}}{2}\frac{n_{F}(\xi_{k})}{E-\frac{1}{2}Q^{2}-k^{2}-Qku+2\mu}$

where $u=\hat{Q}\cdot\hat{k}$, $\xi_{k}=k^{2}/2-\mu$, and $n_{F}(x)=(e^{\beta
x}+1)^{-1}$ is the Fermi-Dirac distribution function. In writing Eq.
(LABEL:app:Amb_1), we have neglected terms that produce a contribution to the
T-matrix of the form: $\left({\bf p\cdot Q}\right)\left({\bf Q\cdot
q}\right)\tilde{T}(E)$ as these vanish when the ultraviolet cutoff is taken to
infinity.

The angular integration in Eq. (LABEL:app:Amb_1) can be performed analytically
to give:

$\displaystyle A_{\text{mb}}(Q,E)$
$\displaystyle=\int_{0}^{\infty}\frac{dk}{2\pi^{2}}\
k^{4}\frac{e^{-\frac{\beta}{2}k^{2}}z}{1+e^{-\frac{\beta}{2}k^{2}}z}$
$\displaystyle\left[\frac{E+2\mu-k^{2}-Q^{2}/2}{2Q^{2}k^{2}}\right.$
$\displaystyle\left.+\frac{(E+2\mu-k^{2}-Q^{2}/2)^{2}-Q^{2}k^{2}}{4Q^{3}k^{3}}\right.$
$\displaystyle\left.\times\ln\left[\frac{E+2\mu-\frac{1}{2}k^{2}-\frac{1}{2}(Q+k)^{2}}{E+2\mu-\frac{1}{2}k^{2}-\frac{1}{2}(Q-k)^{2}}\right]\right].$
(A.2)

For our purposes we evaluate Eq. (A.2) numerically.

##  Appendix B Derivation of the Shear and Thermal Matrices in Eq. (18)

In this appendix we derive the relevant expressions for the matrices
$\mathcal{A}^{(\eta,\kappa)}$ in Eq. (18), using two basis modes. For both the
shear viscosity and thermal conductivity the matrices
$\mathcal{A}^{(\eta,\kappa)}$ have the form

$\displaystyle\mathcal{A}_{i,j}$
$\displaystyle=\frac{\beta}{\left(U_{1},U_{1}\right)}\int_{0}^{\infty}\frac{dQ}{2\pi^{2}}\int_{0}^{\infty}\frac{dp}{2\pi^{2}}\int\frac{d\Omega_{Q}}{4\pi}\int\frac{d\Omega_{p}}{4\pi}\int\frac{d\Omega_{q}}{4\pi}$
$\displaystyle\frac{Q^{2}p^{7}}{4\pi}\left(\hat{\bf p}\cdot\hat{\bf
q}\right)^{2}|T({\bf Q,p})|^{2}$
$\displaystyle\times\frac{1}{2}\frac{1}{\cosh(a)+\cosh(b\hat{\bf
Q}\cdot\hat{\bf p})}\frac{1}{2}\frac{1}{\cosh(a)+\cosh(b\hat{\bf
Q}\cdot\hat{\bf q})}$
$\displaystyle\times\left[U_{i}({\bf\frac{Q}{2}+p})+U_{i}({\bf\frac{Q}{2}-p})\right]$
$\displaystyle\left[U_{j}({\bf\frac{Q}{2}+p})+U_{j}({\bf\frac{Q}{2}-p})-U_{j}({\bf\frac{Q}{2}+q})-U_{j}({\bf\frac{Q}{2}-q})\right]$
(B.1)

where we have used the identity

$\displaystyle n_{\bf\frac{Q}{2}+p}$ $\displaystyle
n_{\bf\frac{Q}{2}-p}(1-n_{\bf\frac{Q}{2}+q})(1-n_{\bf\frac{Q}{2}-q})=$
$\displaystyle\frac{1}{4}\frac{1}{\cosh\left(a\right)+\cosh\left(b\hat{\bf
Q}\cdot\hat{\bf p}\right)}\frac{1}{\cosh\left(a\right)+\cosh\left(b\hat{\bf
Q}\cdot\hat{\bf q}\right)}$ (B.2)

as well as $a=\beta(Q^{2}/8+p^{2}/2-\mu)$ and $b=\beta Qp/2$. We have also
explicitly suppressed the $\eta$ and $\kappa$ indices in Eq. (B.1) as it is
valid for both the shear viscosity and thermal conductivity.

Let us begin by considering the shear viscosity. The two basis modes we will
consider are:

$\displaystyle U_{1}^{\eta}$ $\displaystyle=p_{x}p_{y},$ $\displaystyle
U_{2}^{\eta}$ $\displaystyle=p_{x}p_{y}\beta\left(p^{2}-w_{\eta}\right).$
(B.3)

The constant $w_{\eta}$ can be determined from the Gram-Schmidt method, and is
such that the two modes are orthogonal with inner product (19). Given these
two modes one can then evaluate Eq. (B.1).

However, special care is needed in performing the angular integrations. As an
example consider the case $i=1$ and $j=1$. In this case Eq. (B.1) becomes

$\displaystyle\mathcal{A}^{\eta}_{1,1}$
$\displaystyle=\frac{\beta}{\left(U_{1},U_{1}\right)}\int_{0}^{\infty}\frac{dQ}{2\pi^{2}}\int_{0}^{\infty}\frac{dp}{2\pi^{2}}\int\frac{d\Omega_{\bf
Q}}{4\pi}$ $\displaystyle\sum_{\mu,\nu=x,y,z}\frac{Q^{2}p^{11}}{\pi}|T({\bf
Q,p})|^{2}\left[F_{\eta,2}^{\mu,\nu}F_{\eta,0}^{\mu,\nu}-F_{\eta,1}^{\mu,\nu}F_{\eta,1}^{\mu,\nu}\right].$
(B.4)

The functions $F_{\eta,n}^{\mu,\nu}$ are the angular averages

$\displaystyle F_{\eta,n}^{\mu,\nu}$
$\displaystyle=\left\langle\frac{1}{2}\frac{1}{\cosh(a)+\cosh(b\cos{\theta_{\bf
p}})}\hat{\bf p}_{\mu}\hat{\bf p}_{\nu}\left(\hat{\bf p}_{x}\hat{\bf
p}_{y}\right)^{n}\right\rangle_{\Omega_{\bf p}}$ (B.5)

where $a=\beta(Q^{2}/8+p^{2}/2-\mu)$ and $b=\beta Qp/2$, while
$\langle\cdot\rangle_{\Omega_{\bf p}}$ denotes angular averages over $\hat{\bf
p}$. As shown in Appendix C, the angular averages can be decomposed into a sum
of terms depending on different components of the center of mass unit vector,
$\hat{\bf Q}$, and the integrals:

$I_{n}(Q,p)=\int_{-1}^{1}\frac{dx}{4}\frac{P_{n}(x)}{\cosh(a)+\cosh(bx)}$
(B.6)

where $P_{n}(x)$ is the Legendre polynomial.

The dominant contribution to Eq. (B.5) comes from the term proportional to
$I_{0}(Q,p)$, which is equivalent to making the substitution:

$\displaystyle F_{\eta,n}^{\mu,\nu}\approx I_{0}(Q,p)\left\langle\hat{\bf
p}_{\mu}\hat{\bf p}_{\nu}\left(\hat{\bf p}_{x}\hat{\bf
p}_{y}\right)^{n}\right\rangle_{\Omega_{\bf p}}.$ (B.7)

This approximation becomes exact in the high-temperature limit where Eq. (B.2)
is angular independent and can be pulled out of the angular average, such that
only the term proportional to $I_{0}$ survives. To test the accuracy of this
approximation at lower temperatures we plotted the full angular average
contained in Eq. (B.4) (solid black line) versus the result after the
approximation in Eq. (B.7) (red-dashed line) for both fixed $a$ and $b$ in
Fig. 6. In the non-degenerate regime, the dominant contribution to
$F_{n}^{\mu,\nu}$ comes from $a,b\lesssim 1$. The approximation of retaining
only the terms proportional to $I_{0}$ is exceptionally good and simplifies
the calculation enormously especially in the required parameter regime, cf.
also Ref. [49].

Figure 6: Comparison of the full angular average of Eq. (B.4) (black line) and
the solution from Eq. (B.7) (red-dashed line) for various values of
$a=\beta(Q^{2}/8+p^{2}/2-\mu)$ and $b=\beta Qp/2$. The approximation is quite
successful at capturing the full structure of the angular average for the
desired parameter regime, $a,b\lesssim 1$.

With these approximations the $\mathcal{A}^{\eta}_{1,1}$ becomes

$\displaystyle\mathcal{A}^{\eta}_{1,1}$
$\displaystyle=\frac{\beta}{\left(U_{1}^{\eta},U_{1}^{\eta}\right)}\int\frac{dQ}{2\pi^{2}}\int\frac{dp}{2\pi^{2}}\frac{Q^{2}p^{11}}{75\pi}\left|T({\bf
Q,p})\right|^{2}I_{0}^{2}(Q,p).$ (B.8)

One can then repeat the following analysis for the other elements of the
matrix $A^{\eta}$. The final result is:

$\displaystyle\mathcal{A}^{\eta}=\frac{\beta}{\left(U_{1}^{\eta},U_{1}^{\eta}\right)}$
$\displaystyle\times\int_{0}^{\infty}\frac{dQ}{2\pi^{2}}\int_{0}^{\infty}\frac{dp}{2\pi^{2}}\frac{Q^{2}p^{11}}{75\pi}\left|T({\bf
Q,p})\right|^{2}I_{0}^{2}(Q,p)$
$\displaystyle\begin{bmatrix}1&\beta\left(\frac{7Q^{2}}{12}+p^{2}-w_{\eta}\right)\\\
\beta\left(\frac{7Q^{2}}{12}+p^{2}-w_{\eta}\right)&\beta^{2}\left(\left(\frac{7Q^{2}}{12}+p^{2}-w_{\eta}\right)^{2}+\frac{7Q^{4}}{90}\right)\end{bmatrix}.$
(B.9)

An identical analysis can be done for the thermal conductivity. The two basis
modes are given by

$\displaystyle U_{1}^{\kappa}$
$\displaystyle=p_{x}\left(\frac{p^{2}}{2}-w\right),$ $\displaystyle
U_{2}^{\kappa}$
$\displaystyle=p_{x}\left(\frac{p^{2}}{2}-w\right)\beta\left(\frac{p^{2}}{2}-w_{\kappa}\right),$
(B.10)

where $w$ is the enthalpy per particle and $w_{\kappa}$ is a constant that
ensures that the two basis modes are orthogonal, similar to the case for shear
viscosity.

We again assume that the angular averages can be approximated in the manner of
Eq. (B.7). The final result for the matrix $\mathcal{A}^{\kappa}$ is:

$\displaystyle\mathcal{A}^{\kappa}=\frac{\beta}{\left(U_{1}^{\kappa},U_{1}^{\kappa}\right)}$
$\displaystyle\times\int_{0}^{\infty}\frac{dQ}{2\pi^{2}}\int_{0}^{\infty}\frac{dp}{2\pi^{2}}\frac{Q^{4}p^{11}}{270\pi}\left|T({\bf
Q,p})\right|^{2}I_{0}(Q,p)^{2}$
$\displaystyle\begin{bmatrix}1&\beta(\frac{7Q^{2}}{20}+p^{2}-w-w_{\kappa})\\\
\beta(\frac{7Q^{2}}{20}+p^{2}-w-w_{\kappa})&\beta^{2}\bigl{(}(\frac{7Q^{2}}{20}+p^{2}-w-w_{\kappa})^{2}+\frac{3Q^{4}}{200}\bigr{)}\end{bmatrix}$
(B.11)

In the high-temperature limit, one can analytically perform the integration
over $Q$ to obtain:

$\displaystyle\mathcal{A}^{\eta}=\frac{2^{3/2}z}{75\pi}$
$\displaystyle\int_{0}^{\infty}\frac{d\epsilon}{(2\pi)^{2}}e^{-\epsilon}\epsilon^{5}|T(0,\sqrt{\epsilon})|^{2}$
$\displaystyle\times\begin{bmatrix}1&\epsilon-\frac{7}{2}\\\
\epsilon-\frac{7}{2}&\left(\epsilon-\frac{7}{2}\right)^{2}+\frac{77}{6}\end{bmatrix},$
(B.12) $\displaystyle\mathcal{A}^{\kappa}=\frac{2}{3}\frac{2^{3/2}z}{75\pi}$
$\displaystyle\int_{0}^{\infty}\frac{d\epsilon}{(2\pi)^{2}}e^{-\epsilon}\epsilon^{5}|T(0,\sqrt{\epsilon})|^{2}$
$\displaystyle\times\begin{bmatrix}1&\epsilon-\frac{7}{2}\\\
\epsilon-\frac{7}{2}&\left(\epsilon-\frac{7}{2}\right)^{2}+7\end{bmatrix}.$
(B.13)

There are several important properties of the matrices
$\mathcal{A}^{\eta,\kappa}$. Consider the relaxation time approximation (RTA).
In the RTA we only retain the first mode in $\mathcal{A}^{\eta,\kappa}$. This
approximation gives the following expression for the scattering times:

$\displaystyle\tau_{\eta}^{(1)}$
$\displaystyle=\frac{\beta}{\mathcal{A}^{\eta}_{1,1}},$
$\displaystyle\tau_{\kappa}^{(1)}$
$\displaystyle=\frac{\beta}{\mathcal{A}^{\kappa}_{1,1}}.$ (B.14)

Beyond the relaxation time approximation (BRTA), we keep both basis modes, and
evaluate the whole matrix, $\mathcal{A}^{\kappa,\eta}$, in order to determine
the scattering times. In general one can write the scattering times in the
BRTA approximation in terms of the RTA approximation:

$\displaystyle\frac{\tau^{(2)}_{\eta}}{\tau^{(1)}_{\eta}}$
$\displaystyle=\frac{1}{1-A}$
$\displaystyle\frac{\tau^{(2)}_{\kappa}}{\tau^{(1)}_{\kappa}}$
$\displaystyle=\frac{1}{1-B}$ (B.15)

where $\tau_{\eta,\kappa}^{(2)}$ are the scattering times in the BRTA
approximation, and the constants, $A$ and $B$, are defined as:

$\displaystyle A$
$\displaystyle=\frac{\left(\mathcal{A}^{\eta}_{1,2}\right)^{2}}{\mathcal{A}^{\eta}_{1,1}\mathcal{A}^{\eta}_{2,2}},$
$\displaystyle B$
$\displaystyle=\frac{\left(\mathcal{A}^{\kappa}_{1,2}\right)^{2}}{\mathcal{A}^{\kappa}_{1,1}\mathcal{A}^{\kappa}_{2,2}}.$
(B.16)

In general the constants $A$ and $B$ are positive because the eigenvalues of
$\mathcal{A}^{\eta,\kappa}$ express the (positive) decay rates of different
perturbations. Hence the BRTA prediction for the scattering times will be
larger than the RTA prediction; this is in line with the variational
formulation that each truncated basis provides a lower bound on the true
transport time $\tau$ [45, 51]. In Fig. 7 we show the scattering rates
calculated in both the RTA and BRTA including finite-temperature effects. At
high temperatures the results for the BRTA are very close to the RTA, making
the RTA a good approximation for the scattering times. As one lowers the
temperature, i.e., increases $z$, the difference between the RTA and BRTA
becomes more pronounced.

In the high-temperature limit we find $B>A$ from the specific structure of the
$\mathcal{A}^{\eta,\kappa}$ matrices in Eqs. (B.12)–(B.12), which implies that
$\tau^{(2)}_{\kappa}/\tau^{(1)}_{\kappa}>\tau^{(2)}_{\eta}/\tau^{(1)}_{\eta}$.
This translates into the following inequality for the Prandtl number
calculated in the BRTA, $\Pr^{(2)}$:

$\Pr^{(2)}=\frac{\tau_{\eta}^{(2)}}{\tau_{\kappa}^{(2)}}=\Pr^{(1)}\frac{1-B}{1-A}\leq\Pr^{(1)},$
(B.17)

where $\Pr^{(1)}$ is the Prandtl number calculated in the RTA.

The BRTA matches the RTA, $\Pr^{(2)}=\Pr^{(1)}$, when $A=B=0$. In general this
occurs when the off-diagonal matrix elements vanish,
$\mathcal{A}^{\eta,\kappa}_{1,2}=0$. In the high-temperature limit,
$\mathcal{A}^{\eta,\kappa}_{1,2}=0$ when

$0=\int_{0}^{\infty}d\epsilon\frac{\epsilon^{5}\left(\epsilon-\frac{7}{2}\right)e^{-\epsilon}}{\beta
R^{2}(\epsilon-\beta E_{b})^{2}+\epsilon^{3}}.$ (B.18)

Given $\beta R^{2}\gg 1$, the approximate solutions for Eq. (B.18) are $\beta
E_{b}\approx 0,7/2$. Numerical evaluation of Eq. (B.18) gives $\beta
E_{b}\approx 0.186,3.447$, which is consistent with the red dashed lines in
Fig. 2.

Besides these points where $\Pr^{(2)}=\Pr^{(1)}$, it is important to note the
behaviour of the Prandtl number in the strongly interacting (“res.”) limit
$\beta^{3/2}v^{-1}=0$ and in the weakly interacting (“w.i.”) limit
$\beta^{3/2}v^{-1}\gg 1$, respectively. A direct evaluation of the scattering
times in these two limits yields

$\displaystyle\left.\tau_{\eta}^{(1)}\right|_{\text{res.}}$
$\displaystyle=\frac{\pi}{\sqrt{2}Tz}\frac{25}{576}\frac{R^{2}}{mT},$
$\displaystyle\left.\tau_{\kappa}^{(1)}\right|_{\text{res.}}$
$\displaystyle=\frac{3}{2}\left.\tau_{\eta}^{(1)}\right|_{\text{res.}}$
$\displaystyle\left.\tau_{\eta}^{(1)}\right|_{\text{w.i.}}$
$\displaystyle=\frac{\pi}{\sqrt{2}Tz}\frac{5}{2304}\frac{1}{v^{2}(mT)^{3}},$
$\displaystyle\left.\tau_{\kappa}^{(1)}\right|_{\text{w.i.}}$
$\displaystyle=\frac{3}{2}\left.\tau_{\eta}^{(1)}\right|_{\text{w.i.}}$

within the RTA, while in the BRTA

$\displaystyle\left.\tau_{\eta}^{(2)}\right|_{\text{res.}}$
$\displaystyle=\left(1+\frac{3}{202}\right)\left.\tau_{\eta}^{(1)}\right|_{\text{res.}}$
$\displaystyle\left.\tau_{\kappa}^{(2)}\right|_{\text{res.}}$
$\displaystyle=\left(\frac{3}{2}+\frac{3}{88}\right)\left.\tau_{\eta}^{(1)}\right|_{\text{res.}}$
$\displaystyle\left.\tau_{\eta}^{(2)}\right|_{\text{w.i.}}$
$\displaystyle=\left(1+\frac{75}{226}\right)\left.\tau_{\eta}^{(1)}\right|_{\text{w.i.}}$
$\displaystyle\left.\tau_{\kappa}^{(2)}\right|_{\text{w.i.}}$
$\displaystyle=\left(\frac{3}{2}+\frac{75}{104}\right)\left.\tau_{\eta}^{(1)}\right|_{\text{w.i.}}.$

This yields the scattering ratios in the strongly and weakly interacting
limits that are quoted in the main text in Eq. (23).

Figure 7: Comparing the predictions of the scattering rates,
$\tau_{\eta,\kappa}$ from the relaxation time approximation (RTA), beyond the
relaxation time approximation (BRTA), and the high temperature approximation
for the BRTA. For this simulation, $\beta E_{b}=0$ and $\beta^{1/2}R=100$. The
scattering times in the BRTA are larger than the RTA values, which are smaller
than the high-temperature results. This is consistent with the basic structure
of the BRTA.

##  Appendix C Angular Averages

We report some general formulas for angular averages of functions
$f(x={\bf\hat{Q}\cdot\hat{p}})$ over a variable number of unit vectors,
$\hat{\bf p}$. The calculations are straight forward, but become rather
tedious. We restrict ourselves to an even number of unit vectors:

$\displaystyle\langle f(x)\rangle_{\hat{\bf p}}=\int\frac{d\Omega_{\bf
p}}{4\pi}f(x)=\int_{-1}^{1}\frac{dx}{2}f(x),$ (C.1)
$\displaystyle\langle\hat{\bf p}_{i}\hat{\bf p}_{j}f(x)\rangle_{\hat{\bf p}}$
$\displaystyle=\int\frac{d\Omega_{\bf p}}{4\pi}\hat{\bf p}_{i}\hat{\bf
p}_{j}f(x)$
$\displaystyle=\int_{-1}^{1}\frac{dx}{2}f(x)\left[\frac{1-x^{2}}{2}\delta_{i,j}+\hat{\bf
Q}_{i}\hat{\bf Q}_{j}P_{2}(x)\right],$ (C.2) $\displaystyle\langle\hat{\bf
p}_{i}$ $\displaystyle\hat{\bf p}_{j}\hat{\bf p}_{k}\hat{\bf
p}_{l}f(x)\rangle_{\hat{\bf p}}=\int\frac{d\Omega_{\bf p}}{4\pi}\hat{\bf
p}_{i}\hat{\bf p}_{j}\hat{\bf p}_{k}\hat{\bf p}_{l}f(x)=$
$\displaystyle\int_{-1}^{1}\frac{dx}{2}f(x)\left[\frac{(1-x^{2})^{2}}{8}\left(\delta_{i,j}\delta_{k,l}+\text{2
permutations}\right)\right.$
$\displaystyle+\frac{-5x^{4}+6x^{2}-1}{8}\left(\hat{\bf Q}_{i}\hat{\bf
Q}_{j}\delta_{k,l}+\delta_{i,j}\hat{\bf Q}_{k}\hat{\bf Q}_{l}+\text{2
perms.}\right)$ $\displaystyle\left.+\hat{\bf Q}_{i}\hat{\bf Q}_{j}\hat{\bf
Q}_{k}\hat{\bf Q}_{l}P_{4}(x)\right],$ (C.3) $\displaystyle\langle\hat{\bf
p}_{i}$ $\displaystyle\hat{\bf p}_{j}\hat{\bf p}_{k}\hat{\bf p}_{l}\hat{\bf
p}_{r}\hat{\bf p}_{s}f(x)\rangle_{\hat{\bf p}}=\int\frac{d\Omega_{\bf
p}}{4\pi}\hat{\bf p}_{i}\hat{\bf p}_{j}\hat{\bf p}_{k}\hat{\bf p}_{l}\hat{\bf
p}_{r}\hat{\bf p}_{s}f(x)=$
$\displaystyle\int_{-1}^{1}\frac{dx}{2}f(x)\left[\frac{(1-x^{2})^{3}}{48}\left(\delta_{i,j}\delta_{k,l}\delta_{r,s}+\text{14
perms.}\right)\right.$ $\displaystyle+\frac{7x^{6}-15x^{4}+9x^{2}-1}{48}$
$\displaystyle\left(\hat{\bf Q}_{i}\hat{\bf
Q}_{j}\delta_{k,l}\delta_{r,s}+\delta_{i,j}\hat{\bf Q}_{k}\hat{\bf
Q}_{l}\delta_{r,s}+\delta_{i,j}\delta_{k,l}\hat{\bf Q}_{r}\hat{\bf
Q}_{s}\right.$ $\displaystyle\left.+\text{14 perms.}\right)$
$\displaystyle+\frac{-21x^{6}+35x^{4}-15x^{2}+1}{48}$
$\displaystyle\left(\hat{\bf Q}_{i}\hat{\bf Q}_{j}\hat{\bf Q}_{k}\hat{\bf
Q}_{l}\delta_{r,s}+\hat{\bf Q}_{i}\hat{\bf Q}_{j}\delta_{k,l}\hat{\bf
Q}_{r}\hat{\bf Q}_{s}\right.$ $\displaystyle\left.+\delta_{i,j}\hat{\bf
Q}_{k}\hat{\bf Q}_{l}\hat{\bf Q}_{r}\hat{\bf Q}_{s}+\text{14 perms.}\right)$
$\displaystyle\left.+\hat{\bf Q}_{i}\hat{\bf Q}_{j}\hat{\bf Q}_{k}\hat{\bf
Q}_{l}\hat{\bf Q}_{r}\hat{\bf Q}_{s}P_{6}(x)\right].$ (C.4)

In these equations $P_{n}(x)$ is the Legendre polynomial. We note that the
only term proportional to $P_{0}(x)$ is the term that is independent of the
center-of-mass momentum, ${\bf Q}$. Thus if $f(x)$ does not depend on $x$,
only the first lines of Eqs. (C.1-C.4) are nonzero, reducing to the standard
formulae for averages of unit vectors.

##  Appendix D The Bose-Fermi Model

In this appendix we review the thermodynamics of the Bose-Fermi model, which
describes a non-interacting mixture of bosons with mass $2m$ and chemical
potential $2\mu-E_{b}\theta(-E_{b})$, and fermions of mass $m$ and chemical
potential $\mu$. The pressure of the system is then given by:

$\displaystyle P=\frac{1}{(2\pi\beta)^{3/2}\beta}$
$\displaystyle\left[-\operatorname{Li}_{5/2}(-z)\right.$
$\displaystyle\left.+3\times
2^{3/2}\operatorname{Li}_{5/2}\left(z^{2}e^{-\beta
E_{b}\theta(-E_{b})}\right)\right],$ (D.1)

where $\operatorname{Li}_{a}(z)$ is the polylogarithm function. In Eq. (D.1),
the first term is related to the noninteracting Fermi gas, while the second
term describes the noninteracting bosons. The factor of $3$ originates from
the $2\ell+1$ degeneracy of the p-wave coupling, i.e., there are $3$ bound
states corresponding to coupling in the $m_{\ell}=-1,0,1$ channels, where
$m_{\ell}$ is the azimuthal quantum number.

In this model the density is given by

$\displaystyle n=\frac{1}{(2\pi\beta)^{3/2}}$
$\displaystyle\left[-\operatorname{Li}_{3/2}(-z)\right.$
$\displaystyle\left.+6\times
2^{3/2}\operatorname{Li}_{3/2}\left(z^{2}e^{-\beta
E_{b}\theta(-E_{b})}\right)\right].$ (D.2)

From Eqs. (D.1-D.2) one can calculate the specific heat at constant pressure
per unit volume and the density for arbitrary temperature in the normal phase
[31]. In this article we primarily focus on the limit where $|\beta E_{b}|$ is
small, hence we ignore the bound state energy and its correction to the
bosonic chemical potential.

## References

  * [1] A. J. Leggett, Rev. Mod. Phys. 47 331 (1975).
  * [2] G. E. Volovik, J. Exp. Theor. Phys. 67, 1804 (1988) [Zh. Eksp. Teor. Fix. 94, 123 (1988)].
  * [3] N. Read, and D. Green, Phys. Rev. B 61, 10267 (2000).
  * [4] V. Gurarie, and L. Radzihovsky, Ann. Phys. (Amsterdam) 322, 2-119 (2007).
  * [5] S.-J. Jiang, and F. Zhou, Phys. Rev. A 97, 063606 (2018).
  * [6] F. Yang, S.-J. Jiang, and F. Zhou, Phys. Rev. Lett. 124, 225701 (2020).
  * [7] H. Suno, B. D. Esry, and C. H. Greene, Phys. Rev. Lett. 90, 0533202 (2003).
  * [8] M. Jona-Lasinio, L. Pricoupenko, and Y. Castin, Phys. Rev. A 77, 043611 (2008).
  * [9] M. Schmidt, H.-W. Hammer, L. Platter. Phys. Rev. A 101, 062702 (2020).
  * [10] S. Zhu, Z. Yu, and S. Zhang, arXiv:2201.00962 (2022).
  * [11] X. Cui, Phys. Rev. A 94, 043636 (2016); X. Cui, and H. Dong, Phys. Rev. A 94, 063650 (2016).
  * [12] Y.-C. Zhang, and S. Zhang, Phys. Rev. A 95, 023603 (2017).
  * [13] L. Pan, S. Chen, and X. Cui, Phys. Rev. A 98, 011603 (2018).
  * [14] X. Yin, X.-W. Guan, Y. Zhang, H. Su, and S. Zhang, Phys. Rev. A 98, 023605 (2018).
  * [15] Y. Sekino, S. Tan, and Y. Nishida, Phys. Rev. A 97, 013621 (2018).
  * [16] J. Maki, Phys. Rev. A 104, 063314 (2021).
  * [17] Y. Sekino, and Y. Nishida, Phys. Rev. A 103, 043307 (2021).
  * [18] C. A. Regal, C. Ticknor, J. L. Bohm, and D. S. Jin, Phys. Rev. Lett. 90, 053201 (2003).
  * [19] J. Zhang, E. G. M. van Kempen, T. Bourdel, L. Khaykovich, J. Cubizolles, F. Chevy, M. Teichmann, L. Tarruell, S. J. J. M. F. Kokkelmans, and C. Salomon, Phys. Rev. A 70, 030702(R) (2004).
  * [20] C. H. Schunck, M. W. Zwierlein, C. A. Stan, S. M. F. Raupach, W. Ketterle, A. Simoni, E. Tiesinga, C. J. Williams, and P. S. Julienne, Phys. Rev. A 71, 045601 (2005).
  * [21] F. Chevy, E. G. M. van Kempen, T. Bourdel, J. Zhang, L. Khaykovich, M. Teichmann, L. Tarruell, S. J. J. M. F. Kokkelmans, and C. Salomon, Phys. Rev. A 71, 062710 (2005).
  * [22] J. P. Gaebler, J. T. Stewart, J. L. Bohn, and D. S. Jin, Phys. Rev. Lett. 98, 200403 (2007).
  * [23] J. Fuchs, C. Ticknor, P. Dyke, G. Veeravalli, E. Kuhnle, W. Rowlands, P. Hannaford, and C. J. Vale, Phys. Rev. A 77, 053616 (2008).
  * [24] Y. Inada, M. Horikoshi, S. Nakajima, M. Kuwata-Gonokami, M. Ueda, and T. Mukaiyama, Phys. Rev. Lett. 101, 100401 (2008).
  * [25] J. Yoshida, T. Saito, M. Waseem, K. Hattori, and T. Mukaiyama, Phys. Rev. Lett. 120, 133401 (2018).
  * [26] Y.-T. Chang, R. Senaratne, D. Cavazos-Cavazos, and R. G. Hulet, Phys. Rev. Lett. 125, 263402 (2020).
  * [27] A. S. Marcum, F. R. Fonta, A. M. Ismail, and K. M. O’Hara, arXiv:2007.15783 (2020).
  * [28] D. J. M. Ahmed-Braun, K. G. Jackson, S. Smale, C. J. Dale, B. A. Olsen, S. J. J. M. F. Kokkelmans, P. S. Julienne, and J. H. Thywissen, Phys. Rev. Research 3, 033269 (2021).
  * [29] K. G. Jackson, C. J. Dale, J. Maki, K. G. S. Xie, B. A. Olsen, D. J. M. Ahmed-Braun, S. Zhang, and J. H. Thywissen, arXiv:2206.10415 (2022).
  * [30] Z. Yu, J. H. Thywissen, and S. Zhang, Phys. Rev. Lett. 115, 135304 (2015).
  * [31] J. Yao, and S. Zhang, Phys. Rev. A 97, 043612 (2018).
  * [32] S. Ding, and S. Zhang, Phys. Rev. Lett. 123, 070404 (2019).
  * [33] A. Schwenk, and C. J. Pethick, Phys. Rev. Lett. 95, 160401 (2005).
  * [34] D. Lacroix, Phys. Rev. A 94, 043614 (2016).
  * [35] L. M. Schonenberg, and G. J. Conduit, Phys. Rev. A 95, 013633 (2017).
  * [36] G. A. Miller, Phys. Lett. B 777, 442 (2018).
  * [37] J. Maki and S. Zhang, Phys. Rev. Lett. 125, 240402 (2020).
  * [38] K. Dusling, and T. Schäfer, Phys. Rev. Lett. 111, 120603 (2013).
  * [39] Y. Nishida, Ann. of Phys. 410, 167949 (2019).
  * [40] T. Enss, Phys. Rev. Lett. 123, 205301 (2019).
  * [41] J. Hofmann, Phys. Rev. A 101, 013620 (2020).
  * [42] K. Fujii, and Y. Nishida, Phys. Rev. A 102, 023310 (2020).
  * [43] K. Fujii, and T. Enss, arXiv:2208.03353 (2022).
  * [44] G. M. Bruun, and H. Smith, Phys. Rev. A 72, 043605 (2005).
  * [45] G. M. Bruun, and H. Smith, Phys. Rev. A 75, 043612 (2007).
  * [46] G. Rupak, and T. Schäfer, Phys. Rev. A 76, 053607 (2007).
  * [47] M. Braby, J. Chao, and T. Schäfer, Phys. Rev. A 82, 033619 (2010).
  * [48] T. Enss, R. Haussmann, and W. Zwerger, Ann. Phys. (NY) 326, 770 (2011).
  * [49] T. Enss, Phys. Rev. A 86, 013616 (2012).
  * [50] M. Bluhm, J. Hou, and T. Schäfer, Phys. Rev. Lett. 119, 065302 (2017).
  * [51] B. Frank, W. Zwerger, and T. Enss, Phys. Rev. Research 2, 023301 (2020).
  * [52] K. Fujii, and Y. Nishida, Phys. Rev. A 103, 053320 (2021).
  * [53] E. M. Lifshitz, and E. M. Pitaevskii, Physical Kinetics (Butterworth-Heinemann, 1981).
  * [54] H. Smith, and H. H. Jensen, Transport Phenomena (Clarendon Press, 1989).
  * [55] In principle the two-body T-matrix can be a rank-2 tensor of the form: $\langle{\bf Q}/2\pm{\bf p}|\mathcal{T}|{\bf Q}/2\pm{\bf q}\rangle=p_{i}q_{i}T_{i,j}({\bf Q},E)$, where $T_{i,j}({\bf Q},E)=\delta_{i,j}T^{(0)}({\bf Q},E)+Q_{i}Q_{j}T^{(1)}({\bf Q},E)$. However, $T^{(1)}$ vanishes like $\Lambda^{-3}$, where the ultraviolet cutoff $\Lambda$ is sent to infinity.
  * [56] Y. Nishida, Phys. Rev. A 86, 012710 (2012).
  * [57] E. Braaten, P. Hagen, H.-W. Hammer, and L. Platter, Phys. Rev. A 86, 012711 (2012).
  * [58] P. K. Kovtun, D. T. Son, and A. O. Starinets, Phys. Rev. Lett. 94, 111601 (2005).
  * [59] P. Nozières, and S. Schmitt-Rink, J. Low Temp. Phys. 59, 195 (1985).