# Stellar-mass microlensing of gravitational waves

Mark H. Y. Cheung,1 Joseph A. J. Gais,1 Otto A. Hannuksela2,3 and Tjonnie G.
F. Li1
1Department of Physics, The Chinese University of Hong Kong, Shatin, N.T.,
Hong Kong.
2Nikhef – National Institute for Subatomic Physics, Science Park, 1098 XG
Amsterdam, The Netherlands
3Institute for Gravitational and Subatomic Physics (GRASP), Department of
Physics, Utrecht University, Princetonplein 1, 3584 CC Utrecht, The
Netherlands E-mail<EMAIL_ADDRESS>
###### Abstract

When gravitational waves pass through the nuclear star clusters of galactic
lenses, they may be microlensed by the stars. Such microlensing can cause
potentially observable beating patterns on the waveform due to waveform
superposition and magnify the signal. On the one hand, the beating patterns
and magnification could lead to the first detection of a microlensed
gravitational wave. On the other hand, microlensing introduces a systematic
error in strong lensing use-cases, such as localization and cosmography
studies. We show that diffraction effects are important when we consider GWs
in the LIGO frequency band lensed by objects with masses $\lesssim 100\,\rm
M_{\odot}$. We also show that the galaxy hosting the microlenses changes the
lensing configuration qualitatively, so we cannot treat the microlenses as
isolated point mass lenses when strong lensing is involved. We find that for
stellar lenses with masses $\sim 1\,\rm M_{\odot}$, diffraction effects
significantly suppress the microlensing magnification. Thus, our results
suggest that gravitational waves lensed by typical galaxy or galaxy cluster
lenses may offer a relatively clean environment to study the lens system, free
of contamination by stellar lenses. We discuss potential implications for the
strong lensing science case. More complicated microlensing configurations will
require further study.

###### keywords:

gravitational waves – gravitational lensing: strong – gravitational lensing:
micro

††pubyear: 2020††pagerange: Stellar-mass microlensing of gravitational
waves–15

##  1 Introduction

Both gravitational waves (GWs) and electromagnetic (EM) waves can be
gravitationally lensed when propagating near massive astrophysical objects,
such as galaxies and galaxy clusters (Ohanian, 1974; Thorne, 1982; Deguchi &
Watson, 1986; Wang et al., 1996; Nakamura, 1998; Takahashi & Nakamura, 2003).
To date, only the lensing of EM waves have definitively been observed, but
several searches for GW lensing have already been performed (Hannuksela et
al., 2019; Li et al., 2019a; McIsaac et al., 2020; Pang et al., 2020; Dai et
al., 2020; Liu et al., 2020).

Suppose the GW is lensed by a galaxy or a galaxy cluster only. In that case,
we can use the geometrical optics approximation to solve for the trajectory,
the arrival time of distinct GW images, and their image types Takahashi &
Nakamura (2003). These GW images will only differ in their amplitudes, arrival
times, and overall phases.111This generally applies to (2,2) modes. However,
if the GWs exhibit higher-order modes or precession, the so-called saddle
point images (type-II) could experience nontrivial waveform changes (Dai &
Venumadhav, 2017; Ezquiaga et al., 2020). Thus, we can identify the strongly
lensed images as repeated events with identical frequency evolution using
matched filtering and Bayesian analysis (Haris et al., 2018; Hannuksela et
al., 2019; Li et al., 2019a; McIsaac et al., 2020; Dai et al., 2020; Liu et
al., 2020).

Several exciting science cases have been proposed for such strongly lensed
events. For example, suppose we observe four GW images of a single event. In
that case, it will allow us to uniquely match the GW image properties to any
given lensed galaxy, which may enable us to localize the host galaxy and
pinpoint the location of the merger event within the galaxy (Hannuksela et
al., 2020). Two images might still allow us to constrain the number of
candidates (Sereno et al., 2011; Yu et al., 2020), but to a lesser degree as
we will need to rely mainly on the magnification ratios to pinpoint the source
location.222Thus, the image properties are largely degenerate with the lens
and source properties (Hannuksela et al., 2020). In particular, we will need
to rely on magnification ratio measurements, subject to much more uncertainty
than the time delays, to identify the source position. Four images allow us to
map the position uniquely based on the time delays alone. However, if a galaxy
cluster lenses the event, even two events might allow localization, owing to
the rarity of such clusters (Smith et al., 2018; Smith et al., 2017, 2019;
Robertson et al., 2020; Ryczanowski et al., 2020).

If we localized such strongly lensed GWs, we might be able to perform
precision cosmography studies at high redshift (Sereno et al., 2011; Liao et
al., 2017; Cao et al., 2019; Li et al., 2019b; Hannuksela et al., 2020), test
the speed of gravity at high precision (Collett & Bacon, 2017; Baker &
Trodden, 2017; Fan et al., 2017), connect binary black holes to their host
galaxy properties and investigate formation channels (Chen & Holz, 2016;
Hannuksela et al., 2020), and conduct improved tests of the GW polarization
content (Goyal et al., 2020).333We note that these latter polarization tests
could be done even for non-localized events (Goyal et al., 2020). Still, the
localization would significantly improve the tests, as we could solve the sky
location and the image properties independently in the electromagnetic band
(Hannuksela et al., 2020).  Other intriguing use-cases in, e.g., cosmography,
have been documented in (Oguri, 2019). However, most of these science cases,
as well as the localization studies, rely on a clean environment devoid of
systematic errors introduced due to microlensing by stars.444For example, a
substantial fraction of the error budget in the localization studies proposed
in Ref. Hannuksela et al. (2020) is dominated by microlensing.

In particular, when GWs are strongly lensed by galaxy and galaxy clusters,
their trajectories might pass through dense nuclear star clusters, introducing
a non-negligible chance of microlensing (Christian et al., 2018). The stellar
lens would then introduce additional microlensing effects. In gravitational-
wave lensing localization studies, they can introduce a sizable error in the
magnification ratio or even the time-delay measurement when the geometrical
optics limit is assumed (Hannuksela et al., 2020). In general, microlensing
could lead to systematic errors which will be counter-beneficial to studies of
lensed transient events (see discussion in Oguri (2019)) and indeed almost any
study relying on gravitational-wave localization (e.g., Sereno et al. (2011);
Collett & Bacon (2017); Baker & Trodden (2017); Fan et al. (2017); Liao et al.
(2017); Cao et al. (2019); Li et al. (2019b); Hannuksela et al. (2020)).

Luckily, as GWs have much longer wavelengths than typical EM waves,
diffraction effects can suppress microlensing when the stellar lens has an
Einstein radius comparable to the wavelength of the GWs (Oguri, 2019).
However, microlenses in the wave optics limit have been mostly studied in
isolation (see, however, Diego et al. (2019) and Diego (2020) for an
interesting exploration of microlensing effects in the extreme high-
magnification regime).

We consider the typical case of a stellar lens on top of a galactic lens model
at moderate magnification, similar to the usual strong lensing system we
expect to observe. We focus on the waveform suppression effects brought forth
by wave optics diffraction and elaborate on its implications. Much of our work
is motivated by the discussion in Ref. Oguri (2019), which discusses some of
these wave optics suppression effects; our work aims to quantify the wave
optics suppression when macromodels are involved and investigate their
dependence on the microlens properties and position in order to obtain a
qualitative understanding.

We first lay out the procedures for treating a compound lens system consisting
of lenses with masses and length scales orders of magnitudes apart. Then, we
show that the geometrical optics approximation is not sufficient to treat the
case of stellar-mass microlensing for LIGO/Virgo. After that, we quantify the
deviations introduced by microlenses to the GW waveforms and show that they
would generally not affect events at LIGO/Virgo frequencies. Finally, we will
put our results in an astrophysical context and elaborate on the implications
and potential use-cases, especially those related to the localization of GW
sources and cosmography.

##  2 Methods

Strongly lensed GWs may be detectable at design sensitivity, based on current
lens and binary models (Ng et al., 2018; Li et al., 2018; Oguri, 2018). The
typical analysis considers gravitational lensing by a macromodel, such as a
galaxy. However, microlensing by, e.g., a point-mass within a larger galaxy
macromodel may significantly affect the lensed GW waveform (Lai et al., 2018;
Christian et al., 2018; Jung & Shin, 2019; Diego et al., 2019; Diego, 2020).
Suppose the wavelength of the GW $\lambda$ is comparable to the Einstein
radius of the point mass $\theta_{m}$. In that case, there will be significant
diffraction effects, in which case the geometrical optics approximation no
longer holds (Takahashi & Nakamura, 2003). In this section, we review the
necessary methods to treat lensing in the full wave-optics regime. Note: in
this paper, we work in units where $c=G=1$.

###  2.1 Computing the Diffraction Integral

Figure 1 shows the geometry of a typical lensing configuration. A GW source is
positioned at an angular diameter distance of $D_{L}$. In contrast, a massive
object, or ‘lens’, is positioned near the straight-line path between the
source and the observer at an angular diameter distance $D_{S}$. As $D_{L}$
and $D_{S}$ are very large compared to the length scale of the lens, we can
apply the thin-lens approximation, where the lensing of GWs occurs in the
plane of the lens.

Figure 1: Schematic diagram of the lensing geometry. We use the thin-lens
approximation and assume that all of the lens mass resides on the lens plane.
This diagram is not drawn to scale; the lengths $D_{LS}$, $D_{L}$ and $D_{S}$
are very large compared to $\|\vec{\eta}\|$ and $\|\vec{\xi}\|$.

When GWs are radiated from the source, the portion that passes near the lens
might be focused toward the observer. As the lens’s length scale is small
compared to the distances between the source, lens, and observer, we can
approximate the paths of lensed ‘rays’ of GWs as straight line segments that
change their direction abruptly at the lens plane. In Figure 1, we show one of
these rays. If we draw all possible rays from the source that would arrive at
the observer, we would observe a one-to-one correspondence of such rays with
points on the lens plane. Moreover, as these rays take different spacetime
paths, they would arrive at the observer at different times. Thus, we can
construct a function $t_{d}$ that maps a point on the lens plane $\vec{\xi}$
to the arrival time of the ray that passes through that point, given a source
position $\vec{\eta}$. For analytical calculations with $\vec{\xi}$ and
$\vec{\eta}$, it is easier to work in units normalized by the Einstein radius

$\theta_{E}=\sqrt{\frac{4M_{L}D_{L}D_{LS}}{D_{S}}},$ (1)

where $M_{L}$ is the total mass of the lens and $D_{LS}$ is the angular
diameter distance between the lens and the source (in systems with multiple
lenses, we will revert to un-normalized units to avoid confusion). In these
normalized units, writing $\vec{x}=\vec{\xi}/\theta_{E}$ and
$\vec{y}=\vec{\eta}/\theta_{E}$, the time delay function $t_{d}$ is given by

$t_{d}(\vec{x},\vec{y})=\frac{D_{S}\theta_{E}^{2}(1+z_{L})}{D_{L}D_{LS}}T(\vec{x},\vec{y}),$
(2)

with

$T(\vec{x},\vec{y})=\frac{1}{2}\|\vec{x}-\vec{y}\|^{2}-\phi(\vec{x})-T_{0}(\vec{y}),$
(3)

where $T_{0}(\vec{y})$ corresponds to the arrival time of the first image, and
$\phi(\vec{x})$ is the Fermat potential which depends on the nature of the
lens (e.g. we have $\phi(\vec{x})=\ln{\|\vec{x}\|}$ for a point mass lens).
Besides, we consider the time delay on cosmological distances, such that eq. 2
includes a factor of $1+z_{L}$, where $z_{L}$ is the lens’s redshift.

Let the GW waveform observed without the presence of the lens be $h(f)$ in the
frequency domain. If the waveform observed with lensing is $h_{L}(f)$, we can
quantify the change in the waveform due to lensing effects by the
amplification factor

$F(f)=\dfrac{h_{L}(f)}{h(f)},$ (4)

which is complex in general. If $\phi$ and hence $t_{d}$ is known, we can
obtain $F(f)$ by the diffraction integral (Schneider et al., 1992)

$F(f)=\frac{D_{S}\theta_{E}^{2}(1+z_{L})}{D_{L}D_{LS}}\frac{f}{i}\iint_{\mathbb{R}^{2}}d^{2}\vec{x}\
e^{2\pi ift_{d}(\vec{x},\vec{y})}.$ (5)

By using a dimensionless frequency scale $w=8\pi M_{L}(1+z_{L})f$, we arrive
at

$\tilde{F}(w)=-iw\iint_{\mathbb{R}^{2}}d^{2}\vec{x}\
e^{iwT(\vec{x},\vec{y})}.$ (6)

The integrand in eq. 6 is highly oscillatory for high $w$, so the integral is
computationally expensive to solve with typical numerical integration
techniques even for simple lens models. To solve this integral, we follow a
similar method mentioned in Ulmer & Goodman (1994) and used in Diego et al.
(2019), where we compute an area integration directly in the lens plane and
take its Fourier transform, which is easier to implement for general lens
models and gives $F(f)$ for a wide frequency spectrum with a single
calculation. The case of microlensing on a macromodel introduces important
subtleties to this method, making the evaluation of the integral not as
straight forward as models containing lenses with masses all of the same
order. In particular, the regions around the time delay function’s saddle
points require careful treatment (see Sec. 4.2). These cases and the detailed
methodology for calculating eq. 6 for general lens models will be discussed in
a later work (Cheung et al., in preparation). In this paper, we focus on the
results from the minimum time delay image of the macromodel, which can already
give us insights into the effects of stellar lensing on GWs.

###  2.2 Geometrical optics approximation

For $w\gg 1$, we can work in the geometrical optics limit, where the integral
in eq. 6 will mostly be contributed by the stationary points of
$T(\vec{x},\vec{y})$ (Fermat’s principle). We call these stationary points the
image positions of the lensing configuration. In such a limit, $\tilde{F}(w)$
is given by (Schneider et al., 1992)

$\tilde{F}(w)=\sum_{j}|\mu(\vec{x_{j}})|^{1/2}\exp{(i\pi(2wT(\vec{x_{j}},\vec{y})-n_{j}))},$
(7)

where $\vec{x_{j}}$ are the image positions,
$\mu(\vec{x})=\det(\partial\vec{y}/\partial\vec{x})$ is the image
magnification, $n_{j}=0,\frac{1}{2},1$ when $\vec{x_{j}}$ is at a minimum,
saddle, and maximum point of $T(\vec{x},\vec{y})$ respectively. We solve the
image properties in the geometrical optics limit using the LensingGW package
(Pagano et al., 2020).

###  2.3 Microlensing

Figure 2: Schematic diagram of an example scenario of microlensing on a
macromodel. (Top) In the geometrical optics limit, one can treat the waves as
rays propagating from the source to the observer, only taking specific routes
through particular points on the lens plane. If there were no microlens, there
would be two of these rays, passing through macroimage points (blue points) on
the lens plane (bottom panel) before arriving at the observer. However, when a
microlens (small black dot on the lens plane) is present near a macroimage,
the beam (shown in dotted lines) might be split into multiple ones and pass
through different points on the lens plane, which correspond to microimage
positions (red points). The distance between the microlens and the macroimage
is written as $\eta_{m}$ and called the micro impact parameter.

Strong lensing occurs when a GW passes by a lens with an impact parameter
$\lesssim$ the Einstein radius. A scenario particularly interesting would be
when a GW passes inside the Einstein ring of two lenses of drastically
different lens scales. It is useful for such cases to introduce the notion of
‘macro’ and ‘micro’ lens models. We refer to the lensing images produced by
the two models as macroimages and microimages, respectively, following Diego
et al. (2019).

Let there be some lenses of mass $\sim M$ and some of mass $\sim m$ in the
lens plane (where the thin-lens approximation is invoked), and let $M\gg m$.
For simplicity, let there be no relevant lenses in the intermediate-mass
range. We can then define the $\sim M$ lenses to be the macrolenses, and those
$\sim m$ to be the microlenses. We can obtain an intuitive understanding of
the behavior of GWs lensed by such a system by thinking of how the lenses
would split the GW beams. We can think of the macrolenses splitting the beams
into separate rays (macroimages), and then the microlenses split the beams
further (into microimages). An example is shown in Figure 2.

However, there is a caveat: when thinking of lensed GWs as propagating towards
us by taking certain discrete ray-like paths, we have implicitly invoked the
geometrical optics approximation. As mentioned previously in this section, if
the wavelength of the GWs is comparable to the Einstein radius of the lenses,
diffraction effects will occur, and the GWs will no longer take only a
discrete number of paths. The next section will show that the geometrical
optics approximation is insufficient for treating microlenses of stellar mass
for GWs detectable by LIGO. Nonetheless, the microimage positions and the
macroimage positions are always the stationary points of the time delay
function $t_{d}$, and it is convenient to keep track of their positions even
if we are working with a full wave-optics treatment.

The Einstein radii of the microlenses are orders of magnitudes smaller than
those of the macrolenses. The difference in arrival time between a group of
microimages near a particular macroimage is also orders of magnitude below the
difference in arrival time between macroimages. Therefore, we can carry out
our integration in Eq. 6 for the microimages near the minimum macroimage of
the macromodel (galactic lens), that is, within a small domain on the lens
plane near the macroimage. We will then obtain an amplification factor $F(f)$
that corresponds to the interferences between the relevant microimages with
diffraction included. For saddle point macroimages, we will have to cover a
much more extensive integration domain. This is to be performed in a later
work (see also Sec. 4.2).

##  3 Results

This section quantifies the effects of stellar microlenses on GWs (within the
LIGO frequency band) lensed by a galaxy. We place the source at a redshift of
$z_{S}=1$ and place a galactic lens at $z_{L}=0.5$. For our galactic
macromodel, we use a singular isothermal sphere (SIS) lens of mass
$10^{10}\,\rm M_{\odot}$. We will only consider cases corresponding to strong
lensing. The macro impact parameter $y$ (position of the source in the sky
with respect to the center of the SIS macromodel) is less than
$\theta_{\mathit{SIS}}$, the Einstein radius of the SIS macromodel. For such
cases, we will see two macroimages, one corresponding to a minimum point in
the time delay function, and the other corresponding to the saddle point.
Then, we introduce point-mass microlenses in the vicinity of the minimum
macroimage. We can naively treat the minimum macroimage position as the
‘source position’ for the microlensing configuration and define the micro
impact parameter $\eta_{m}$ as the distance of the microlens from the
macroimage position.

We assume that the microlens is hosted by the macromodel galaxy, which will
likely be true for most cases of microlensing where the microlens is a star
from the nuclear cluster of the galaxy, so we can treat the macrolens and the
microlens as residing on the same plane in the sky (Diego et al., 2019). If
the microlenses are of stellar mass, say $m<100M_{\odot}$, the arrival time of
the microimages will reside within a time window which is less than the period
of GWs measured by LIGO. Thus, we will have to use the methods discussed in
the previous section to recover the amplification factor for the signal
measured that corresponds to the lensed GW rays near the macroimage.

In this section, as we are working with two very different length scales
corresponding to the macromodel and microlenses, we will revert to
unnormalized units of length, i.e., using $\vec{\xi}$ and $\vec{\eta}$ instead
of $\vec{x}$ and $\vec{y}$. We denote the norms of the vectors
$\xi=\|\vec{\xi}\|$ and $\eta=\|\vec{\eta}\|$ etc. Relevant lengths for the
SIS macromodel are the macro impact parameter (distance from the center of the
SIS lens to the source, projected onto the lens plane) $\eta_{\mathit{SIS}}$,
and the Einstein radius of the SIS lens $\theta_{\mathit{SIS}}$. For
microlenses near a parent macroimage, we denote the micro impact parameter
(distance between the parent macroimage and the microlens on the lens plane)
as $\eta_{m}$, and the Einstein radius of the microlens as $\theta_{m}$.

###  3.1 The necessity of a full wave-optics analysis

Figure 3: The lens plane of an SIS macromodel of mass $M=10^{10}\,\rm
M_{\odot}$ with a microlens of mass $m=30\,\rm M_{\odot}$ near the minimum
macroimage. The source position measured from the center of the SIS lens and
projected onto the lens plane (macro impact parameter) is
$\eta_{\mathit{SIS}}=0.1\,\theta_{\mathit{SIS}}$, while the microlens is
placed $\eta_{m}=1.0\,\theta_{m}$ away from the macroimage. The contour lines
of the time delay function $t_{d}$ are shown in grey. (Left panel) view of the
complete model, including the source position and macroimages. The macroimage
on the left is of type-II (negative parity) as it is a saddle point of
$t_{d}$, while that on the right is of type-I (positive parity). (Center
panel) zoom in on the region around the type-I minimum macroimage. The
microlens splits the macroimage into four microimages, two of type-I (both
minimum points) and two of type-II. (Right panel) Further zoom in onto the
macroimage position. Note that as the macroimage is split into microimages,
the point corresponding to the macroimage position is not a stationary point
of $t_{d}$ anymore; it is an actual image only when the microlens is absent.

Figure 3 shows the configuration of the full lens model and the contours of
the time delay function for the case of a macro impact parameter
$\eta_{\mathit{SIS}}=0.1\,\theta_{\mathit{SIS}}$, micro impact parameter
$\eta_{m}=1.0\,\theta_{m}$, and microlens mass $m=30\,\rm M_{\odot}$. The
microlens is placed on the straight line between the SIS center and the
minimum macroimage position. More general cases are studied in Appendix B and
C. The macroimage and microimage positions are found with the help of the
LensingGW package (Pagano et al., 2020).

We can work in the geometrical optics limit for the macromodel (the SIS lens)
due to the lens’s considerable mass. We will observe two GW signals of the
same phase but different amplitudes and arrival times (see the two macroimages
in the left panel of Figure 3). In this case, the macroimage only induces an
overall magnification on the waveform, such that $F(f)$ is a constant.

However, the addition of a stellar-mass lens near the macroimage disturbs the
time delay function near the arrival time of the macroimage, causing $F(f)$ to
deviate from the constant value. Figure 4 shows $F(f)$ for the GW rays that
come from around the minimum macroimage. Note that the amplification factor is
normalized by $\sqrt{\mu}$, where $\mu$ is the magnification of the macroimage
in the geometrical optics limit. If we work with the geometrical optics
approximation and use eq. 7, we obtain the dashed curve, which deviates
significantly from the results without approximations for lower frequencies,
showing that a full wave-optics analysis is necessary when dealing with
microlenses.

We also demonstrate that if we naively treated the microlenses independently
of the macromodel (in isolation), we obtain the dotted curve in Figure 4
instead, corresponding to the standard analytical expression for the
amplification factor of the point mass lens (see Peters (1974)). This
amplification factor is significantly different from the full non-approximate
amplification factor. In fact, by observing that the microlens splits the
macroimage into four microimages instead of two in the center panel of Figure
3, we can already infer that the macromodel would introduce qualitative
differences from the case of an isolated point mass. This further confirms the
results in Diego et al. (2019); Diego (2020); Pagano et al. (2020), that it is
necessary to consider the macromodel when treating microlensing in the
presence of strong lensing.

Figure 4: The normalized norm (top, left axis) and phase (bottom, left axis)
of the frequency domain amplification factor $F(f)$ near the minimum
macroimage for the lensing set up described in Figure 3 ($30\,\rm M_{\odot}$
microlens). The deviation of the geometrical optics approximation (red dashed
line) and the full wave optics treatment (red solid line) demonstrates that
wave optics effects are important for stellar lenses of lower masses, and they
suppress the amplification factor. The disagreement between treating the
microlens in isolation (red dotted line) and considering both the macro and
microlenses (red solid line) shows that we have to also consider the
macromodel’s contribution to the time delay function in our analysis. The PSD
of LIGO is plotted in grey dotted lines with the axis on the right of each
panel.

###  3.2 Waveform suppression due to diffraction

We saw that the normalized amplification factor in our full wave-optics
analysis reduces to unity as $f\rightarrow 0$, which corresponds to the limit
when the wavelength of the GW $\lambda>\theta_{m}$, where the waves would
propagate through the microlens unimpeded as diffraction effects dominate.
This effectively means that GWs will be largely unaffected by the microlenses
in the very low-mass regime. We examine the level of waveform suppression for
different microlens masses in Figure 5. With $m\in\\{1,5,10,30\\}\,\rm
M_{\odot}$, we see that lower microlens masses correspond to amplification
factors approaching $1$ at most relevant LIGO/Virgo frequencies, implying that
the GW signal is amplified only by the macroimage magnification. We also plot
the frequency domain and time domain waveforms for $m\in\\{0,10\\}\,\rm
M_{\odot}$ in Figure 6. Here, an IMRPhenomD waveform of a $30+30\,\rm
M_{\odot}$ binary black hole merger generated by pycbc (Nitz et al., 2020) is
used. The result obtained from the geometrical optics approximation for the
$m=10\,\rm M_{\odot}$ waveform is also shown in dashed lines, showing that we
would overestimate the deviations of the waveform due to microlensing if we
work in the geometrical optics limit.

Figure 5: Same as Figure 4, but with a point mass microlens of mass
$m\in\\{1,5,10,30\\}\,\rm M_{\odot}$, placed accordingly at a distance of
$\eta_{m}=1.0\,\theta_{m}$ away from the macroimage position similar to the
set up in Figure 3. $F/\sqrt{\mu}$ converges to the constant limit of $1$ for
low frequencies and low microlens masses. We show the geometrical optics
approximation results for a $1\,\rm M_{\odot}$ microlens in the green dashed
lines. The wave optics diffraction effects suppress the waveform relative to
the geometrical optics approximation.  Figure 6: The frequency domain (top)
and time domain (bottom) waveforms for the cases $m\in\\{0,10\\}\rm M_{\odot}$
for the set up shown in Figure 3. The waveform obtained for $m=10\,\rm
M_{\odot}$ from the geometrical optics limit is plotted in dashed lines.
$m=0\,\rm M_{\odot}$ corresponds to the case where there is no microlens. The
PSD of aLIGO is also plotted in the top panel in black dotted lines for
reference.

Then, we investigate how the microlens location relative to the macroimage
affects waveform suppression. Fixing $m=10\,\rm M_{\odot}$, we vary the micro
impact parameter $\eta_{m}\in\\{0.1,0.3,0.5,0.7,1.0\\}\,\theta_{m}$ and
calculate the amplification factors at different image positions. As shown in
Figure 7, the amplification factor converges to the same curve at low
frequencies, suggesting that the effects of waveform suppression are similar
regardless of the microlens location. This behavior is also observed in the
case of isolated point mass lenses (Takahashi & Nakamura, 2003; Nakamura,
1998).

We then quantify the differences introduced by microlensing to the lensed
waveforms. We first define

$\displaystyle h_{L,\mathit{full}}(f)=F(f)h(f),$ (8) $\displaystyle
h_{L,\mathit{macro}}(f)=\sqrt{\mu}\,h(f).$ (9)

Then, we define the effective micro-magnification by

$\mu_{m}=\frac{(h_{L,\mathit{full}}\,|\,h_{L,\mathit{full}})}{(h_{L,\mathit{macro}}\,|\,h_{L,\mathit{macro}})},$
(10)

where

$(a|b)=4\operatorname{Re}\int^{f_{\mathit{high}}}_{f_{\mathit{low}}}\frac{a(f)b^{*}(f)}{S_{n}(f)}$
(11)

is the inner product weighted by the power-spectral density (PSD) $S_{n}(f)$
of LIGO. By taking the noise weighted inner product, we remove the frequency
dependence, so we can now approximate the change induced by the microlens on
the amplification of the GW waveform. In Figure 8, we plot $\sqrt{\mu_{m}}$
for $m\in(1,100)\,\rm M_{\odot}$ and $\eta_{m}\in(0.1,1.0)\,\theta_{m}$. For
$h(f)$, we again use an IMRPhenomD waveform of a $30+30\,\rm M_{\odot}$ binary
black hole merger, generated with pycbc (Nitz et al., 2020). We took the
square root because the amplification to the waveform is $\sim\sqrt{\mu_{m}}$.
If the microlens were absent, $h_{L,\mathit{full}}=h_{L,\mathit{macro}}$ and
eq. 10 reduces to $\mu_{m}=1$. For lower $m$, $\mu_{m}$ does not significantly
deviate from unity. For $m\lesssim 3\,\rm M_{\odot}$, $\sqrt{\mu_{m}}$ is
larger than $1$ by $\lesssim 5\%$.

To further quantify the deviation in the waveforms, we define the waveform
overlap, which is given by

$\mathcal{O}[h_{L,\mathit{macro}},h_{L,\mathit{full}}]=\frac{(h_{L,\mathit{full}}\,|\,h_{L,\mathit{macro}})}{\sqrt{(h_{L,\mathit{full}}\,|\,h_{L,\mathit{full}})(h_{L,\mathit{macro}}\,|\,h_{L,\mathit{macro}})}},$
(12)

We then define the match
$\mathcal{M}[h_{L,\mathit{macro}},h_{L,\mathit{full}}]$ to be the maximum of
$\mathcal{O}[h_{L,\mathit{macro}},h_{L,\mathit{full}}]$ over time and phase,
and $1-\mathcal{M}[h_{L,\mathit{macro}},h_{L,\mathit{full}}]$ to be the
microlens induced mismatch. The match is also computed using the pycbc
software (Nitz et al., 2020).

Figure 9 shows a contour plot of the mismatch
$1-\mathcal{M}[h_{L,\mathit{macro}},h_{L,\mathit{full}}]$ between the
waveforms with and without microlensing as a function of $m$ and $\eta_{m}$.
For most of the parameter space concerned, the mismatch $\lesssim 1\%$ except
for $m,\eta_{m}\sim 100\,\rm M_{\odot},0.1\,\theta_{m}$. For $m\lesssim
10\,\rm M_{\odot}$, the mismatch $\lesssim 0.5\%$. We note that the waveform’s
mismatch does not encode any information about the macroimage magnification
$\mu$, which scales the signal to noise ratio (SNR) of the signal as
$\sqrt{\mu}$. However, as the mismatch is at the sub-percent level for typical
stellar masses of $m\lesssim 10\,\rm M_{\odot}$, it is expected that
microlenses would not significantly affect the measurement of other GW
parameters for typical magnifications $\mu\lesssim 10$. We note that more
complicated scenarios with a field of microlenses still exist and there are
some rare cases when we could observe extremely large macromodel
magnifications ($\mu\sim 100$). These scenarios may still induce more
substantial stellar microlensing features on the waveform (Diego et al., 2019;
Diego, 2020).

Figure 7: Same as Figure 5, but this time fixing the microlens mass at
$m=10\,\rm M_{\odot}$ while varying the micro impact parameter
$\eta_{m}\in\\{0.1,0.3,0.5,0.7,1.0\\}\,\theta_{m}$. Due to diffraction
effects, the amplification factors in the lower frequency region are
suppressed to the same limit regardless of the micro impact parameter. Figure
8: Contour plot of the effective micro-amplificaiton $\sqrt{\mu_{m}}$ for our
lens set up as a function of the microlens mass $m$ and the micro impact
parameter $\eta_{m}$. The behavior of $\sqrt{\mu_{m}}$ almost does not depend
on $\eta_{m}$ for a low $m$.  Figure 9: The mismatch between the waveforms
with and without microlensing for our lens set up as a function of the
microlens mass $m$ and the micro impact parameter $\eta_{m}$. The contour
lines of mismatch $\in\\{0.1,0.5,1.0\\}\%$ are shown. Stellar mass microlenses
of $m\lesssim 10\,\rm M_{\odot}$ have a mismatch in the subpercent level.

##  4 Discussion

###  4.1 Implications of waveform suppression

In most micro- and millilensing situations above $100\,\rm M_{\odot}$, the
geometrical optics limit describes the physical configuration well at the most
sensitive LIGO/Virgo frequencies. However, for microlenses below $100\,\rm
M_{\odot}$, wave diffraction effects become non-negligible. Wave optics
effects suppress the lensing magnification due to microlenses, especially in
common astrophysical strong lensing scenarios with microlenses $\sim 1\,\rm
M_{\odot}$ (Figure 5). We note that there are still rarer astrophysical
situations where two or more microlenses could plausibly conspire together to
mimic larger lenses (Diego et al., 2019; Diego, 2020), which is discussed in
Appendix B. However, our results are nevertheless promising, and may have
exciting implications for studies of fundamental physics and cosmology through
localization studies.

In particular, as pointed out by Oguri (2019), typical microlensing may
significantly and adversely affect the accuracy and precision of time delay
and magnification measurements of lensing transients. As the typical radii of
the objects emitting gravitational waves are small (the orbital radius of
binary compact object mergers can be as low as $R\sim 10\,\rm km$), GWs can in
principle be subject to considerable microlensing in the usual geometrical
optics limit (Oguri, 2019). However, our results confirm that the wave optics
effects suppress stellar microlensing (see also discussion on isolated lenses
in Oguri (2019)). This suggests that gravitational waves may allow for
relatively clean measurements where we need to worry less about microlensing
systematics in, e.g., localization studies (Hannuksela et al., 2020) and
possible follow-up studies thereafter (Sereno et al., 2011; Collett & Bacon,
2017; Baker & Trodden, 2017; Fan et al., 2017; Liao et al., 2017; Cao et al.,
2019; Li et al., 2019b; Hannuksela et al., 2020). Indeed, for example, the
lensed GWs may be a promising means for us to constrain cosmological
parameters, such as the Hubble constant $H_{0}$, by accurately measuring the
time delay of different strongly lensed GW images. Intriguingly, GWs are not
as vulnerable to microlensing as the EM waves from supernovae and other
explosive transients (see discussion in Oguri (2019)). Systematic errors due
to microlensing can potentially play a significant role in these measurements.
Thus, wave optics suppression can be seen as an advantage in this case.

Moreover, Hannuksela et al. (2020) demonstrated the use of quadruply lensed
GWs to localize the host galaxy by matching the GW image properties with the
lensed galaxy systems found in the direction of the GW in an electromagnetic
follow-up. Microlensing could induce a significant error on the lensing time-
delay and magnification, hampering the match between the image properties
(magnifications, time delays) found in the EM and the GW channels. Our work
suggests that these localization studies may suffer from stellar microlensing
in the common astrophysical scenarios to a much lesser degree, which may spell
improvements for these localization studies. Subsequently, microlensing may
pose less of a problem in measuring $H_{0}$ through methods based on the
distance-ladder or the time-delay distance. However, there is still more work
in investigating more complex microlensing scenarios, such as fields of
microlenses near the saddle point macroimage (type-II image). Nevertheless,
our results are very promising.

Finally, our results illustrate the mass ranges for which wave optics effect
becomes significant. As shown in Figure 5, both the magnitude and argument of
the amplification factor deviate significantly from the geometric optics model
even at a microlens mass of $30\,\rm M_{\odot}$. In Appendix A, we show that
the validity of the geometrical optics approximation depends at least on the
mass of the lens and the frequency range in question, and the approximation
might be good enough for LIGO detections for microlenses of mass $m\gtrsim
300\,\rm M_{\odot}$.

###  4.2 Non-minimum macroimages

In this paper, we focused on the amplification factors of microlenses near the
macromodel minimum point (a type-I image) of the time delay function. This is
for computational simplicity, as only the area immediately surrounding the
minimum contributes to the diffraction integral around the macroimage time
delay.

However, for a type-II images (saddle points), the regions of the lens plane
with time delay
$t_{d}(\vec{\xi},\vec{\eta})\in[t_{\mathit{im}}-\tau,t_{\mathit{im}}+\tau]$
correspond to thin strips centering around the contour lines
$t_{d}(\vec{\xi},\vec{\eta})=t_{\mathit{im}}$. As the macroimage is a saddle
point, such a contour line can have an arbitrary length that depends on the
lens model. In SIS lensing with typical galactic masses, the relevant contours
at the saddle point will be vast, making the calculation of eq. 6 very
computationally expensive. It has been shown in Diego et al. (2019) that
microlensing around type-II images have qualitative differences from that
around type-I ones. We will develop methods to tackle type-II images of
general lens configurations in future work.

##  5 Conclusions

We have shown that a wave optics treatment of the full lens model is required
for stellar-mass microlensing on a galactic macromodel. Due to diffraction
effects, microlensing is suppressed, making the contributions by stellar
lenses of mass $\sim 1\,\rm M_{\odot}$ small regardless of their proximity to
the macrolensed GW beams. This allows for more accurate measurements of the
magnification and time delay of lensed GW images and better confidence in
localization studies. On the flip side, we might require third-generation
detectors to use GW lensing to probe objects of masses $\mathcal{O}(1\,\rm
M_{\odot})$. Future work will include the cases of microlensing near type-II
macroimages.

## Acknowledgements

We would like to thank Apratim Ganguly, Ajit Mehta, Anupreeta More, Eungwang
Seo, and K. Haris for their insightful comments. Furthermore, we would like to
thank Miguel Zumalacarregui for helping with the calculation of the
diffraction integral, and also Thomas Collett for related discussions. OAH is
supported by the research program of the Netherlands Organization for
Scientific Research (NWO). TGFL is partially supported by grants from the
Research Grants Council of Hong Kong (Project No. 14306218), Research
Committee of the Chinese University of Hong Kong, and the Croucher Foundation
of Hong Kong. We are grateful for computational resources provided by the LIGO
Laboratory and supported by the National Science Foundation Grants PHY-0757058
and PHY-0823459.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding authors.

## References

  * Baker & Trodden (2017) Baker T., Trodden M., 2017, Phys. Rev. D, 95, 063512
  * Cao et al. (2019) Cao S., Qi J., Cao Z., Biesiada M., Li J., Pan Y., Zhu Z.-H., 2019, Sci. Rep., 9, 11608
  * Chen & Holz (2016) Chen H.-Y., Holz D. E., 2016, arXiv preprint arXiv:1612.01471
  * Christian et al. (2018) Christian P., Vitale S., Loeb A., 2018, Phys. Rev. D, 98, 103022
  * Collett & Bacon (2017) Collett T. E., Bacon D., 2017, Phys. Rev. Lett., 118, 091101
  * Dai & Venumadhav (2017) Dai L., Venumadhav T., 2017, arXiv preprint arXiv:1702.04724
  * Dai et al. (2020) Dai L., Zackay B., Venumadhav T., Roulet J., Zaldarriaga M., 2020, arXiv preprint arXiv:2007.12709
  * Deguchi & Watson (1986) Deguchi S., Watson W. D., 1986, Phys. Rev. D, 34, 1708
  * Diego (2020) Diego J. M., 2020, Phys. Rev. D, 101, 123512
  * Diego et al. (2019) Diego J., Hannuksela O., Kelly P., Broadhurst T., Kim K., Li T., Smoot G., Pagano G., 2019, Astron. Astrophys., 627, A130
  * Ezquiaga et al. (2020) Ezquiaga J. M., Holz D. E., Hu W., Lagos M., Wald R. M., 2020, arXiv preprint arXiv:2008.12814
  * Fan et al. (2017) Fan X.-L., Liao K., Biesiada M., Piorkowska-Kurpas A., Zhu Z.-H., 2017, Phys. Rev. Lett., 118, 091102
  * Goyal et al. (2020) Goyal S., Haris K., Mehta A. K., Ajith P., 2020, arXiv preprint arXiv:2008.07060
  * Hannuksela et al. (2019) Hannuksela O., Haris K., Ng K., Kumar S., Mehta A., Keitel D., Li T., Ajith P., 2019, Astrophys. J. Lett., 874, L2
  * Hannuksela et al. (2020) Hannuksela O. A., Collett T. E., Çalışkan M., Li T. G. F., 2020, Monthly Notices of the Royal Astronomical Society, 498, 3395
  * Haris et al. (2018) Haris K., Mehta A. K., Kumar S., Venumadhav T., Ajith P., 2018, arXiv preprint arXiv:1807.07062
  * Jung & Shin (2019) Jung S., Shin C. S., 2019, Phys. Rev. Lett., 122, 041103
  * Lai et al. (2018) Lai K.-H., Hannuksela O. A., Herrera-Martín A., Diego J. M., Broadhurst T., Li T. G., 2018, Phys. Rev. D, 98, 083005
  * Li et al. (2018) Li S.-S., Mao S., Zhao Y., Lu Y., 2018, Mon. Not. Roy. Astron. Soc., 476, 2220
  * Li et al. (2019a) Li A. K., Lo R. K., Sachdev S., Chan C., Lin E., Li T. G., Weinstein A. J., 2019a, arXiv preprint arXiv:1904.06020
  * Li et al. (2019b) Li Y., Fan X., Gou L., 2019b, Astrophys. J., 873, 37
  * Liao et al. (2017) Liao K., Fan X.-L., Ding X.-H., Biesiada M., Zhu Z.-H., 2017, Nature Commun., 8, 1148
  * Liu et al. (2020) Liu X., Hernandez I. M., Creighton J., 2020, arXiv preprint arXiv:2009.06539
  * McIsaac et al. (2020) McIsaac C., Keitel D., Collett T., Harry I., Mozzon S., Edy O., Bacon D., 2020, Phys. Rev. D, 102, 084031
  * Nakamura (1998) Nakamura T. T., 1998, Phys. Rev. Lett., 80, 1138
  * Ng et al. (2018) Ng K. K., Wong K. W., Broadhurst T., Li T. G., 2018, Phys. Rev. D, 97, 023012
  * Nitz et al. (2020) Nitz A., et al., 2020, gwastro/pycbc: PyCBC release v1.16.11, doi:10.5281/zenodo.4134752, https://doi.org/10.5281/zenodo.4134752
  * Oguri (2018) Oguri M., 2018, Mon. Not. Roy. Astron. Soc., 480, 3842
  * Oguri (2019) Oguri M., 2019, Reports on Progress in Physics, 82, 126901
  * Ohanian (1974) Ohanian H., 1974, Int. J. Theor. Phys., 9, 425
  * Pagano et al. (2020) Pagano G., Hannuksela O. A., Li T. G., 2020, Astron. Astrophys., 643, A167
  * Pang et al. (2020) Pang P. T., Hannuksela O. A., Dietrich T., Pagano G., Harry I. W., 2020, Mon. Not. Roy. Astron. Soc.
  * Peters (1974) Peters P. C., 1974, Phys. Rev. D, 9, 2207
  * Robertson et al. (2020) Robertson A., Smith G. P., Massey R., Eke V., Jauzac M., Bianconi M., Ryczanowski D., 2020, Mon. Not. Roy. Astron. Soc.
  * Ryczanowski et al. (2020) Ryczanowski D., Smith G. P., Bianconi M., Massey R., Robertson A., Jauzac M., 2020, Mon. Not. Roy. Astron. Soc., 495, 1666
  * Schneider et al. (1992) Schneider P., Ehlers J., Falco E. E., 1992, Gravitational Lenses. Springer
  * Sereno et al. (2011) Sereno M., Jetzer P., Sesana A., Volonteri M., 2011, Mon. Not. Roy. Astron. Soc., 415, 2773
  * Smith et al. (2017) Smith G., et al., 2017, IAU Symp., 338, 98
  * Smith et al. (2018) Smith G. P., Jauzac M., Veitch J., Farr W. M., Massey R., Richard J., 2018, Mon. Not. Roy. Astron. Soc., 475, 3823
  * Smith et al. (2019) Smith G. P., Robertson A., Bianconi M., Jauzac M., 2019, arXiv preprint arXiv:1902.05140
  * Takahashi & Nakamura (2003) Takahashi R., Nakamura T., 2003, Astrophys. J., 595, 1039
  * Thorne (1982) Thorne K., 1982, in Les Houches Summer School on Gravitational Radiation. pp 1–57
  * Ulmer & Goodman (1994) Ulmer A., Goodman J., 1994, arXiv preprint astro-ph/9406042
  * Wang et al. (1996) Wang Y., Stebbins A., Turner E. L., 1996, Phys. Rev. Lett., 77, 2875
  * Yu et al. (2020) Yu H., Zhang P., Wang F.-Y., 2020, Mon. Not. Roy. Astron. Soc., 497, 204

##  Appendix A Accuracy of the geometrical optics approximation

It has generally been accepted that the geometrical optics limit approximates
the exact results better for higher frequencies $f$ and higher (micro)lens
masses $m$. In Figure 10, we compare the full wave optics calculation (solid
lines) and geometrical optics approximate (dashed lines) of the amplification
factor $F$ for the set up in Figure 3 but with $m=100\,\rm M_{\odot}$.
Comparing Figures 4 and 10, we see that the geometrical optics approximation
approaches the wave optics calculations at a lower frequency for the more
massive case $m=100\,\rm M_{\odot}$. More explicitly, we can quantify the
error of the geometrical optics limit by taking the norm of the difference
between the wave optics amplification factor $F_{\mathit{wave}}$ and the
geometrical optics amplification factor $F_{\mathit{geom}}$. In Figure 11 we
plot the results for $m\in\\{1,30,100,300,500\\}\,\rm M_{\odot}$ (normalized
by $\sqrt{\mu}$). For $m\in\\{300,500\\}\,\rm M_{\odot}$, the normalized error
in the geometrical optics approximation is $<0.1$ for most of the frequency
range sensitive for LIGO, meaning that the approximation might be good enough
for $m\gtrsim 300\,\rm M_{\odot}$ depending on the accuracy required.

Figure 10: Same as Figure 4, but with a microlens mass of $m=100\,\rm
M_{\odot}$ instead. There is still a clear mismatch between the geometrical
and wave optics limits for $f\lesssim 10^{2}\,\rm Hz$. Figure 11: Deviation
between the wave optics and geometrical optics results of the amplification
factor for $m\in\\{1,30,100,300,500\\}\,\rm M_{\odot}$. The shaded region
corresponds to $|F_{\mathit{wave}}-F_{\mathit{geom}}|/\sqrt{\mu}<0.1$. The
curves for $m\in\\{300,500\\}\,\rm M_{\odot}$ lie wtihin this region for the
most sensitive frequency range $f\sim\mathcal{O}(10^{2}\,\rm Hz)$ of LIGO.

##  Appendix B Two microlenses near the macroimage

Although microlensing due to a single star might not affect GWs significantly,
we are interested in the case of GW microlensing by multiple stars as well.
Here we study the case of a GW ray passing through the Einstein radii of two
stars. We use the same $10^{10}\,\rm M_{\odot}$ SIS macromodel as that in the
main text, and restrict ourselves to the minimum macroimage. As shown in
Figure 12, we place two microlenses, both of mass $1\,\rm M_{\odot}$, at
random positions $\vec{\eta}_{m1}$ and $\vec{\eta}_{m2}$ within an Einstein
radius $\theta_{m}$ of the macroimage ($\theta_{m}$ corresponds to the
Einstein radius of a lens with mass $m=1\,\rm M_{\odot}$). We simulate $1000$
of these microlens pairs, calculate the amplification factor for each
configuration, and compute the micro magnification factor as in eq. 10 and the
microlens induced mismatch as in Figure 9. The results are shown in Figure 13.
The effective micro amplification $\sqrt{\mu_{m}}$ increases when compared to
the single microlens case (see Figure 8), but it still deviates little from
unity, with $\sqrt{\mu_{m}}<1.04$. Likewise, the waveform mismatch increases
only slightly due to the additional microlens, remaining in the $<0.1\%$
range. This shows that even for microlensing by two stars, the effects are
still relatively small. In general, however, we find that both the
microlensing magnification and the induced mismatch increase with additional
microlenses.

Figure 12: Microlensing set up for the two microlenses case. We place two
microlenses, both of mass $m=1\,\rm M_{\odot}$, at random positions
$\vec{\eta}_{m1}$ and $\vec{\eta}_{m2}$ within an Einstein radius $\theta_{m}$
from the macroimage.

Figure 13: The effective micro amplification $\sqrt{\mu_{m}}$ and microlens
induced mismatch for the $1000$ simulated double microlens configurations.
$\sqrt{\mu_{m}}<1.04$ and mismatch $<0.1\%$ for all of these configurations.

##  Appendix C Microlenses placed at other angles from the macroimage

In the main text, we only considered the case where the microlenses are placed
on the straight line between the SIS lens center and the minimum macroimage.
Here, with the same macromodel set up, we also vary the direction we place the
microlenses from the macroimage, with
$\phi_{m}\in\\{45\degree,90\degree,135\degree,180\degree\\}$, where $\phi_{m}$
is the relative angle from the SIS center to the microlens measured at the
macroimage as shown in Figure 14. The case examined in the main text
corresponds to $\phi_{m}=0\degree$, while the results for the range
$180\degree<\phi_{m}<360\degree$ are similar to those examined here by
symmetry of the model. The deviation in micro magnification and the mismatch
for the above cases are plotted in Figure 15, showing that the microlens
displacement relative to the macroimage does not affect the results
significantly for $m\lesssim 10\,\rm M_{\odot}$.

Figure 14: Varying the angular direction with respect to the macroimage for
which we place the microlens. The SIS center is to the left of the figure and
$\phi_{m}$ is the angle between the microlens and the SIS center measured at
the macroimage. The dashed and solid circles correspond to
$\phi_{m}\in\\{0\degree,45\degree,90\degree,135\degree,180\degree\\}$. Figure
15: Micro amplification factor $\mu_{m}$ (left column) and mismatch (right
column) for $\phi_{m}\in\\{45\degree,90\degree,135\degree,180\degree\\}$ (from
top to bottom). For $m\lesssim 10\,\rm M_{\odot}$, the features are
qualitatively similar to those in the case discussed in the main text (see
Figures 8 and 9). For $m\gtrsim 10$, varying $\phi_{m}$ does introduce some
differences, meaning that the ability to detect microlenses in this higher
mass range would at least depend on the three parameters $m$, $\eta_{m}$, and
$\phi_{m}$


and in the single-component limit ($N_{\text{spec}}\rightarrow 1$), where
$P_{i}\equiv P_{0}=n_{0}/\beta$, we reproduce the result from Ref. [21]:

$\displaystyle\mathcal{A}^{(2)}_{11,00}\rightarrow\mathcal{A}^{(2)}_{00}=\frac{3}{5}\sigma_{\text{tot}}n_{0}\;.$
(198)

Here, Eq. (197) defines the entries of an $N_{\text{spec}}$-dimensional
rectangular, regular matrix. The elements of its inverse are:

$\displaystyle\tau^{(2)}_{si,00}=\frac{5}{\sigma_{\text{tot}}\beta}\frac{1}{12P^{2}}\left(3\delta_{si}P+P_{s}\right)\;.$
(199)

#### Vector moments ($\ell=1$):

The relevant contractions of the $\mathcal{X}$-tensor for the vector moments
are

$\displaystyle\sum_{j,a,b\,=\,1}^{N_{\text{spec}}}\Delta_{\mu_{1}\nu_{1}}\mathcal{X}_{sijab,r}^{\mu_{1}\nu_{1}}=$
$\displaystyle\;\sigma_{\text{tot}}\sum_{j\,=\,1}^{N_{\text{spec}}}\left(\delta_{js}-\delta_{is}\right)\Bigg{(}\frac{1}{3}I_{i,r+1,0}I_{j,20}+I_{i,r+2,0}I_{j,10}\Bigg{)}\;,$
(200)
$\displaystyle\sum_{j,a,b\,=\,1}^{N_{\text{spec}}}\Delta_{\mu_{1}\nu_{1}}u_{\nu_{2}}\mathcal{X}_{sijab,r}^{\mu_{1}\nu_{1}\nu_{2}}=$
$\displaystyle\;\sigma_{\text{tot}}\sum_{j\,=\,1}^{N_{\text{spec}}}\left[\frac{2}{9}\left(2\delta_{js}-\delta_{is}\right)I_{i,r+1,0}I_{j,30}\right.$
$\displaystyle\left.+\frac{4}{9}\left(\delta_{js}+\delta_{is}\right)I_{i,r+2,0}I_{j,20}+\frac{6}{9}\left(\delta_{js}-2\delta_{is}\right)I_{i,r+3,0}I_{j,10}\right]\;,$
(201)

thus the elements of the collision matrix (176) read in the
$(10+4N_{q})$-moment approximation (where $N_{1}=1$):

$\displaystyle\mathcal{A}^{(1)}_{is,r0}$
$\displaystyle=-\frac{1}{2}\sum_{n^{\prime}\,=\,0}^{1}\sum_{m\,=\,0}^{n^{\prime}}\frac{a^{(1)}_{s,n^{\prime}0}a^{(1)}_{s,n^{\prime}m}}{3\,I_{s,21}}\sum_{j,a,b\,=\,1}^{N_{\text{spec}}}u_{\nu_{2}}\cdots
u_{\nu_{1+m}}\Delta^{\mu_{1}\nu_{1}}\mathcal{X}_{sijab,r}^{\mu_{1}\nu_{1}\ldots\nu_{1+m}}$
$\displaystyle=-\frac{\sigma_{\text{tot}}}{6\,I_{s,21}}\sum_{j\,=\,1}^{N_{\text{spec}}}\Bigg{\\{}\left(a^{(1)}_{s,00}a^{(1)}_{s,00}+a^{(1)}_{s,10}a^{(1)}_{s,10}\right)\left(\delta_{js}-\delta_{is}\right)\Bigg{(}\frac{1}{3}I_{i,r+1,0}I_{j,20}+I_{i,r+2,0}I_{j,10}\Bigg{)}$
$\displaystyle+a^{(1)}_{s,10}a^{(1)}_{s,11}\left[\frac{2}{9}\left(2\delta_{js}-\delta_{is}\right)I_{i,r+1,0}I_{j,30}+\frac{4}{9}\left(\delta_{js}+\delta_{is}\right)I_{i,r+2,0}I_{j,20}+\frac{6}{9}\left(\delta_{js}-2\delta_{is}\right)I_{i,r+3,0}I_{j,10}\right]\Bigg{\\}}\;,$
(202)

and

$\displaystyle\mathcal{A}^{(1)}_{is,r1}=$
$\displaystyle\;-\frac{1}{2}\sum_{m\,=\,0}^{1}\frac{a^{(1)}_{s,11}a^{(1)}_{s,1m}}{3\,I_{s,21}}\sum_{j,a,b\,=\,1}^{N_{\text{spec}}}u_{\nu_{2}}\cdots
u_{\nu_{1+m}}\Delta^{\mu_{1}\nu_{1}}\mathcal{X}_{sijab,m}^{\mu_{1}\nu_{1}\ldots\nu_{1+m}}$
$\displaystyle=$
$\displaystyle\;-\frac{\sigma_{\text{tot}}}{6\,I_{s,21}}\sum_{j\,=\,1}^{N_{\text{spec}}}\Bigg{\\{}a^{(1)}_{s,11}a^{(1)}_{s,10}\left(\delta_{js}-\delta_{is}\right)\Bigg{(}\frac{1}{3}I_{i,r+1,0}I_{j,20}+I_{i,r+2,0}I_{j,10}\Bigg{)}$
$\displaystyle+a^{(1)}_{s,11}a^{(1)}_{s,11}\left[\frac{2}{9}\left(2\delta_{js}-\delta_{is}\right)I_{i,r+1,0}I_{j,30}+\frac{4}{9}\left(\delta_{js}+\delta_{is}\right)I_{i,r+2,0}I_{j,20}+\frac{6}{9}\left(\delta_{js}-2\delta_{is}\right)I_{i,r+3,0}I_{j,10}\right]\Bigg{\\}}\;.$
(203)

After some calculations, the relevant matrix elements read in terms of
pressure:

$\displaystyle\mathcal{A}^{(1)}_{is,00}$
$\displaystyle=-\frac{4}{9}\sigma_{\text{tot}}\beta\sum_{j\,=\,1}^{N_{\text{spec}}}\frac{P_{i}P_{j}}{P_{s}}\left(\delta_{js}-2\delta_{is}\right)=\frac{4}{9}\sigma_{\text{tot}}\beta\left(2\delta_{is}P-P_{i}\right)\;,$
(204) $\displaystyle\mathcal{A}^{(1)}_{is,10}$
$\displaystyle=-\frac{1}{2}\sigma_{\text{tot}}\sum_{j\,=\,1}^{N_{\text{spec}}}\frac{P_{i}P_{j}}{P_{s}}\left(\delta_{js}-\delta_{is}\right)=\frac{1}{2}\sigma_{\text{tot}}\left(\delta_{is}P-P_{i}\right)\;,$
(205) $\displaystyle\mathcal{A}^{(1)}_{is,01}$
$\displaystyle=-\frac{1}{18}\sigma_{\text{tot}}\beta^{2}\sum_{j\,=\,1}^{N_{\text{spec}}}\frac{P_{i}P_{j}}{P_{s}}\left(\delta_{js}+\delta_{is}\right)=-\frac{1}{18}\sigma_{\text{tot}}\beta^{2}\left(\delta_{is}P+P_{i}\right)\;,$
(206) $\displaystyle\mathcal{A}^{(1)}_{is,11}$
$\displaystyle=-\frac{1}{2}\sigma_{\text{tot}}\beta\sum_{j\,=\,1}^{N_{\text{spec}}}\frac{P_{i}P_{j}}{P_{s}}\left(\delta_{js}-\delta_{is}\right)=\frac{1}{2}\sigma_{\text{tot}}\beta\left(\delta_{is}P-P_{i}\right)\;.$
(207)

Equations (204) – (207) are the elements of the four
$N_{\text{spec}}$-dimensional rectangular block matrices of the
$(2N_{\text{spec}}\times 2N_{\text{spec}})$-matrix $\mathcal{A}^{(1)}$. Its
single-component limit ($N_{\text{spec}}=1$) is consistent with Ref. [21]:

$\displaystyle\mathcal{A}^{(1)}\equiv\begin{pmatrix}\mathcal{A}^{(1)}_{11,00}&\mathcal{A}^{(1)}_{11,01}\\\
\mathcal{A}^{(1)}_{11,10}&\mathcal{A}^{(1)}_{11,11}\end{pmatrix}=\frac{1}{9}\sigma_{\text{tot}}\beta
P_{0}\begin{pmatrix}4&-\beta\\\ 0&0\end{pmatrix}\;.$ (208)

We observe that the $\mathcal{A}^{(1)}$ matrix is singular even in the single-
component limit. This is due to the momentum-conservation equation. In order
to construct the inverse matrix $\tau^{(1)}$, we follow the steps presented in
Appendix G by introducing the reduced matrix $\tilde{\mathcal{A}}^{(1)}$
(therefore effectively removing the irreducible moment $\rho^{\mu}_{1,1}$),
inverting it, and adding zero elements corresponding to the originally removed
rows and columns to that inverse, yielding $\tau^{(1)}$. We illustrate this
procedure in the single-component limit. The reduced matrix and its inverse
then just consist of one entry,

$\displaystyle\tilde{\mathcal{A}}^{(1)}=\frac{4}{9}\sigma_{\text{tot}}\beta
P_{0}\;,\quad\text{and}\quad\tilde{\mathcal{\tau}}^{(1)}=\frac{9}{4}\frac{1}{\sigma_{\text{tot}}\beta
P_{0}}\;.$ (209)

Adding zero elements yields the final inverse

$\displaystyle\mathcal{\tau}^{(1)}=\frac{9}{4}\frac{1}{\sigma_{\text{tot}}\beta
P_{0}}\begin{pmatrix}1&0\\\ 0&0\end{pmatrix}\;.$ (210)

In the multicomponent case, we find the following entries of the
$2N_{\text{spec}}$-dimensional rectangular inverse matrix:

$\displaystyle\tau^{(1)}_{si,00}$
$\displaystyle=\frac{9}{68\beta\sigma_{\text{tot}}P^{2}}(8\delta_{si}P+9P_{s})\;,$
(211) $\displaystyle\tau^{(1)}_{si,01}$
$\displaystyle=\frac{2}{17\sigma_{\text{tot}}P}(\delta_{si}-\delta_{sN_{\text{spec}}})\;,$
(212) $\displaystyle\tau^{(1)}_{si,10}$
$\displaystyle=\frac{18}{17\sigma_{\text{tot}}\beta^{2}P^{2}}(1-\delta_{sN_{\text{spec}}})\left(P_{s}-\delta_{si}P\right)\;,$
(213) $\displaystyle\tau^{(1)}_{si,11}$
$\displaystyle=\frac{32}{17\sigma_{\text{tot}}\beta
P}\delta_{si}(1-\delta_{sN_{\text{spec}}})\;.$ (214)

We note that the elements $\tau^{(1)}_{si,nr}$ are indeed constructed in a way
that they vanish in the cases $n=1$ and $s=N_{\text{spec}}$, or $r=1$ and
$i=N_{\text{spec}}$ (i.e. the row and column which was originally removed from
$\mathcal{A}^{(1)}$). We remind the reader that adding these zeros simplifies
our notation in this work.

###  F.2 Transport coefficients

Now that the collision matrix has been determined, we can proceed to calculate
the transport coefficients of the theory. We remark that the scalar moments
have not been discussed since the bulk viscous pressure vanishes identically
in the ultrarelativistic (massless) case,
$\Pi\equiv\sum\limits_{s\,=\,1}^{N_{\text{spec}}}\frac{m_{s}^{2}}{3}\rho_{s,0}=0$.
The coefficients in the $(10+4N_{q})$-moment approximation ($N_{0}=2$,
$N_{1}=1$, $N_{2}=0$) have been defined in Sec. III.4 and Appendix C. We
evaluate them in the ultrarelativistic scenario, where all mass terms vanish.
Since the transport coefficients are defined via the coefficients listed in
Eqs. (98) – (100), these have to be evaluated first. For this, we introduce
short-hand notations for the charge concentration of type $q$, and the
concentration of the charge combination $qq^{\prime}$ in the system,
respectively:

$\displaystyle\mathbf{c}_{q}\equiv\sum_{j\,=\,1}^{N_{\text{spec}}}q_{j}\frac{P_{j}}{P}\;,\quad\text{and}\quad\mathbf{c}_{qq^{\prime}}\equiv\sum_{j\,=\,1}^{N_{\text{spec}}}q_{j}q^{\prime}_{j}\frac{P_{j}}{P}\;.$
(215)

Further, the derivatives in temperature and chemical potential of the weighted
partial pressures of a classical gas read

$\displaystyle\frac{\partial}{\partial\beta}\left(\frac{P_{i}}{P}\right)=0\quad\text{and}\quad\frac{\partial}{\partial\alpha_{q}}\left(\frac{P_{i}}{P}\right)=\frac{P_{i}}{P}\left(q_{i}-\mathbf{c}_{q}\right),$
(216)

and from this the derivatives in the charge concentration follow

$\displaystyle\frac{\partial\mathbf{c}_{q}}{\partial\beta}=0\quad\text{and}\quad\frac{\partial\mathbf{c}_{q}}{\partial\alpha_{q^{\prime}}}=\mathbf{c}_{qq^{\prime}}-\mathbf{c}_{q}\mathbf{c}_{q^{\prime}}\;.$
(217)

respectively. The relevant expressions for the vector and tensor moments are
then obtained from Eqs. (98) – (100) as:

$\displaystyle\eta_{s,0}$
$\displaystyle\equiv\sum_{i\,=\,1}^{N_{\text{spec}}}\sum_{r\,=\,0}^{N_{2}}\tau_{si,0r}^{\left(2\right)}\alpha_{s,r}^{\left(2\right)}=\frac{4}{3}\frac{P_{s}}{P}\frac{1}{\sigma_{\text{tot}}\beta}\;,$
(218) $\displaystyle\kappa_{s,0,q}$
$\displaystyle\equiv\sum_{i\,=\,1}^{N_{\text{spec}}}\sum_{r\,=\,0}^{N_{1}}\tau_{si,0r}^{\left(1\right)}\alpha_{i,r,q}^{\left(1\right)}=\frac{8}{17\,\sigma_{\text{tot}}}\frac{P_{s}}{P}\left(q_{s}-\frac{77}{128}\mathbf{c}_{q}\right)\;,$
(219) $\displaystyle\kappa_{s,1,q}$
$\displaystyle\equiv\sum_{i\,=\,1}^{N_{\text{spec}}}\sum_{r\,=\,0}^{N_{1}}\tau_{si,1r}^{\left(1\right)}\alpha_{i,r,q}^{\left(1\right)}=\frac{26}{17\sigma_{\text{tot}}\beta}(1-\delta_{sN_{\text{spec}}})\frac{P_{s}}{P}\left(q_{s}-\mathbf{c}_{q}\right)\;.$
(220)

From this, the shear viscosity and the diffusion-coefficient matrix
immediately follow:

$\displaystyle\eta$
$\displaystyle\equiv\sum_{s\,=\,1}^{N_{\text{spec}}}\eta_{s,0}=\frac{4}{3}\frac{1}{\sigma_{\text{tot}}\beta}\;,$
(221) $\displaystyle\kappa_{qq^{\prime}}$
$\displaystyle\equiv\sum_{s\,=\,1}^{N_{\text{spec}}}q_{s}\kappa_{s,0,q^{\prime}}=\frac{8}{17\,\sigma_{\text{tot}}}\left(\mathbf{c}_{qq^{\prime}}-\frac{77}{128}\mathbf{c}_{q}\mathbf{c}_{q^{\prime}}\right)\;.$
(222)

It should be noted that the derivative of the above diffusion coefficients in
inverse temperature, and therefore also from its inverse, vanishes,

$\displaystyle\frac{\partial}{\partial\beta}\kappa_{qq^{\prime}}=0\quad\text{and}\quad\frac{\partial}{\partial\beta}\left(\kappa^{-1}\right)_{qq^{\prime}}=0\;.$
(223)

We note that the result for the shear viscosity is the same as in Ref. [21],
while the diffusion coefficients are different since they depend on the
various charges in the system, which are not taken into account in a single-
component approach. However, in the limit where there is only one conserved
particle species in the system (i.e. $N_{q}=1$ and $q_{s}=1$ and therefore
$\mathbf{c}_{q}=\mathbf{c}_{qq^{\prime}}=1$), the obtained expression is
equivalent to the value derived in Ref. [21],
$\kappa=3/(16\sigma_{\text{tot}})$. We also remark that from the above
equation we read off that the diffusion-coefficient matrix is symmetric,
$\kappa_{qq^{\prime}}=\kappa_{q^{\prime}q}$ as shown by the Onsager reciprocal
relations[12, 13]. The relevant weighted base coefficients then read

$\displaystyle\bar{\eta}_{i,0}$
$\displaystyle\equiv\frac{\eta_{i,0}}{\eta}=\frac{P_{i}}{P}\;,$ (224)
$\displaystyle\bar{\eta}_{i,-1}$
$\displaystyle\equiv\mathcal{F}^{(2)}_{i,1,0}\bar{\eta}_{i,0}=\frac{\beta}{5}\frac{P_{i}}{P}\;,$
(225) $\displaystyle\bar{\kappa}^{(q)}_{i,0}$
$\displaystyle\equiv\sum_{q^{\prime}}^{\left\\{B,Q,S\right\\}}\kappa_{i,0,q^{\prime}}\left(\kappa^{-1}\right)_{q^{\prime}q}=\frac{8}{17\sigma_{\text{tot}}}\sum_{q^{\prime}}^{\left\\{B,Q,S\right\\}}\left(\kappa^{-1}\right)_{q^{\prime}q}\frac{P_{i}}{P}\left(q^{\prime}_{i}-\frac{77}{128}\mathbf{c}_{q^{\prime}}\right)\;,$
(226) $\displaystyle\bar{\kappa}^{(q)}_{i,1}$
$\displaystyle=\frac{26}{17\beta\sigma_{\text{tot}}}\sum_{q^{\prime}}^{\left\\{B,Q,S\right\\}}\left(\kappa^{-1}\right)_{q^{\prime}q}\frac{P_{i}}{P}\left(q^{\prime}_{i}-\mathbf{c}_{q^{\prime}}\right)\left(1-\delta_{iN_{\text{spec}}}\right)\;.$
(227)

For the relaxation times defined in Eqs. (116) and (119) we get:

$\displaystyle\tau_{\pi}$
$\displaystyle=\frac{5}{3}\frac{1}{\sigma_{\text{tot}}\beta P}\;,$ (228)
$\displaystyle\tau_{q^{\prime}q}$
$\displaystyle=\frac{9}{68\beta\sigma^{2}_{\text{tot}}P}\sum_{q^{\prime\prime}}^{\left\\{B,Q,S\right\\}}\left(\kappa^{-1}\right)_{q^{\prime\prime}q}\left(\frac{784}{153}\mathbf{c}_{q^{\prime\prime}q^{\prime}}-\frac{4741}{2448}\mathbf{c}_{q^{\prime\prime}}\mathbf{c}_{q^{\prime}}\right)\overset{\text{single}}{\rightarrow}\frac{9}{4}\frac{1}{\beta\sigma_{\text{tot}}P}\;.$
(229)

We again note that also the relaxation-time matrix is symmetric,
$\tau_{qq^{\prime}}=\tau_{q^{\prime}q}$. It is further apparent that the shear
relaxation time is equivalent to the value derived in Ref. [21]. In the case
where only one conserved particle species is present (see above), such an
equivalence is also recovered for the diffusive relaxation time,
$\tau_{q}=27/(64\beta\sigma_{\text{tot}}^{2}P\kappa)=9/(4\beta\sigma_{\text{tot}}P)$,
as indicated with the notation “$\overset{\text{single}}{\rightarrow}$”. In
the following we will keep this notation, and continue with the second-order
coefficients. For the coefficients in the shear-stress tensor equation,
defined in Appendix C.3, we get:

$\displaystyle\delta_{\pi\pi}$ $\displaystyle=\frac{4}{3}\tau_{\pi}\;,$ (230)
$\displaystyle\tau_{\pi\pi}$ $\displaystyle=\frac{10}{7}\tau_{\pi}\;,$ (231)
$\displaystyle\ell_{\pi V}^{(q)}$
$\displaystyle=\frac{52}{51}\frac{1}{\beta^{2}\sigma_{\text{tot}}^{2}P}\sum_{q^{\prime}}^{\\{B,Q,S\\}}(\kappa^{-1})_{q^{\prime}q}\left(\mathbf{c}_{q^{\prime}}-q^{\prime}_{N_{\text{spec}}}\right)\frac{P_{N_{\text{spec}}}}{P}\overset{\text{single}}{\rightarrow}0\;,$
(232) $\displaystyle\tau_{\pi V}^{(q)}$
$\displaystyle=4\beta\frac{\partial}{\partial\beta}\ell^{(q)}_{\pi
V}-7\ell^{(q)}_{\pi V}\overset{\text{single}}{\rightarrow}0\;,$ (233)
$\displaystyle\lambda_{\pi V}^{(q,q^{\prime})}$
$\displaystyle=\frac{\partial}{\partial\alpha_{q^{\prime}}}\ell^{(q)}_{\pi
V}+\frac{\beta}{4}\mathbf{c}_{q^{\prime}}\frac{\partial}{\partial\beta}\ell^{(q)}_{\pi
V}-2\mathbf{c}_{q^{\prime}}\ell^{(q)}_{\pi
V}\overset{\text{single}}{\rightarrow}0\;.$ (234)

For the coefficients in the vector equations of motion (see Appendix C.2) we
derive:

$\displaystyle\delta_{VV}^{(q^{\prime},q)}=$
$\displaystyle\;\frac{1}{\beta\sigma^{2}_{\text{tot}}P}\sum_{q^{\prime\prime}}^{\\{B,Q,S\\}}\left(\kappa^{-1}\right)_{q^{\prime\prime}q}\left[\frac{640}{867}\mathbf{c}_{q^{\prime}q^{\prime\prime}}-\frac{17551}{55488}\mathbf{c}_{q^{\prime}}\mathbf{c}_{q^{\prime\prime}}+\frac{52}{289}\left(\mathbf{c}_{q^{\prime}q^{\prime\prime}}-\mathbf{c}_{q^{\prime}}\mathbf{c}_{q^{\prime\prime}}\right)\frac{4P}{\beta}\left(\mathcal{T}_{00}+\frac{\beta}{4}\sum_{q^{\prime\prime\prime}}^{\\{B,Q,S\\}}\mathcal{T}_{0q^{\prime\prime\prime}}\mathbf{c}_{q^{\prime\prime\prime}}\right)\right]$
$\displaystyle\quad-\frac{4}{\beta\sigma^{2}_{\text{tot}}}\sum_{q^{\prime\prime}}^{\\{B,Q,S\\}}\left\\{\frac{\partial}{\partial\alpha_{q^{\prime\prime}}}\left[\sum_{q^{\prime\prime\prime}}^{\\{B,Q,S\\}}\left(\kappa^{-1}\right)_{q^{\prime\prime\prime}q}\left(\frac{196}{289}\mathbf{c}_{q^{\prime}q^{\prime\prime\prime}}-\frac{1109}{2312}\mathbf{c}_{q^{\prime}}\mathbf{c}_{q^{\prime\prime\prime}}\right)\right]\right.$
$\displaystyle\left.\quad+\frac{162}{289}\mathbf{c}_{q^{\prime}}\frac{\partial}{\partial\alpha_{q^{\prime\prime}}}\left[\sum_{q^{\prime\prime\prime}}^{\\{B,Q,S\\}}\left(\kappa^{-1}\right)_{q^{\prime\prime\prime}q}\left(\mathbf{c}_{q^{\prime}q^{\prime\prime\prime}}-\mathbf{c}_{q^{\prime}}\mathbf{c}_{q^{\prime\prime\prime}}\right)\right]\right\\}\left(\mathcal{T}_{q^{\prime\prime}0}+\frac{\beta}{4}\sum_{q^{\prime\prime\prime}}^{\\{B,Q,S\\}}\mathcal{T}_{q^{\prime\prime}q^{\prime\prime\prime}}\mathbf{c}_{q^{\prime\prime\prime}}\right)\overset{\text{single}}{\rightarrow}\tau_{q}\;.$
(235) $\displaystyle\lambda_{VV}^{(q^{\prime},q)}$
$\displaystyle=\frac{9}{68\beta\sigma^{2}_{\text{tot}}P}\sum_{q^{\prime\prime}}^{\left\\{B,Q,S\right\\}}\left(\kappa^{-1}\right)_{q^{\prime\prime}q}\Bigg{(}\frac{2768}{765}\mathbf{c}_{q^{\prime\prime}q^{\prime}}-\frac{20879}{12240}\mathbf{c}_{q^{\prime\prime}}\mathbf{c}_{q^{\prime}}\Bigg{)}\overset{\text{single}}{\rightarrow}\frac{3}{5}\tau_{q}\;,$
(236) $\displaystyle\ell_{V\pi}^{(q^{\prime})}$
$\displaystyle=\frac{9}{80\sigma_{\text{tot}}P}\mathbf{c}_{q^{\prime}}\overset{\text{single}}{\rightarrow}\frac{\beta}{20}\tau_{q}\;,$
(237) $\displaystyle\tau_{V\pi}^{(q^{\prime})}$
$\displaystyle=\ell_{V\pi}^{(q^{\prime})}\overset{\text{single}}{\rightarrow}\frac{\beta}{20}\tau_{q}\;,$
(238) $\displaystyle\lambda_{V\pi}^{(q^{\prime},q)}$
$\displaystyle=\frac{28}{85}\frac{1}{\sigma_{\text{tot}}P}\left(\mathbf{c}_{qq^{\prime}}-\frac{295}{448}\mathbf{c}_{q}\mathbf{c}_{q^{\prime}}\right)\overset{\text{single}}{\rightarrow}\frac{\beta}{20}\tau_{q}\;.$
(239)

Note that the coefficient $\tau_{n\pi}$ from Denicol (2012) [21] was defined
as
$\tau_{n\pi}\overset{\text{single}}{=}\frac{1}{\epsilon+P}\tau_{V\pi}^{(q^{\prime})}$,
which then yields $\tau_{n\pi}=\frac{\beta}{80}\tau_{q}$. Therefore, in the
single-component limit we retrieve the same coefficients as in Ref. [21].

##  Appendix G Inverting the collision matrix

For the calculation of the transport coefficients, the inverse of the
collision matrix
$\mathcal{A}^{(\ell)}\equiv\left(\mathcal{A}^{(\ell)}_{ij,rn}\right)$ from Eq.
(90) must be calculated. In the tradition of Ref. [21], we provide a detailed
discussion of the derivation of the linearized collision term (see Appendix E)
and its inverse. In this section, we show that the collision matrix is
singular in the cases $\ell=0$ and $\ell=1$ due to the conservation of energy-
momentum and charge. While in the single-component system the construction of
the inverse was immediately clear [21], in the case of a multicomponent system
such a construction is not obvious.

The conservation equations (73) for the various charges imply that certain
moments of the Boltzmann equation vanish,

$\displaystyle\sum_{i\,=\,1}^{N_{\text{spec}}}q_{i}C_{i,0}$
$\displaystyle=\sum_{i\,=\,1}^{N_{\text{spec}}}q_{i}\int\mathrm{d}K_{i}k^{\nu}_{i}\partial_{\nu}f_{i,\mathbf{k}}\equiv
0\;.$ (240)

Similarly, projections of the conservation law (74) for energy and momentum
give

$\displaystyle\sum_{i\,=\,1}^{N_{\text{spec}}}C_{i,1}$
$\displaystyle=\sum_{i\,=\,1}^{N_{\text{spec}}}\int\mathrm{d}K_{i}E_{i,\mathbf{k}}k^{\nu}_{i}\partial_{\nu}f_{i,\mathbf{k}}\equiv
0,$ (241)
$\displaystyle\sum_{i\,=\,1}^{N_{\text{spec}}}C^{\langle\mu\rangle}_{i,0}$
$\displaystyle=\sum_{i\,=\,1}^{N_{\text{spec}}}\int\mathrm{d}K_{i}k^{\langle\mu\rangle}_{i}k^{\nu}_{i}\partial_{\nu}f_{i,\mathbf{k}}\equiv
0\;.$ (242)

These relations imply that a subset of row (or column) vectors of the
collision matrices $\mathcal{A}^{(\ell)}$ are linearly dependent. Since the
irreducible moments $\rho^{\mu_{1}\cdots\mu_{\ell}}_{i,r}$ are in principle
independent of each other, Eq. (90) implies that the following linear
combinations of the elements of the collision matrix must vanish:

$\displaystyle 0$
$\displaystyle=\sum_{i\,=\,1}^{N_{\text{spec}}}q_{i}\mathcal{A}^{(0)}_{is,1n}\;,$
(243) $\displaystyle 0$
$\displaystyle=\sum_{i\,=\,1}^{N_{\text{spec}}}\mathcal{A}^{(0)}_{is,2n}\;,$
(244) $\displaystyle 0$
$\displaystyle=\sum_{i\,=\,1}^{N_{\text{spec}}}\mathcal{A}^{(1)}_{is,1n}\;.$
(245)

This means that the conservation laws render the collision matrix for the
scalar moments, $\mathcal{A}^{(0)}$, and the vector moments,
$\mathcal{A}^{(1)}$, singular. We note that the tensor moments ($\ell=2$) are
not affected, and thus $\mathcal{A}^{(2)}$ is in principle regular. For a
meaningful fluid-dynamical theory, an equivalent description of the above
discussed collision matrices has to be found such that they are rendered
invertible, and at the same time their microscopic information is not altered.
The linear dependence further implies that $N_{q}+4$ equations need to be
removed from the set (83), (85), and (87) of equations of motion, i.e.
$N_{q}+1$ scalar moments $\rho_{i,r}$ and one vector moment $\rho^{\mu}_{i,r}$
(three equations). In the case of the vector moments, the choice of the frame
provides a relationship between the vector moments, and it allows us to
eliminate one of them from the equations of motion entirely.

For the Landau frame, from Eq. (70) we have

$\rho^{\mu}_{N_{\text{spec}},1}=-\sum\limits_{i\,=\,1}^{N_{\text{spec}}-1}\rho^{\mu}_{i,1}\;,$
(246)

while in the Eckart frame we could impose via Eq. (138)
$\rho^{\mu}_{N_{\text{spec}},1}=-h_{q}\sum\limits_{i\,=\,1}^{N_{\text{spec}}}q_{i}\rho^{\mu}_{i,0}-\sum\limits_{i\,=\,1}^{N_{\text{spec}}-1}\rho^{\mu}_{i,1}$.
As before, we proceed in the Landau frame. With the help of Eq. (246), we can
write:

$\displaystyle\sum_{s\,=\,1}^{N_{\text{spec}}}\sum_{n\,=\,0}^{N_{1}}\mathcal{A}^{(1)}_{is,rn}\rho^{\mu}_{s,n}$
$\displaystyle=\sum_{n\,=\,0,\neq
1}^{N_{1}}\sum_{s\,=\,1}^{N_{\text{spec}}}\mathcal{A}^{(1)}_{is,r0}\rho^{\mu}_{s,0}+\sum_{s\,=\,1}^{N_{\text{spec}}-1}\mathcal{A}^{(1)}_{is,r1}\rho^{\mu}_{s,1}+\mathcal{A}^{(1)}_{iN_{\text{spec}},r1}\rho^{\mu}_{N_{\text{spec}},1}$
$\displaystyle=\sum_{n\,=\,0,\neq
1}^{N_{1}}\sum_{s\,=\,1}^{N_{\text{spec}}}\mathcal{A}^{(1)}_{is,r0}\rho^{\mu}_{s,0}+\sum_{s\,=\,1}^{N_{\text{spec}}-1}\underbrace{\left(\mathcal{A}^{(1)}_{is,r1}-\mathcal{A}^{(1)}_{iN_{\text{spec}},r1}\right)}_{(\star)}\rho^{\mu}_{s,1}\overset{!}{=}\sum_{n\,=\,0}^{N_{1}}\sum_{s\,=\,1}^{N_{\text{spec}}}\tilde{\mathcal{A}}^{(1)}_{is,rn}\rho^{\mu}_{s,n}\;.$
(247)

First, we note that the term marked with $(\star)$ vanishes for
$s=N_{\text{spec}}$. Therefore, we can extend the last sum to run up to
$s=N_{\text{spec}}$. In the last step of the above equation, we introduced the
matrix $\tilde{\mathcal{A}}^{(1)}$, which is a reduced version of the matrix
$\mathcal{A}^{(1)}$ with the row corresponding to $r=1$ and
$i=N_{\text{spec}}$ and the column corresponding to $n=1$ and
$s=N_{\text{spec}}$, respectively, removed, i.e. its elements read:

$\displaystyle\tilde{\mathcal{A}}^{(1)}_{is,rn}=\mathcal{A}^{(1)}_{is,rn}~{}~{}\text{for}~{}~{}n\neq
1\;,\quad\text{and}\quad\tilde{\mathcal{A}}^{(1)}_{is,r1}=\mathcal{A}^{(1)}_{is,r1}-\mathcal{A}^{(1)}_{iN_{\text{spec}},r1}\;.$
(248)

Then, $\tilde{\mathcal{A}}^{(1)}$ is an $(N_{\text{spec}}\cdot
N_{1}-1)$-dimensional, rectangular, regular matrix.

An explicit example for the construction of the corrected matrix is given in
Appendix F in the case of an ultrarelativistic gas mixture. We further note
that we are free to choose which irreducible moment $\rho^{\mu}_{i,1}$ we
remove from the set of equations of motion, and thus we could have chosen any
line and associated column (corresponding to that particular moment) to be
removed (e.g. $\rho^{\mu}_{\lambda,1}$ instead of
$\rho^{\mu}_{N_{\text{spec}},1}$).

Once the reduced collision matrix $\tilde{\mathcal{A}}^{(1)}$ is obtained it
can be inverted. This yields the reduced inverse
$\tilde{\mathcal{\tau}}^{(1)}$ of dimension $N_{\text{spec}}\cdot N_{1}-1$. In
order to make the inverse $\mathcal{\tau}^{(1)}$ equivalent to the one of
dimension $N_{\text{spec}}\cdot N_{1}$ introduced in Eq. (93), and in order to
keep a simple notation regarding the summations over the indices $r$ and $s$
in all equations following that definition (e.g. Eqs. (95) ff.), we add zero-
element row(s) and column(s), which correspond to the ones originally removed
from the matrix $\tilde{\mathcal{A}}^{(1)}$. For instance, in the case of the
vector moments this means that we add a zero row for $r=1$ and
$i=N_{\text{spec}}$, and a zero column for $n=1$ and $s=N_{\text{spec}}$. Due
to the zero-element row and column, the removed irreducible moment
$\rho^{\mu}_{N_{\text{spec}},1}$, even though it formally still appears in the
equations following Eq. (95), effectively does not contribute anymore.

The procedure is analogous for the removal of the $N_{q}+1$ scalar moments.
The energy- and charge-conservation laws provide $N_{q}+1$ relations for the
linear dependence of the row vectors of matrix $\mathcal{A}^{(0)}$. These are
given by Eqs. (243) and (244). In order to remove the corresponding moments,
we impose the Landau matching conditions,

$\displaystyle\sum_{i\,=\,1}^{N_{\text{spec}}}\rho_{i,2}=0\;,\quad\text{and}\quad\sum_{i\,=\,1}^{N_{\text{spec}}}q_{i}\rho_{i,1}=0\;.$
(249)

The corrected matrix $\tilde{\mathcal{A}}^{(0)}$ is then an
$(N_{\text{spec}}\cdot N_{0}-N_{q}-1)$-dimensional, rectangular, regular
matrix and its elements can be written as:

$\displaystyle\tilde{\mathcal{A}}^{(0)}_{is,rn}=\mathcal{A}^{(0)}_{is,rn}~{}~{}\text{for}~{}~{}n\neq
1,2~{},\quad\tilde{\mathcal{A}}^{(0)}_{is,r1}=\mathcal{A}^{(0)}_{is,r1}-\frac{q_{s}}{q_{\lambda}}\mathcal{A}^{(0)}_{i\lambda,r1},\quad\text{and}\quad\tilde{\mathcal{A}}^{(0)}_{is,r2}=\mathcal{A}^{(0)}_{is,r2}-\mathcal{A}^{(0)}_{i\lambda,r2}\;,$
(250)

where again $i\neq\lambda$ for $r=1,2$, and $s\neq\lambda$ for $n=1,2$. It
should be noted that the $N_{q}$ charge-conservation laws and the $N_{q}$
matching conditions associated to the net-charge densities only allow for the
removal of moments corresponding to species $\lambda$ with non-vanishing
charge of type $q$, $q_{\lambda}\neq 0$. This can be understood by noting that
the above relations (250) are not well-defined when $q_{\lambda}=0$. In order
to simplify the notation in this work, we construct the inverse $\tau^{(0)}$
as an $(N_{\text{spec}}\cdot N_{0})$-dimensional, rectangular matrix, which is
the inverse of $\tilde{\mathcal{A}}^{\left(0\right)}$ and contains zero
elements for the rows and columns which were removed from
$\mathcal{A}^{\left(0\right)}$.

## References

  * [1] J. E. Bernhard, J. S. Moreland and S. A. Bass, Nature Phys. 15, no.11, 1113-1117 (2019) doi:10.1038/s41567-019-0611-8
  * [2] J. Auvinen, K. J. Eskola, P. Huovinen, H. Niemi, R. Paatelainen and P. Petreczky, Phys. Rev. C 102, no.4, 044911 (2020) doi:10.1103/PhysRevC.102.044911 [arXiv:2006.12499 [nucl-th]].
  * [3] G. Nijs, W. van der Schee, U. Gürsoy and R. Snellings, Phys. Rev. C 103, no.5, 054909 (2021) doi:10.1103/PhysRevC.103.054909 [arXiv:2010.15134 [nucl-th]].
  * [4] D. Everett et al. [JETSCAPE], Phys. Rev. C 103, no.5, 054904 (2021) doi:10.1103/PhysRevC.103.054904 [arXiv:2011.01430 [hep-ph]].
  * [5] J. E. Parkkila, A. Onnerstad and D. J. Kim, Phys. Rev. C 104, no.5, 054904 (2021) doi:10.1103/PhysRevC.104.054904 [arXiv:2106.05019 [hep-ph]].
  * [6] W. Israel and J.M. Stewart, Annals Phys. 118, 341 (1979).
  * [7] H. Grad, Comm. Pure Appl. Math. 2, 331 (1949); ibid. 325 (1949).
  * [8] I. Müller, Z. Physik 198, 329 (1967); Living Rev. Relativity 2, 1 (1999); http://www.livingreviews.org/lrr-1999-1.
  * [9] M. Prakash, M. Prakash, R. Venugopalan and G. Welke, Phys. Rept. 227, 321-366 (1993) doi:10.1016/0370-1573(93)90092-R
  * [10] A. Monnai and T. Hirano, Nucl. Phys. A 847, 283-314 (2010) doi:10.1016/j.nuclphysa.2010.08.002 [arXiv:1003.3087 [nucl-th]].
  * [11] A. Monnai and T. Hirano, J. Phys. Conf. Ser. 270, 012042 (2011) doi:10.1088/1742-6596/270/1/012042 [arXiv:1008.3298 [nucl-th]].
  * [12] L. Onsager, Phys. Rev. 37, 405 (1931) doi:10.1103/physrev.37.405 .
  * [13] L. Onsager, Phys. Rev. 38, 2265 (1931) doi:10.1103/physrev.38.2265 .
  * [14] Y. Kikuchi, K. Tsumura and T. Kunihiro, Phys. Rev. C 92, no.6, 064909 (2015) doi:10.1103/PhysRevC.92.064909 [arXiv:1507.04894 [hep-ph]].
  * [15] M. Greif, J. A. Fotakis, G. S. Denicol and C. Greiner, Phys. Rev. Lett. 120, no.24, 242301 (2018) doi:10.1103/PhysRevLett.120.242301 [arXiv:1711.08680 [hep-ph]].
  * [16] J. A. Fotakis, M. Greif, C. Greiner, G. S. Denicol and H. Niemi, Phys. Rev. D 101, no.7, 076007 (2020) doi:10.1103/PhysRevD.101.076007 [arXiv:1912.09103 [hep-ph]].
  * [17] J.-B. Rose, M. Greif, J. Hammelmann, J. A. Fotakis, G. S. Denicol, H. Elfner and C. Greiner, Phys. Rev. D 101, no.11, 114028 (2020) doi:10.1103/PhysRevD.101.114028 [arXiv:2001.10606 [nucl-th]].
  * [18] J. A. Fotakis, O. Soloveva, C. Greiner, O. Kaczmarek and E. Bratkovskaya Phys. Rev. D 104, no.3, 034014 (2021) doi:10.1103/PhysRevD.104.034014 [arXiv:2102.08140 [hep-ph]].
  * [19] G. S. Denicol and H. Niemi, Nucl. Phys. A 904-905, 369c-372c (2013) doi:10.1016/j.nuclphysa.2013.02.026 [arXiv:1212.1473 [nucl-th]].
  * [20] Z. Wolff and D. Molnar, Phys. Rev. C 96, no.4, 044909 (2017) doi:10.1103/PhysRevC.96.044909 [arXiv:1611.09185 [nucl-th]].
  * [21] G.S. Denicol, H. Niemi, E. Molnar and D.H. Rischke, Phys. Rev. D 85, 114047 (2012) [Erratum-ibid. D 91, no. 3, 039902 (2015)] [arXiv:1202.4551 [nucl-th]].
  * [22] S.R. de Groot, W.A. van Leeuwen and Ch.G. van Weert, Relativistic Kinetic Theory - Principles and applications, (North Holland, Amsterdam, 1980).
  * [23] C. Cercignani and G.M. Kremer, The Relativistic Boltzmann Equation: Theory and Applications, (Birkhäuser, Basel, 2002).
  * [24] P. Chakraborty and J. I. Kapusta, Phys. Rev. C 83, 014906 (2011) doi:10.1103/PhysRevC.83.014906 [arXiv:1006.0257 [nucl-th]].
  * [25] M. Albright and J. I. Kapusta, Phys. Rev. C 93, no.1, 014903 (2016) doi:10.1103/PhysRevC.93.014903 [arXiv:1508.02696 [nucl-th]].
  * [26] F. Jüttner, Ann. Phys. 339, 856 (1911); Z. Phys. 47, 542 (1928).
  * [27] G.S. Denicol, E. Molnár, H. Niemi and D.H. Rischke, Eur. Phys. J. A 48, 170 (2012) [arXiv:1206.1554 [nucl-th]].
  * [28] L.D. Landau and E.M. Lifshitz, Fluid Dynamics, Second Edition, (Butterworth-Heinemann, Oxford, 1987).
  * [29] C. Eckart, Phys. Rev. 58, 919 (1940).
  * [30] A. El, I. Bouras, F. Lauciello, Z. Xu and C. Greiner, [arXiv:1103.4038 [hep-ph]].
  * [31] A. El, I. Bouras, C. Wesp, Z. Xu and C. Greiner, Eur. Phys. J. A 48, 166 (2012) doi:10.1140/epja/i2012-12166-6 [arXiv:1206.3465 [hep-th]].
  * [32] G. S. Denicol, H. Niemi, I. Bouras, E. Molnar, Z. Xu, D. H. Rischke and C. Greiner, Phys. Rev. D 89, no.7, 074005 (2014) doi:10.1103/PhysRevD.89.074005 [arXiv:1207.6811 [nucl-th]].
  * [33] H. Struchtrup, Physics of Fluids 16, 3921 (2004) https://doi:10.1063/1.1782751
  * [34] A. Monnai, B. Schenke and C. Shen, “Equation of state at finite densities for QCD matter in nuclear collisions,” Phys. Rev. C 100, no.2, 024907 (2019) doi:10.1103/PhysRevC.100.024907 [arXiv:1902.05095 [hep-ph]].
  * [35] J. Noronha-Hostler, P. Parotto, C. Ratti and J. M. Stafford, Phys. Rev. C 100, no.6, 064910 (2019) doi:10.1103/PhysRevC.100.064910 [arXiv:1902.06723 [hep-ph]].
  * [36] E. Molnár, H. Niemi, G. S. Denicol and D. H. Rischke, Phys. Rev. D 89, no.7, 074010 (2014) doi:10.1103/PhysRevD.89.074010 [arXiv:1308.0785 [nucl-th]].