\section{Introduction}
The pressure (or equation of state) of Quantum Chromodynamics (QCD) is the most fundamental equilibrium quantity in the study of strongly interacting matter. Knowledge of the pressure is of fundamental importance in studies of the early universe \cite{Schwarz:2003du,fodor}, heavy-ion collisions \cite{Busza:2018rrf}, and neutron stars (NSs) \cite{Lattimer:2000nx}. At sufficiently large temperatures $T$ or chemical potentials $\mu$, asymptotic freedom allows one to expand the pressure in a weak-coupling series in the strong coupling $g$. Perturbative QCD (pQCD) calculations making use of this property have been applied to compute the pressure to partial next-to-next-to-next-to-leading order (N3LO) at $T>0,\mu=0$ \cite{Kajantie:2002wa} and at $T=0,\mu>0$ \cite{Gorda:2021znl,Gorda:2021kme}, and to complete next-to-next-to-leading order (N2LO) at both $T>0,\mu>0$ \cite{Kurkela:2016was}. 

However, in many such calculations, the masses of quarks $m$ are neglected as small. While there is evidence suggesting that the effect of quark masses to the pressure are indeed small \cite{Kurkela:2009gj}, there are \emph{also} indications that they could play a role in certain systems, such as hypothetical strange quark stars \cite{Fraga:2004gz}. 
In fact, even in NSs the quark-mass effects have particularly interesting qualitative implications:  Without masses, QCD matter with equal densities of up, down, and strange quarks is both $\beta$-equilibriated and charge neutral, which are precisely the conditions experienced for cold, isolated neutron stars. Thus, in order to include leptons in $\beta$-equilibrium at high densities, the strange-quark-mass effects are necessary. 

Including quark masses in full generality within pQCD calculations is a tedious process. It has been achieved at $T=0,\mu>0$ to N2LO \cite{Kurkela:2009gj} and at $T>0,\mu>0$ to next-to-leading order (NLO) \cite{Laine:2006cp}, but even in these relatively low-order cases the resulting expressions are cumbersome and involve complicated numerical integrals: Indeed, even in the simplest leading order (LO) and NLO cases the massive, nonzero-temperature one-loop integrals cannot be expressed in terms of standard special functions and their derivatives.

Here, we choose a different approach, and examine a method of obtaining vastly simpler expressions by sacrificing only very little accuracy. We start with the observation that at the values of $T$ and $\mu$ where pQCD results are reliable, the ratio $m/(\pi T)$ or $m/\mu$ is still small for the three lightest quarks.\footnote{Note that with regards to temperature, the relevant scale is the energy of the first fermionic  Matsubara mode, which is $\pi T$.} For the more massive quarks, this approximation breaks down, and in particular one may not use this approximation when studying quark decoupling \cite{Laine:2006cp}. However, in many situations of interest, only the three lightest quarks are active.
For instance, at high densities the pQCD results are convergent for $\mu_q \gtrsim 0.9$~GeV, which is larger than the mass of the strange quark $m_s\approx 0.1$~GeV, and the charm quark is not yet active at this density \cite{Jimenez:2019kji}. A similar situation exists for temperatures $T \gtrsim 0.3$~GeV, where $m_s/(\pi T)$ is small, which approaches the regime of relevance for heavy-ion collisions. 

In this paper, we compute the pQCD pressure at nonvanishing chemical potentials, temperatures, and quark masses to $O(g^5)$ by expanding it in the quark masses, under the simplifying assumption that the masses are \emph{soft} --- here defined as $O(g\mu)$ or $O(g T)$. While there is no fundamental mechanism to connect the quark mass and coupling in this way, we observe that the mass of the strange quark is indeed a soft quantity\footnote{Strictly speaking, we demand the ratio $m_q/m_\text{E}$ to be at most order unity, where $m_{\text{E}}$ [defined below in Eq.~\eqref{eq:leading_me}] is the scale of Debye screening of gluons, but we opt to use the conventional notation for soft scales instead.} for the values of $\mu$ and $T$ discussed in the previous paragraph. Moreover, by considering a more general scaling of the quark masses in powers of the coupling, we can organize the relative importance of the mass corrections and coupling corrections within the same series expansion. Our result then brings the massive pQCD result at both $T, \mu > 0$ in line with the $m=0$ counterpart. 

\begin{figure*}[hbt]
    \centering
    \includegraphics[width=0.96\textwidth]{dpOp0_T_progression_v2.pdf}
   
    \caption{A comparison of the mass correction to the pressure $\Delta p \equiv - \Delta \Omega$ calculated with the full quark-mass dependence (dashed lines) and within our perturbative scheme (solid lines) at $T = 0$ (left panel), $T = 100$~MeV (middle panel), and $T = 200$~MeV (right panel)  for three quarks with the same chemical potential $\mu_q$. We take the up and down quarks as massless, and the one-loop running strange quark mass is used. 
    In the leftmost panel, the comparison to the full results from \Ref\cite{Fraga:2004gz} is used, while in all other panels, the full result from \Ref\cite{Laine:2006cp} plus the full mass correction to the DR term, computed in the appendices is shown. 
    All panels show the mass correction divided by the free massless pressure, and the shaded regions show the renormalization-scale variation, as explained in the main text.}
    \label{fig:pres_comparisons_prog}
\end{figure*}


\begin{figure}[hbt]
    \centering
   
    \includegraphics[width=0.4\textwidth]{dpOp0_muQ_0_v2.pdf}
    \caption{A comparison of the mass corrections to the pressure $\Delta p \equiv - \Delta \Omega$ calculated with the full quark-mass dependence (dashed lines) and within our perturbative scheme (solid lines) as a function of $T$ at $\mu_q = 0$. See Fig.~\ref{fig:pres_comparisons_prog} for more details.}
    \label{fig:pres_comparisons_mu0}
\end{figure}

The outline of our paper is as follows. 
In Sec.~\ref{sec:methodology}, we start by explaining how the assumption of soft masses allows us to treat the quark-mass term in the Dirac Lagrangian as a perturbation and expand expressions about the massless limit, as well as outline the structure of the free energy after an expansion in soft quark masses. Expanding in the masses allows us to consider significantly simpler loop integrals and obtain closed-form results in $3+1$ dimensions, which we explicitly demonstrate in Sec.~\ref{sec:calc}, giving additional details in the appendices.
We next in Sec.~\ref{sec:results} examine the effect of these mass corrections to the pressure by comparing our simple results with the exact numerical expressions, finding good agreement. We then provide some applications of our results to the astrophysically relevant regime of high density and small to moderate temperatures, evaluating quantities which all become trivial in the limit of vanishing quark masses. Finally, we conclude with a small discussion.

\begin{figure*}
    \centering
    \includegraphics[width=0.98\textwidth]{xs_T_progression.pdf}
    \caption{Left: A comparison of the strangeness fraction $n_s / n_{\rm B}$ at $T = 0$ (left panel), $T = 100$~MeV (middle panel), and $T = 200$~MeV (right panel)  for charge-neutral matter in beta equilibrium, ignoring neutrinos. Here, $n_0 = 0.16\text{ fm}^{-3}$ corresponds to nuclear saturation density. Note that ignoring the strange quark mass would set $n_s / n_{\rm B} = 1/3$ under these conditions. }
    \label{fig:xs_comparisons}
\end{figure*}


\section{Methodology} 
\label{sec:methodology}
Conceptually, our method to compute the (vacuum-subtracted) QCD grand potential $\Omega_\text{QCD}$ of matter with nonzero $\mu$, $T$, and $m$ is as follows. For simplicity in this section, let us assume only one quark flavor with mass $m$. We consider the massive QCD Lagrangian, and separate off the quark-mass part
\begin{equation}
    \mathcal{L}_{\text{QCD,} m} = \mathcal{L}_{\text{QCD,} m = 0} + m \bar{\psi} \psi.
\end{equation}
We will consider the quark-mass term in the above as a perturbation $\delta \mathcal{L}_m \equiv m \bar{\psi} \psi$ and calculate corrections to the grand potential as insertions of this mass term.
Such a calculational scheme is justified for massive particles for all but the $n = 0$ Matsubara mode of a massive boson; in particular, it is always justified for fermions as long as $m/\mu$ or $m/(\pi T)$ is small. This condition is always fulfilled in the case of the three lightest quark flavors for the region of the $T,\mu$-plane where pQCD exhibits good convergence. Thus, it is only in the vacuum subtractions where we must compute expressions unexpanded in the mass, and everywhere else we can use simpler, expanded expressions.\footnote{These vacuum subtractions will all contribute at $O(g^n m^4)$ on dimensional grounds (with the powers of the coupling present for loop corrections).} In fact, only even powers of $\delta \mathcal{L}_m$ will lead to nonzero contributions, since one insertion of $\delta \mathcal{L}_m$ adds one more fermion propagator to a fermion loop without another $\gamma^\mu$ from a new vertex, leading to a zero trace.

Let us now determine precisely what terms contribute to the grand potential of QCD matter under the conditions described above. The grand potential has contributions from three dynamical scales: the hard [$O(T,\mu)$], soft [$O(g \mu, g T)$] and ultrasoft [$O(g^2 T)$] scales. Contributions from the hard scale can be evaluated using unresummed perturbation theory, while the soft and ultrasoft contributions for $T>0$ can be treated using the dimensionally reduced (denoted ``DR'')  effective field theories of electrostatic QCD (EQCD) and magnetostatic QCD (MQCD), the latter of which does not enter until $O(g^6)$ and as such will not be needed explicitly here. Lastly, the soft contributions from $\mu > 0$ are associated with hard-thermal loop (HTL) perturbation theory. 

Gathering together the contributions from different sources, we have the following expression, wherein we must subtract out the unresummed (``naive'') contributions from each of these theories to avoid potential double counting \cite{Kurkela:2016was}
\begin{equation}
    \Omega_{\rm QCD} = \Omega^{\text{naive}}_\text{QCD} + \underbrace{\left( \Omega_\text{DR}^\text{res} - \Omega_\text{DR}^\text{naive} \right)}_{\Omega_\text{DR}^\text{corr}} + \underbrace{\left ( \Omega_\text{HTL}^\text{res} - \Omega_\text{HTL}^\text{naive} \right )}_{\Omega_\text{HTL}^\text{corr}},
\end{equation}
where here $\Omega^{\text{naive}}_\text{QCD}$ is our notation for the unresummed full-theory diagrams. These contributions have the following expansions in $g$ up to $O(g^6)$ \cite{Vuorinen:2003fs,Kurkela:2016was}\footnote{We have been notified by the authors of \cite{Kurkela:2016was} that the zero-temperature limit given in the publication has a minor error, and the first term of the second line of their Eq.~(B7) should have an overall factor of $1/48$ instead of $1/72$. This only applies to their zero-temperature results, which we have not directly used.} 
\begin{align}
    \label{eq:Omega_naive}
    \Omega_\text{QCD}^\text{naive} =& \; \Omega_\text{free} + g^2 \Omega_{2,\text{QCD}} + g^4 \Omega_{3,\text{QCD}} + O(g^6) \displaybreak[0]\\
    \label{eq:Omega_DR}
    \Omega_\text{DR}^\text{corr}/T =& \; m_{\text{E}}^3 \Omega_{1,\text{DR}} + g_{\text{E}}^2 m_{\text{E}}^2 \Omega_{2,\text{DR}} + g_{\text{E}}^4 m_{\text{E}} \Omega_{3,\text{DR}} \nonumber \\ &+ O(g^6) \displaybreak[0]\\
    \Omega_\text{HTL}^\text{corr} =& \; m_{\text{E}}^4 \Omega_{1,\text{HTL}} + O(g^6).
\end{align}
For a general description of the dimensionally reduced sector, see \Ref\cite{Kajantie:2002wa}. Above, the perturbative expansions of the EQCD parameters are given (in $d$ = 3) by
\begin{align}
    m_{\text{E}}^2 =& g^2 \alpha_\text{E4} + \frac{g^4}{(4\pi)^2} \alpha_\text{E6} + O(g^6) \\
    g_{\text{E}}^2 =& g^2 T + O(g^4),
\end{align}
where the $\alpha_{\mathrm{E}i}$ are matching parameters of EQCD, often written at $\mu=0$ but obtainable for $\mu>0$ via their general definitions \cite{Laine:2005ai}. For massive quarks, the $\Omega$ and $\alpha_\text{E}$ functions appearing above will in general become functions of this quark mass as well. Let us define 
\begin{equation}
    \Delta \Omega(m,\mu,T) = \Omega_{\rm QCD}(m,\mu,T) - \Omega_{\rm QCD}(0,\mu,T)
\end{equation} 
as the contribution from the quark masses, with the massless part in this expression being given in \Ref\cite{Kurkela:2016was}. Due to the structure of the loop integrals at nonzero $T,\mu$, the nonzero-mass expressions are expected to be analytic in the squared masses (beyond logarithms related to renormalization). We use this assumption schematically in our power-counting, and observe it \emph{a posteriori} at low orders.  Hence, after expanding for soft $m$, we have the following contributions to $\Omega$ up to $O(g^5)$. For the unresummed QCD pressure:
\begin{enumerate}
    \item $\Omega_\text{free}$ to $O(g^0 m^4)$,
    \item $g^2 \Omega_{2,\text{QCD}}$ to $O(g^2 m^2)$,
\end{enumerate}
and for the DR term:
\begin{enumerate}
    \item $m_{\text{E}}^3 \Omega_{1,\text{DR}}$ to $O(g^3 m^2)$.
\end{enumerate}
As the mass corrections from the HTL sector will not contribute until $O(g^4 m^2)$, they are not necessary here.

\section{Evaluating the massive QCD free energy to \texorpdfstring{$O(g^5)$}{O(g\^{}5)}}
\label{sec:calc}

We start by evaluating the unresummed contributions. These come in two types. First, there are medium-corrections to the one- and two-loop hard-theory diagrams, which can be calculated perturbatively in the quark mass $m$. Second, there are vacuum corrections from the quark mass, which cannot be expanded in the quark mass $m$. The latter is, by dimensional arguments, $O(m^4)$, so only the free contribution is needed here.

At this point, we may also immediately generalize the earlier discussion to $N_f$ quarks with different masses $ m_f $ and chemical potentials $ \mu_f$. Using the superscripts ``M'' and ``V'' to distinguish the matter and vacuum parts, we find the following corrections in the quark masses, leaving the details for Appendix~\ref{app:unres_QCD}:  
\begin{align}
    \label{eq:Domega_MFD}
    \Delta\Omega_\text{free}^\text{M} \simeq& -2 N_c \sum_f \left ( 
    m_f^2 \widetilde{\mathcal{I}}^0_1 
    - \frac{m_f^4}{2} \widetilde{\mathcal{I}}^0_2 
    \right ) \displaybreak[0]\\
    \label{eq:Domega_VFD}
    \Delta\Omega_\text{free}^\text{V} =& 
    -2 N_c \sum_f \int_P \ln (P^2 + m_f^2) \nonumber \\
    =& - \frac{4 N_c}{D} \left( \frac {e^{\gamma_{\rm E}} \bar{\Lambda}^2}{4 \pi} \right)^\epsilon \frac{\Gamma ( 1 - \frac{D}{2} )}{(4 \pi)^{D/2}} \sum_f (m_f^2)^{\frac{D}{2}} \nonumber \\
    \simeq& \sum_f \frac{N_c m_f^4}{(4 \pi)^2} \left[ \frac{1}{\epsilon} + \frac{3}{2} + 2 \ln \frac{\bar{\Lambda}}{m_f} + O(\epsilon) \right], \displaybreak[0]\\
    \label{eq:Domega_2}
    g^2 \Delta\Omega_\text{2,QCD} =& 
    -2 d_{\text{A}} g^2 \sum_f m_f^2 \left[ 
        \tilde{\mathcal{T}} 
        + (1 - \epsilon ) (\widetilde{\mathcal{I}}^0_1 - \mathcal{I}^0_1) \widetilde{\mathcal{I}}^0_2 
    \right ],
\end{align}
where $\gamma_{\rm E} = 0.5772\ldots$ is the Euler--Mascheroni constant and $\bar{\Lambda}$ is the $\overline{\text{MS}}$ renormalization scale. The sums are over all quark flavors, and we note that the integrals $\widetilde{\mathcal{I}}_n^k$ as well as the two-loop integral $\widetilde{\mathcal{T}}$, whose explicit expressions are given in Appendix~\ref{app:ints}, each depend on a (single) chemical potential $\mu_f$ as well as the temperature $T$. The one remaining contribution at $O(g^5)$ is the first correction to $m_{\text{E}}$ from a nonzero quark mass, or equivalently the first mass correction to the matching coefficient $\alpha_\text{E4}$.


\begin{figure*}
    \centering
    \includegraphics[width=0.98\textwidth]{mulep_T_progression_v2.pdf}
    \caption{Left: A comparison of lepton chemical potential $\mu_\text{lep}$ at $T = 0$ (left panel), $T = 100$~MeV (middle panel), and $T = 200$~MeV (right panel)  for charge-neutral matter in beta equilibrium, ignoring neutrinos. Note that ignoring the strange quark mass would set $\mu_\text{lep} = 0$ under these conditions. }
    \label{fig:mu_lep}
\end{figure*}

We leave details of this computation to Appendix~\ref{app:mass_exp_DR}, including writing down expressions that can be used to compute numerical results for any $m_f$. In the end the computation results in
\begin{equation}\label{eq:mEcorr}
    \Delta(m_{\text{E}}^3 \Omega_\text{1,DR}) 
    = \frac{2 d_{\text{A}}}{(4 \pi)^3} m_{\text{E,0}} g^2 \sum_f m_f^2, 
\end{equation}
with the leading-order $m_{\text{E,0}}^2$ given by \cite{Laine:2005ai}
\begin{equation}
    m_{\text{E,0}}^2 = (d-1) g^2 \left \{ (d-1) N_c \mathcal{I}_1^0 + \sum_f \left[ 2\widetilde{\mathcal{I}}_2^2 -\widetilde{\mathcal{I}}_1^0 \right] \right \}.
\label{eq:leading_me}
\end{equation}
Note that we have written this definition of $m_{\text{E,0}}$ in a form that is conveniently generalized to nonzero quark mass, instead of using zero-mass recurrence relations for the scalar integrals $\widetilde{\mathcal{I}}_n^k$, which are commonly applied in similar contexts.


\section{Results and Discussions}
\label{sec:results}

The sum of Eqs.~\eqref{eq:Domega_MFD}-\eqref{eq:mEcorr} constitute our unrenormalized result, with the masses $m_f$ appearing in these equations being the bare quark masses $m_{f,\text{B}}$. To carry out renormalization and cancel the divergences present in $\Delta \Omega_\text{free}^\text{V}$ in Eq.~\eqref{eq:Domega_VFD}, as well as in $\widetilde{\mathcal{T}}$ and $\widetilde{\mathcal{I}}_2^0$, which can be derived from Eqs.~\eqref{eq:Ttilde}-\eqref{eq:Itilde}, we use the standard one-loop relation 
\begin{equation}
    m_{f,\text{B}} = \left( 1 + \frac{\delta_1 g^2}{(4 \pi)^2} + \cdots \right) m_{f,\text{ren}},
\end{equation}
with $\delta_1 = -3 C_\text{F} / \epsilon$, in the above equations, which is only necessary in the terms proportional to $m_f^2$. This removes the remaining divergences, and we can then use our results with one-loop running quark masses $m_f(\bar{\Lambda})$ and two-loop running $\alpha_s(\bar{\Lambda})$. Two-loop running of the coupling is necessary since our results are accurate beyond $O(g^4)$, while one-loop running of the mass is sufficient since we are only working to $O(m^4)$. We henceforth work with the three lightest quarks only, taking $m_u = m_d = 0$, and using the one-loop running of the strange quark mass, with $m_s(2 \text{ GeV}) = 93.6$~MeV, which is the central value from \Ref\cite{Chakraborty:2014aca}, and we take $\Lambda_\text{QCD} = 378$~MeV in the two-loop running of the coupling. This quark content is relevant for the environment connecting to NS cores and the mergers of NS-NS binaries.

We first compare our renormalized results to the full results obtained without expanding in the quark mass. Since we are only calculating contributions up to $O(g^5)$ assuming soft $m$, we can take as the full result the first two terms of the unresummed $\Omega_\text{QCD}^\text{naive}$ in Eq.~\eqref{eq:Omega_naive} and the leading term in $\Omega_\text{DR}/T$ in Eq.~\eqref{eq:Omega_DR} both with the full quark-mass dependence, since including the quark mass anywhere else will only result in higher-order corrections. The unresummed contributions to this order were calculated for all $T$ in \Ref\cite{Laine:2006cp}, with the simpler $T = 0$ limit calculated in \Ref\cite{Fraga:2004gz}. The DR correction is calculated below in Appendix~\ref{app:mass_exp_DR}. The comparison with our results is shown in Fig.~\ref{fig:pres_comparisons_prog}, for quark matter at high quark chemical potential and increasing temperatures, and in Fig.~\ref{fig:pres_comparisons_mu0} for a hot quark-gluon plasma at zero quark chemical potential. In these figures, we show the usual renormalization-scale variation in $\bar{\Lambda} = X \sqrt{(2\mu_q)^2+ (0.723 \times 4 \pi T)^2}$ \cite{Kurkela:2016was,Kurkela:2009gj,Kajantie:2002wa}, with $X \in \{1/2, 1, 2\}$, with the central dashed or solid lines corresponding to $X = 1$.  As these figures illustrate, our analytic perturbative results are in excellent agreement with the full, numerical expressions at all temperatures and densities where pQCD is applicable.

We next apply our results to investigate three quantities of approximate astrophysical interest. When plotting these quantities, we show results down to the density where the calculated pQCD pressure diverges and becomes negative. 
In Fig.~\ref{fig:xs_comparisons}, we show the strangeness fraction $n_s/n_B$ as a function of baryon chemical potential ${\mu_B = \mu_u+\mu_d+\mu_s}$ for charge-neutral matter at high density in beta equilibrium at three different temperatures. Here, and below $n_i$ is the number density of the relevant particle $i$, and $n_B$ is the baryon density ${n_B = n_u + n_d + d_s}$. We take as the possible matter content the three lightest quarks $u, d, s$ and the two lightest leptons $e, \mu$, with the leptons being taken as non-interacting beyond the condition of beta equilibrium. We neglect the mass of the electron in the following, but we include the mass of the muon (and do not expand in it). We implement beta equilibrium and charge neutrality by the usual two conditions
\begin{align}
    \label{eq:beta_equil}
    \mu_d &= \mu_s = \mu_u + \mu_\text{lep}, \quad \mu_e = \mu_\mu = \mu_\text{lep} \\
    &\frac{2}{3} n_u- \frac{1}{3} n_d - \frac{1}{3} n_s - n_e - n_\mu = 0,
\end{align}
and so we explicitly ignore (i) the presence of neutrinos and neutrino trapping, and (ii) any possible modifications to the $T = 0$ equilibrium conditions that may arise from the blurring of the Fermi surfaces as the temperature is increased \cite{Alford:2018lhf}. The former condition can be easily weakened by including neutrino chemical potentials in Eq.~\eqref{eq:beta_equil}, while the latter requires a careful calculation of rates. In general, both of these conditions must be included to correctly model the full astrophysical environment of binary NS-NS mergers. We see from this figure that the strangeness fraction deviates from its high-density limit of $1/3$ by less than $2\%$ at all of the densities and temperatures shown. Additionally, the central value of the quantity increases only very slightly as the temperature increases. While quite small for these values of $T$, the effect has a physical explanation in the increase of the thermal population of the particles.

\begin{figure*}
    \centering
    \includegraphics[width=0.98\textwidth]{xlep_T_progression.pdf}
    \caption{Left: A comparison of lepton fraction $n_\text{lep} / n_B$ at $T = 0$ (left panel), $T = 10$~MeV (middle panel), and $T = 100$~MeV (right panel)  for charge-neutral matter in beta equilibrium, ignoring neutrinos. Note that ignoring the strange quark mass would set $n_\text{lep} / n_{B} = 0$ under these conditions.}
    \label{fig:xlep_comparisons}
\end{figure*}
In Fig.~\ref{fig:mu_lep}, we show the progression of the lepton chemical potential as a function of baryon density under these conditions. As expected, the lepton chemical potential increases with decreasing  baryon density, as the quark chemical potentials decrease and the difference between $d$ and $s$ quarks become more pronounced. However, like the strangeness fraction above, $\mu_\text{lep}$ is only very modestly dependent on the temperature for $T \lesssim 200$~MeV at these densities, and for $T < 100$~MeV, relevant for NS mergers, thermal corrections to these quantities are negligible at these densities.  We note here that since the ratio $m_\mu / \mu_\text{lep}$ is not small, the muon mass must indeed be fully included within these computations, as has been done.


Finally, in Fig.~\ref{fig:xlep_comparisons}, we show the lepton fraction $n_\text{lep}/n_B$, with ${n_\text{lep} = n_e + n_\mu}$ as a function of baryon density under the same conditions. In this case, we see that the quantity deviates from its high-density limit of $0$ by less than $10^{-6}$ at all densities shown at $T = 0$, and as the temperature increases it grows quite strongly. This growth can be explained by the form of the number density of a massless lepton (for simplicity of explanation):
\begin{equation}
    \label{eq:massless_lepton_density}
    n_e(\mu_e, T) = \frac{\mu_e^3}{3 \pi ^2}+\frac{\mu_e T^2}{3}.
\end{equation}
The growth with temperature stems from the second term in Eq.~\eqref{eq:massless_lepton_density}, coupled with the fact that the lepton chemical potential is approximately temperature independent. We note that the crossover between $T = 0$ and $T$-dominated behavior for the lepton fraction occurs when $\pi T \approx \mu_\text{lep}$, which occurs for $T \approx 2$~MeV for the values of the $\mu_\text{lep}$ obtained here.

We turn now to a discussion of the results and methods presented in this work.
The methodology used above has some significant advantages over fully including the quark-mass effects. 

First, the relevant loop integrals are significantly simpler: Rather than containing the full mass dependence, they are zero-mass integrals with increasing powers of momenta in the denominator. This holds true to any finite-order expansion in the masses. With this in mind, considering masses larger in comparison to the density or temperature poses no fundamental problems so long as they still scale according to the coupling. For example, assuming the masses to be of the order $O(g^{1/2}\mu, g^{1/2} T)$ instead of the usual soft scale $O(g\mu, g T)$ considered here is a straightforward extension of our results. These improvements could for instance be carried out in the context of beyond-the-standard-model theories with additional heavy particles at high temperatures. Similarly, extension to higher orders in the coupling is relatively simple, as no fundamentally new integrals are required: expanded massive corrections to $n$-loop integrals are simply massless $n$-loop integrals with denominators raised to a higher power. For example, the 1-loop integrals $\widetilde{\mathcal{I}}^k_n$ in Eq.~\eqref{eq:Itilde} have a known closed-form expression for all $n$ and $k$, and in Appendix~\ref{app:mass_exp_DR} we give results for contributions to $m_{\text{E}}^3$ in a form that immediately allows for evaluation to any order in the masses. Hence, as further perturbative corrections beyond $O(g^5)$ are calculated in the massless case, such as those following the program at $T=0$ recently advanced in Refs.~\cite{Gorda:2021kme,Gorda:2021znl}, corrections from nonzero quark masses can then be included following our approach.

Second, our results themselves can be presented in a very simple form. Rather than involving complicated numerical integrals which are often badly behaved in certain regions of the parameter space, the results can be given in a closed analytic form, or in a simple numerical approximation. For this reason, our results have great potential for the community when it comes to applications to dense systems in astrophysics and other fields: The results incorporate the effects of a nonzero temperature and nonzero quark masses, allowing the study of phenomenologically interesting problems where these effects are necessary, yet they still maintain a simple and convenient form. For instance, with these results, one can extend the matching program of Refs.~\cite{Annala:2017llu,Annala:2019puf} to constrain the NS equation of state to nonzero temperatures while including the effects of the strange quark mass.

As our numerical comparisons show, the advantages of our methodology come with only a minimal drop in accuracy compared to the full result, which we believe to be a more than acceptable tradeoff. 

\section*{Acknowledgements} We thank Aleksi Kurkela and Aleksi Vuorinen for helpful comments and suggestions.
TG was partially supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -- Project-Id 279384907 -- SFB 1245.

