\section{The renormaalization effect of local impurity states}
Starting from Eq.(6)-(8) in the main text, we can write down the functional path-integral description of the quantum magnet, and its coupling to the local impurity states. We construct the path integral using the spin coherent state  \cite{efradkan}, which, for $S=1/2$, is denoted by $|n\rangle=\mathrm{exp}[{i\theta}(\hat{z}\times\mathbf{n})\cdot\mathbf{S}]|1/2,1/2\rangle$ and $\langle n| \mathbf{S}|n\rangle=\frac{1}{2}\mathbf{n}$. The zero temperature partition function reads as $Z=\int D\mathbf{n}\int D\overline{c}Dc e^{i(S_0+S^{\prime})}$. Here, $S_0$ is the action describing the frustrated quantum magnet:
\begin{equation}\label{eqsad1}
 S_0[\mathbf{n}]=S\sum_{\mathbf{r}}S_{WZ}[\mathbf{n}(\mathbf{r})]-\int dt\sum_{\langle \mathbf{r},\mathbf{r}^{\prime}\rangle}JS^2\mathbf{n}^{\parallel}(\mathbf{r},t)\cdot\mathbf{n}^{\parallel}(\mathbf{r}^{\prime},t)-...,
\end{equation}
where the $S_{WZ}$ term is the Wess-Zumino term describing the Berry phase accumulation due to the dynamics of each spin operator on the lattice. $\mathbf{n}^{\parallel}=(n^x,n^y)$ because we consider the XY model. Here, the NN term with coupling $J$ is written explicitly and the further neighbor terms are denoted by the ellipsis.  $S^{\prime}$ depicts the local impurity states as well as its local coupling to the nearby spins, which read as,
\begin{equation}\label{eqsad2}
 S^{\prime}=\int dt[\sum_{l,\sigma}\overline{c}_{\mathbf{r}_0,l,\sigma}(i\partial_t-\epsilon_l)c_{\mathbf{r}_0,l,\sigma}-V\sum_{\mathbf{r},l,l^{\prime},\sigma,\sigma^{\prime}}c^{\dagger}_{\mathbf{r}_0,l,\sigma}\boldsymbol{\sigma}_{\sigma,\sigma^{\prime}}c_{\mathbf{r}_0,l^{\prime},\sigma^{\prime}}\cdot\mathbf{n}(\mathbf{r},t)]
 \end{equation}
 where the sum of $\mathbf{r}$ is over the sites nearest to $\mathbf{r}_0$. For simplicity, considering a bipartite lattice with sublattice A and B, and the impurity $\mathbf{r}_0$ being located at the center of the lattice bond, then the sites $\mathbf{r}$ to be summed in the second term can be denoted as $(\mathbf{r}_0,A)$ and $(\mathbf{r}_0,B)$. Now we can integrate out the local quantum dot operators in $S^{\prime}$, which generates the following renormalization to the quantum magnet $S_0$ as,
\begin{equation}\label{eqsad3}
 \Delta S=-i\mathrm{ln}\int D\overline{c}Dce^{i\int^{\infty}_0dt\overline{c}G^{-1}[\mathbf{n}]c},
 \end{equation}
 where we introduced the Grassman field in the basis $c=[c_{\mathbf{r}_0,l,\uparrow},c_{\mathbf{r}_0,l,\downarrow}]^T$ and $\overline{c}=[\overline{c}_{\mathbf{r}_0,l,\uparrow},\overline{c}_{\mathbf{r}_0,l,\downarrow}]$. We defined the inverse of the Green's function (GF), $G^{-1}[\mathbf{n}]=G^{-1}_0-\Sigma[\mathbf{n}]$, where $G_0$, the GF of the local quantum dot, reads in frequency and orbital representation as, $G_{0,l,l^{\prime}}=[\delta_{l,l^{\prime}}/(\omega-\epsilon_l)]\sigma^0$, and $\Sigma_{l,l^{\prime}}[\mathbf{n}]=V(\widetilde{n}^z\sigma^z+\widetilde{n}^-\sigma^++\widetilde{n}^+\sigma^-)$, with $\sigma^{\pm}=\sigma^x\pm i\sigma^y$ and $\widetilde{n}^z=\sum_{\mathbf{r}}n^z(\mathbf{r},\omega)$, $\widetilde{n}^{\pm}=\sum_{\mathbf{r}}n^x(\mathbf{r},\omega)\pm in^y(\mathbf{r},\omega)$. Then, we integrate over the Grassman field $c$, $\overline{c}$, and make expansion the with respect to $\Sigma$, which is proportional to $V$ with $V\ll J,|\epsilon_l|$.  To the leading
 order of correction (second order), one obtains:
 \begin{equation}\label{eqsad4}
\Delta S=-V^2\int \frac{d\nu}{2\pi}\Pi(\nu)\sum_{\langle \mathbf{r},\mathbf{r}^{\prime}\rangle}\mathbf{n}(\mathbf{r},\nu)\cdot\mathbf{n}(\mathbf{r}^{\prime},-\nu)
  \end{equation}
where the sum of the NN sites is automatically generated.
$\Pi(\nu)$ is the polarization function due to the fermionic loop of quantum impurity states.  Since we are only interested in the low-energy window in consistent with low-energy CS fermion description of the Dirac QSL, we take the low-frequency approximation of the polarization, which well preserves the low-energy physics for $\nu\ll |\epsilon_l|$. The polarization function is then derived as,
 \begin{equation}\label{eqsadnew4}
\Pi(0)=i\sum_{l,l^{\prime}}\int\frac{d\omega}{2\pi}\frac{1}{(\omega-\epsilon_l e^{-i0^+})(\nu+\omega-\epsilon_{l^{\prime}}e^{-i0^+})}=-\sum_{l,l^{\prime}}\frac{n_F(\epsilon_l)-n_F(\epsilon_{l^{\prime}})}{\epsilon_l-\epsilon_{l^{\prime}}},
  \end{equation}
where $n_F(x)$ is the Fermi-Dirac distribution function. It is obvious that $\Pi(0)>0$ as long as there exists orbitals above and below the Fermi energy. The total renormalized action for the quantum magnet is then obtained as,
\begin{equation}\label{eqsad5}
 S_0[\mathbf{n}]=S\sum_{\mathbf{r}}S_{WZ}[\mathbf{n}(\mathbf{r})]-\int^{\infty}_0dt\sum_{\langle \mathbf{r},\mathbf{r}^{\prime}\rangle}JS^2\mathbf{n}^{\parallel}(\mathbf{r},t)\cdot\mathbf{n}^{\parallel}(\mathbf{r}^{\prime},t)-...-V^2\Pi(0)\int dt\sum_{\langle \mathbf{r},\mathbf{r}^{\prime}\rangle}\mathbf{n}(\mathbf{r},t)\cdot\mathbf{n}(\mathbf{r}^{\prime},t).
\end{equation}
After returning back to the Hamiltonian formalism with respect to the spin operators, the system is then described by:

\begin{equation}\label{eqsad6}
H=\sum_{\mathbf{r},\mathbf{r}^{\prime}}J_{\mathbf{r},\mathbf{r}^{\prime}}(S^x_{\mathbf{r}}S^x_{\mathbf{r}^{\prime}}+S^y_{\mathbf{r}}S^y_{\mathbf{r}^{\prime}})+J^{\prime}\sum_{\langle \mathbf{r},\mathbf{r}^{\prime}\rangle}\mathbf{S}_{\mathbf{r}}\cdot\mathbf{S}_{\mathbf{r}^{\prime}}.
\end{equation}
The quantum impurity states therefore renormalize the frustrated magnets by locally generating a Heisenberg-type exchange coupling of the spin operators with $J^{\prime}=V^2\Pi(0)$. This leads to an additional Ising interaction that is absent in the original XY model. As mentioned, we are considering bipartite lattice with A,B sublattice, thereby the induced Heisenberg term is cast into:
\begin{equation}\label{eqsad7}
H^{\prime\prime}=J^{\prime}\sum_{\langle \mathbf{r},\mathbf{r}^{\prime}\rangle}\mathbf{S}_{\mathbf{r}_0,A}\cdot\mathbf{S}_{\mathbf{r}_0,B}
\end{equation}
where $J^{\prime}\ll J_{\mathbf{r},\mathbf{r}^{\prime}}$ is satisfied since $V\ll  J_{\mathbf{r},\mathbf{r}^{\prime}}$. This term, describing the major effect of the local impurity states, is studied in detail in the remaining part of the supplemental material. In fact, this term is also indicated by the CS fermion description of the Dirac QSLs. In the CS fermion language, the $\mathrm{U}(1)$ gauge field is formed by combination of local string operators defined on neighbor sites \cite{ruia,tsedrakyana,ruib}. The $\mathrm{U}(1)$ gauge field is separated into a non-fluctuating and fluctuating sector \cite{ruia,ruib}, where the latter generates either the long-range spin orderings or the deconfined phases \cite{ruic}.  As required by the CS fermionization method, the non-fluctuating gauge field has to emerge from at least two nearby sites, therefore the local perturbation of the impurity should at least involve two spins on different sublattices, as is shown explicitly by Eq.\eqref{eqsad7} and the above derivations. Noting that the (emergent) sublattice is an indispensable degrees of freedom in the CS fermion description of Dirac spin-liquid, the above derivation is not unique to the studied model but is general, and is applicable when  the quantum impurity states are coupled to a Dirac spin-liquid.

\section{The impurity effect on the Dirac CS fermions}
We assume the formation and stability of a Dirac QSL, so that the first term of Eq.\eqref{eqsad6} is described by deconfined CS fermions as discussed in the main text.   Now we study $H^{\prime\prime}$, for which, we separate explicitly the Ising term and the XY term as,  $H^{\prime\prime}=J^{\prime}_zS^z_{\mathbf{r}_0,A}\cdot S^z_{\mathbf{r}_0,B}+J^{\prime}_{\parallel}(S^x_{\mathbf{r}_0,A}\cdot S^x_{\mathbf{r}_0,B}+S^y_{\mathbf{r}_0,A}\cdot S^y_{\mathbf{r}_0,B})$. Then, it should be noted that the quantum spins at $\mathbf{r}_0$ are also interacting with spins with neighbor sites through $H_0$, from which the relevant term can be extracted  as, $H_{hyb}=\sum_{\mathbf{r}^{\prime}}[J_{\mathbf{r}^{\prime}}(S^x_{\mathbf{r}_0,A}S^x_{\mathbf{r}^{\prime},B}+S^y_{\mathbf{r}_0,A}S^y_{\mathbf{r}^{\prime},B})+J_{\mathbf{r}^{\prime}}(S^x_{\mathbf{r}_0,B}S^x_{\mathbf{r}^{\prime},A}+S^y_{\mathbf{r}_0,B}S^y_{\mathbf{r}^{\prime},A})]$, where $J_{\mathbf{r}^{\prime}}$ accounts for the major contribution from the nearest neighbor exchange coupling. The bath is described by the rest terms, which, in thermal dynamic limit, still stabilizes a Dirac QSL in the long-wave length limit.  Since $J^{\prime}_z$ and $J^{\prime}_{\parallel}$ are much weaker than $J_{\mathbf{r},\mathbf{r}^{\prime}}$, they are treated as local perturbation to the Dirac CS fermions. Then impurity effect is then cast into the Hamiltonian $H^{\prime}=H^{\prime\prime}+H_{hyb}$.  We note that $J^{\prime}_{\parallel}$ terms are weak modifications and can be absorbed into the bath terms, whereas the  Ising interaction $J^{\prime}_z$ term is absent in the XY model.  Then, after representing the quantum spin using CS fermions, $H^{\prime}$ is transformed into:
\begin{equation}\label{eqs1}
  H^{\prime}=J^{\prime}_z(1/2-f^{\dagger}_{\mathbf{r}_0,A}f_{\mathbf{r}_0,A})(1/2-f^{\dagger}_{\mathbf{r}_0,B}f_{\mathbf{r}_0,B})
  +\sum_{\mathbf{r}}J_{\mathbf{r}}(f^{\dagger}_{\mathbf{r}_0,A}e^{iA_{\mathbf{r}_0,\mathbf{r}}}f_{\mathbf{r},B}+f^{\dagger}_{\mathbf{r}_0,B}e^{iA_{\mathbf{r}_0,\mathbf{r}}}f_{\mathbf{r},A}+h.c.),
\end{equation}
where $A_{\mathbf{r}_0,\mathbf{r}}$ is the non-fluctuating lattice gauge field arising from combination of the nearby string operators, which can be taken simply as complex phases and absorbed into $J_{\mathbf{r}}$. let us rename the first and second term in Eq.\eqref{eqs1} as $H_{imp}$ and $H_{hyb}$, respectively, and introduce a different notation for the CS fermions at $\mathbf{r}_0$ as, $c_{\mathbf{r}_0,\overline{\alpha}}=f_{\mathbf{r}_0,\alpha}$, where $\alpha=A,B$ and $\overline{\alpha}$ denotes the opposite sublattice index of $\alpha$. Then, $H^{\prime}$ is cast into,
\begin{eqnarray}
  H_{imp} &=& \sum_{\alpha} \epsilon_f c^{\dagger}_{\mathbf{r}_0,\alpha}c_{\mathbf{r}_0,\alpha}+\frac{U}{2}\sum_{\alpha}c^{\dagger}_{\mathbf{r}_0,\alpha}c_{\mathbf{r}_0,\alpha}c^{\dagger}_{\mathbf{r}_0,\overline{\alpha}}c_{\mathbf{r}_0,\overline{\alpha}}, \\
  H_{hyb} &=& \sum_{\mathbf{r}\alpha}J_{\mathbf{r}}(c^{\dagger}_{\mathbf{r}_0,\alpha}f_{\mathbf{r},\alpha}+h.c.),
\end{eqnarray}
where $\epsilon_f=-U/2=-J^{\prime}_z/2$. Therefore the local perturbation is transformed to a symmetric Anderson impurity with the local effective CS fermions subjected to a Hubbard interaction, which are coupled to the nearby CS fermions of the bath via hopping terms. After a Schrieffer-Wolf transformation, $H^{\prime}$ can be further written in the Kondo regime simply as the exchange coupling to an effective local spin-half impurity $\mathbf{S}_{imp}$ formed by the local CS fermions at $\mathbf{r}_0$, which reads as,
\begin{equation}\label{eqs2}
  H^{\prime}=\sum_{\mathbf{r}}\lambda({\mathbf{r}})f^{\dagger}_{\mathbf{r},\alpha}\boldsymbol{\tau}_{\alpha\beta}f_{\mathbf{r},\beta}\cdot\mathbf{S}_{imp},
\end{equation}
where $\boldsymbol{\tau}$ is the Pauli matrix defined in the pseudospin (sublattice) space. To be more general, we allow $\mathbf{r}$-dependence of the exchange coupling $\lambda({\mathbf{r}})$, describing the scattering potential to the effective quantum impurity. We have shown that the effect of the local quantum states on the frustrated magnet  can be cast into an effective Kondo-exchange in the pseudospin space. This leads to Eq.(9) of the main text.

\section{mapping to the impurity model coupled to 1D soft fermions}
In the main text, we have shown the general procedure how to map from the impurity model in 2D Dirac fermions to that in 1D soft modes. In this section, we illustrate the detailed transformations using the honeycomb lattice model as the example.

After projection to long-wave regime near the Dirac CS valleys $\mathbf{K}$ and $\mathbf{K}^{\prime}$, the bath is described as
\begin{equation}\label{eqs3}
  H_0=v_F \sum_{\mathbf{k}}f^{(+)\dagger}_{\mathbf{k},\alpha}(\boldsymbol{\tau}_{\alpha\beta}\cdot\mathbf{k}-\mu)f^{(+)}_{\mathbf{k},\beta}
  +v_F\sum_{\mathbf{k}}f^{(-)\dagger}_{\mathbf{k},\alpha}(-\boldsymbol{\tau}^{\mathrm{T}}_{\alpha\beta}\cdot\mathbf{k}-\mu)f^{(-)}_{\mathbf{k},\beta},
\end{equation}
with valley-dependent PSML. $H^{\prime}$ is projected into the long-wave regime as,
\begin{equation}\label{eqs4}
H^{\prime}_p=g_d\sum_{\mathbf{k},\mathbf{k}^{\prime},a=\pm}f^{(a)\dagger}_{\mathbf{k},\alpha}\boldsymbol{\tau}_{\alpha\beta}\cdot \mathbf{S}_{imp},f^{(a)}_{\mathbf{k}^{\prime},\beta}
+g_t\sum_{\mathbf{k},\mathbf{k}^{\prime},a=\pm}f^{(a)\dagger}_{\mathbf{k},\alpha}\boldsymbol{\tau}_{\alpha\beta}\cdot \mathbf{S}_{imp},f^{(\overline{a})}_{\mathbf{k}^{\prime},\beta}.
\end{equation}
where $a$ denotes the two valleys. One firstly make a unitary transformation to diagonalize CS fermions at each Dirac valley, $H^{\prime}_p$ transform accordingly under the unitary transformation, leading to
\begin{equation}\label{eqs4}
 H^{\prime}_p=g_d\sum_{\mathbf{k},\mathbf{k}^{\prime},a=\pm} c^{(a)\dagger}_{\mathbf{k}} U^{(a)}(\theta_{\mathbf{k}})\tau^iU^{(a)\dagger}(\theta_{\mathbf{k}^{\prime}}) c^{(a)}_{\mathbf{k}^{\prime}}S^i_{imp}
+g_t\sum_{\mathbf{k},\mathbf{k}^{\prime},a=\pm} c^{(a)\dagger}_{\mathbf{k}} U^{(a)}(\theta_{\mathbf{k}})\tau^iU^{(\overline{a})\dagger}(\theta_{\mathbf{k}^{\prime}}) c^{(\overline{a})}_{\mathbf{k}^{\prime}}S^i_{imp},
\end{equation}
where $c^{(a)}_{\mathbf{k}}$ is the transformed spinor in band (sublattice) space at valley $a$, $U^{(a)}(\theta_{\mathbf{k}})$ the unitary rotation matrix applied for fermions at valley $a$ which is only dependent on the angle of momentum $\theta_{\mathbf{k}}$. Then, utilizing the rotational symmetry of the impurity scattering, we transform the fermions to the orbital angular momentum partial waves using $c^{(a)}_{\mathbf{k}}=\sum_l e^{il\theta}c^{(a)}_{l,k}/\sqrt{2\pi k}$, where $l$ is the partial wave index. After insertion of the specific form of the unitary rotation matrix $U^{(a)}(\theta_\mathbf{k})$, the integral over the polar angle automatically picks up several different partial waves $l$, generating the following coupling as,
\begin{equation}\label{eqs5}
 H^{\prime}_p=g_d\sum_{a} \int dkdk^{\prime}\sqrt{kk^{\prime}}
 c^{(a)\dagger}_{l,k} U^{(a)}(l)\boldsymbol{\tau}\cdot \mathbf{S}_{imp} U^{(a)\dagger}_l c^{(a)}_{l,k^{\prime}}
+g_t\sum_{a} \int dkdk^{\prime}\sqrt{kk^{\prime}}
 c^{(a)\dagger}_{l,k} U^{(a)}(l)\boldsymbol{\tau}\cdot \mathbf{S}_{imp} U^{(\overline{a})\dagger}_l c^{(\overline{a})}_{l,k^{\prime}},
\end{equation}
where constants have been absorbed into the tuning parameter $g_d$ and $g_t$, the sum over repeated notations such as $l$ is implicit. $U^{(a)}_l$ is the rotation matrix again transformed to the angular orbital momentum space, whose components are delta functions that select the channel $l$ relevant to the impurity, i.e.,
\begin{equation}\label{eqs6}
U^{(a)}_l=\frac{1}{\sqrt{2}}\left(
\begin{array}{cc}
\delta_{l,0} & a\delta_{l,-a}\\
\delta_{l,0} & -a\delta_{l,-a}\\
\end{array}
\right),
\end{equation}
where $a=\pm$ denotes the two valleys. Eq.\eqref{eqs5} implies that the impurity is coupled to an effective CS fermions $d^{(a)}_{\mathbf{k}}=\sum_l U^{(a)\dagger}(\theta_{\mathbf{k}})c^{(a)}_{l,k}$, which is combinations of  1D CS fermions with different index $l$ for different valleys. $l=0,-1$ are coupled to the impurity at $\mathbf{K}$ whle $l=0,1$ are involved at $\mathbf{K}^{\prime}$ valley. Therefore, the impurity only picks up these relevant $l$ channels. Since the bath, after rotation to the angular orbital momentum space, enjoy independent $l$ components with $l$ being good quantum number due to the rotational invariance of the problem, we can select from the bath these relevant channels, leading to,
\begin{equation}\label{eqs7}
H_0=\sum_{l=-1,0}\int^{\infty}_0 dk(\epsilon_{k,\alpha}-\mu)c^{(+)\dagger}_{k,l,\alpha}c^{(+)}_{k,l,\alpha}+\sum_{l=0,1}\int^{\infty}_0 dk(\epsilon_{k,\alpha}-\mu)c^{(-)\dagger}_{k,l,\alpha}c^{(-)}_{k,l,\alpha},
\end{equation}
where $\epsilon_{k,\alpha}=\alpha v_Fk$. It is convenient to introduce the energy representation for the impurity problem \cite{zitko}, and define the effective CS fermions with combination of operators for the conduction and valence Dirac band as,
\begin{eqnarray}
  d^{(+)}_{\epsilon} &=& \frac{1}{\sqrt{2}}[c^{(+)}_{\epsilon,0,+}\theta(\epsilon)+c^{(+)}_{\epsilon,0,-}\theta(-\epsilon), c^{(+)}_{\epsilon,-1,+}\theta(\epsilon)-c^{(+)}_{\epsilon,-1,-}\theta(-\epsilon)]^{\mathrm{T}}, \\
  d^{(-)}_{\epsilon}&=& \frac{1}{\sqrt{2}}[c^{(-)}_{\epsilon,0,+}\theta(\epsilon)+c^{(-)}_{\epsilon,0,-}\theta(-\epsilon), -c^{(-)}_{\epsilon,1,+}\theta(\epsilon)+c^{(-)}_{\epsilon,1,-}\theta(-\epsilon)]^{\mathrm{T}}.
\end{eqnarray}
Using $d^{(a)}_{\epsilon}$, $H_0$ is cast into a simple form as,
\begin{equation}\label{eqs8}
H_0=\sum_{a,\sigma}\int^{\infty}_{-\infty} d\epsilon(\epsilon-\mu)d^{(a)\dagger}_{\epsilon,\sigma}d^{(a)}_{\epsilon,\sigma},
\end{equation}
where $v_F$ is set to 1. $d^{(a)}_{\epsilon,\sigma=1,2}$ are the two entries of the spinor defined in Eq.(19) and (20). Accordingly, the hybridization term $H^{\prime}_p$ is reduced to the following form as,
\begin{equation}\label{eqs9}
 H^{\prime}_p=g_d\sum_{a} \int^{+\infty}_{-\infty} d\epsilon d\epsilon^{\prime} [\rho(\epsilon)\rho(\epsilon^{\prime})]^{1/2}
 d^{(a)\dagger}_{\epsilon} \boldsymbol{\tau}\cdot \mathbf{S}_{imp} d^{(a)}_{\epsilon^{\prime}}
+g_t\sum_{a} \int^{+\infty}_{-\infty} d\epsilon d\epsilon^{\prime} [\rho(\epsilon)\rho(\epsilon^{\prime})]^{1/2}
 d^{(a)\dagger}_{\epsilon} \boldsymbol{\tau}\cdot \mathbf{S}_{imp} d^{(\overline{a})}_{\epsilon^{\prime}},
\end{equation}
where $\rho(\epsilon)=|\epsilon|/2\pi v^2_F$ is the density of states of Dirac CS fermions, leading to a pseudogap in the above hybridizations. Detailed studies on the pseudogapped cases have shown that the strong coupling fixed points at zero temperature are not modified by approximating the density of states by that of the Fermi energy \cite{Ruid}, as long as $\mu\neq0$. With this approximation, one can absorb the density of states into the couplings and rename the fermionic field as $\psi^a$. This leads to an Kondo-exchange model coupled to 1D chiral soft modes, in consistent with the general form, i.e., Eq.(3),(4) in the main text.

\section{derivation of the decoupled Wess-Zumino-Witten CFT using non-Abelian gauge invariance}
The infrared fixed point of the reduced model (impurity coupled to the 1D soft modes) is described by a Wess-Zumino-Witten (WZW) CFT. We now show in this section that the underlying CFT has a decoupled multichannel structure and can be derived simply from the principle of gauge invariance.
The following contents are separated into three steps including the derivation of Ward identities from the non-Abelian gauge symmetry, the chiral symmetry, and the deduction of the exact functional free energy.
\subsection{Ward Identity from non-abelian gauge transformations}

From the mapped 1D model of soft modes, we can start from a Dirac field in a representation $r$ of a Lie group $G$ coupled with a given gauge field $A$,namely $\mathcal{L}=\bar{\psi} (i\slashed D-m)\psi$ with $D_\mu=\partial_\mu-igA_\mu$.
We define the free energy $W$ as
\begin{equation}
e^{-iW[A]}=Z[A]=\int  \mathcal{D}\psi\mathcal{D}\bar{\psi} ~e^{iS}.
\end{equation}
The classical theory is invariant under the gauge transformations, $\psi \rightarrow U\psi$,$\bar{\psi} \rightarrow \bar{\psi}U^{-1}$, and $A_\mu \rightarrow A^{U}=U A_\mu U^{-1}+\frac{i}{g} U\partial_\mu U^{-1}$, whose infinitesimal version is
$ \psi \rightarrow (1+i\alpha)\psi$, $\bar{\psi} \rightarrow \bar{\psi} (1-i\alpha)$, and
$A_\mu \rightarrow A_\mu +\frac{1}{g}\mathcal{D}_\mu \alpha$. Assuming that the functional measurement $\mathcal{D}\bar{\psi} \mathcal{D}\psi$ is also gauge invariant, the free energy satisfies
\begin{equation}
W[A]=W[A^U]+2\pi n[U]
\end{equation}
with $n$ being an integer determined by $U$, which vanishes for infinitesimal transformations. Accordingly,
\begin{equation}
\begin{split}
0 &= W[A_\mu +\mathcal{D}_\mu \alpha]-W[A_\mu]= \int dx \frac{\delta W}{\delta A_\mu^a}(\mathcal{D}_\mu \alpha)^a= -\int dx ~ \mathrm{tr}\alpha\mathcal{D}_\mu\frac{\delta W}{\delta A_\mu},
\end{split}
\end{equation}
which implies
\begin{equation}
\mathcal{D}_\mu\frac{\delta W}{\delta A_\mu}=0.
\end{equation}
The variation of the free energy to the gauge field is then calculated as
\begin{equation}
\begin{split}
-i\frac{\delta W}{\delta A_\mu}&=\frac{1}{Z}\frac{\delta Z}{\delta A_\mu}=\frac{1}{Z} \int \mathcal{D}\psi\mathcal{D}\bar{\psi} ~ i\frac{\delta S}{\delta A_\mu}e^{iS}= i\frac{1}{Z} \int \mathcal{D}\psi\mathcal{D}\bar{\psi} ~ J^\mu e^{iS}= i\langle J^\mu \rangle_{A}.
\end{split}
\end{equation}
Thus we prove the Ward identity
\begin{equation}
\mathcal{D}_{\mu}\langle J^\mu \rangle=0.
\end{equation}

\subsection{Ward Identity from chiral invariance}
Now we study the chiral gauge transformations given by $\psi \rightarrow \psi'=(1+i\alpha \gamma^5)\psi$m
$\bar{\psi} \rightarrow \bar{\psi}'=\bar{\psi}(1+i\alpha \gamma^5)$ and
$A_\mu \rightarrow A_\mu^\prime =A_\mu + \mathcal{D}_\mu \alpha\gamma^5$.
It is straightforward to check that the classical theory is invariant under the gauge transformations. However the functional measurement does not respect the transformations, leading to a Jacobian determinant $\mathcal{J}$. Thus the Ward identity should be modified because of $W[A^\prime]\ne W[A]$. From the partition function, we obtain,
\begin{equation}
\begin{split}
Z[A]&=\int \mathcal{D}\psi\mathcal{D}\bar{\psi}~ e^{iS[\psi,\bar{\psi},A]}=\int \mathcal{D}\psi'\mathcal{D}\bar{\psi'}~ \mathcal{J}~ e^{iS[\psi',\bar{\psi}',A']}= Z[A']+\int \mathcal{D}\psi\mathcal{D}\bar{\psi}\int dx~\frac{\delta \mathcal{J}}{\delta \alpha^a}|_{\alpha(x)=0}\alpha^a(x)e^{iS[\psi,\bar{\psi},A']},
\end{split}
\end{equation}
leading to
\begin{eqnarray}
\frac{Z[A']-Z[A]}{Z[A]}=-\langle \int dx~\frac{\delta \mathcal{J}}{\delta \alpha^a}|_{\alpha(x)=0}\alpha^a(x)\rangle.
\end{eqnarray}
Besides, we have
\begin{eqnarray}
\frac{Z[A']-Z[A]}{Z[A]}=-ig \langle \int dx (\mathcal{D}_\mu J^{5\mu})^a \alpha^a\rangle.
\end{eqnarray}
Thus, conservation equation for the axial current is obtained as,
\begin{equation}
\mathcal{D}_\mu J^{5\mu}=-i\frac{1}{g} \frac{\delta\mathcal{J}}{\delta\alpha}|_{\alpha=0}.
\end{equation}
The remaining task is then to evaluate the Jacobian determinant.  This can be readily done using the method developed by Fujikawa, which is also utilized in a similar situation of 3+1D with the chiral anomaly. A straightforward calculation in 1+1D then generates the Ward identity for the axial current as,
\begin{equation}
\mathcal{D}_\mu J^{5\mu}=-\frac{C(r)}{2\pi}\epsilon^{\mu\nu}F_{\mu\nu}.
\end{equation}
where in the derivation we have defined the Dirac matrices $\gamma^0=\sigma^2$, $\gamma^1=i\sigma^1$ and $\gamma^3=\gamma^0\gamma^1=\sigma^3$ and used  $\mathrm{tr}(t^at^b)=C(r)\delta^{ab}$ with $C(r)$ a constant for eacg representation $r$ with $t^a$ the representation matrix.

\subsection{The exact functional determinant in two dimensions}
Noting  that there exists a unique relation only in 1+1D dimensions,
$\gamma^\mu \gamma^3=-\epsilon^{\mu\nu}\gamma_\nu$, which enables us to rewrite the chiral current as
$J^{3\mu}=-\epsilon^{\mu\nu}J_\nu$.
Therefore, the two Ward identities derived above are collected into a united form of the CS fermion current as,
\begin{eqnarray}
\mathcal{D}_\mu J^\mu&=& 0\\
\epsilon^{\mu\nu} \mathcal{D}_\mu J_\nu &=& \frac{C(r)}{2\pi}\epsilon^{\mu\nu}F_{\mu\nu}.
\end{eqnarray}
Now the uniqueness of dimension two, compared with higher dimensions, lies in that the current $J^\mu$ is completely determined by the two Ward identities. Before solving the equations we first introduce the chiral coordinates,
$x^{+}=x^0+x^1$, $x^{-}=x^0-x^1$.
In the chiral coordinates the metric $\eta$ and total anti-symmetric tensor $\epsilon$ are represented, respectively, as
\begin{equation}
\eta_{\mu\nu}=\begin{pmatrix}
0 & \frac{1}{2}\\
\frac{1}{2} & 0
\end{pmatrix}, \quad \epsilon^{\mu\nu}=\begin{pmatrix}
0 & 2\\
-2 & 0
\end{pmatrix}.
\end{equation}
Accordingly we define $J^+=J^0+J^1$, $J^-=J^0-J^1$,  and $A^+=A^0+A^1$, $ A^-=A^0-A^1$. In these notations the two identities can be cast into the following symmetric form,
\begin{eqnarray}
\partial_{+} J_{-}-i[A_{+},J_{-}]&=&\frac{C(r)}{2\pi}F_{+-}\\
\partial_{-} J_{+}-i[A_{-},J_{+}]&=&\frac{C(r)}{2\pi}F_{-+}.
\end{eqnarray}
To obtain the explicit form of the solution, we introduce the expression for the gauge fields,
$A_{+}=ig^{-1}\partial_{+}g$, $A_{-}=ih^{-1}\partial_{-}h$.
with $g$ and $h$ being fields of group elements in $G$. Then it is straightforward to check that
\begin{eqnarray}
J_{+}&=&\frac{C(r)}{2\pi}(ig^{-1}\partial_{+}g-ih^{-1}\partial_{+}h)\\
J_{-}&=&\frac{C(r)}{2\pi}(ih^{-1}\partial_{-}h-ig^{-1}\partial_{-}g),
\end{eqnarray}
are the solutions of the equations.

As promised we shall work out an explicit expression of the free energy $W[A]$, which is gauge invariant, using the fields $g$ and $h$. The field $t(x)\in G$ gives the gauge transformations,
\begin{eqnarray}
A_{+}=ig^{-1}\partial_{+}g &\longrightarrow& itg^{-1}\partial_{+}gt^{-1}+it\partial_{+}t^{-1}=i(gt^{-1})^{-1}\partial_{+}(gt^{-1})\\
A_{-}=ih^{-1}\partial_{-}h &\longrightarrow& ith^{-1}\partial_{-}ht^{-1}+it\partial_{-}t^{-1}=i(ht^{-1})^{-1}\partial_{-}(ht^{-1}),
\end{eqnarray}
which are translated to $g$ and $h$ as
\begin{equation}
(g,h)\longrightarrow (g,h)t^{-1}.
\end{equation}
The gauge invariance of $W[A]$ is now expressed as
\begin{equation}
W[g,h]=W[gt^{-1},ht^{-1}],
\end{equation}
for any field $t(x)$. So it is sufficient to work with the gauge $A_{-}=0$, or equivalently $h$ constant.
\begin{equation}
\delta W=-\frac{1}{\pi}\int dx~\mathrm{tr}(g^{-1}\partial_{-}g~\delta (g^{-1}\partial_{+}g))=-\frac{1}{\pi}\int dx~ \left(\mathrm{tr}(\partial_{+}\partial_{-}g^{-1}\delta g)-\mathrm{tr}(g^{-1}\partial_{+}g~g^{-1}\partial_{-}g~g^{-1}\delta g)\right)
\end{equation}
Noting that
\begin{equation}
\begin{split}
\delta \int dx~\mathrm{tr\partial_{-}g^{-1}\partial_{+}g}&=-2\int dx~ \mathrm{tr}(\partial_{+}\partial_{-}g^{-1})\delta g + \int dx~ \mathrm{tr}(g^{-1}\partial_{+}g~g^{-1}\partial_{-}g~g^{-1}\delta g)\nonumber\\
&+\int dx~ \mathrm{tr}(g^{-1}\partial_{-}g~g^{-1}\partial_{+}g~g^{-1}\delta g),
\end{split}
\end{equation}
 we have
\begin{equation}
\delta W=-\frac{1}{8\pi}\delta\int dx~ \mathrm{tr}(g^{-1}\partial^\mu g~ g^{-1}\partial_\mu g)+\frac{1}{4\pi}\int dx~\epsilon^{\mu\nu}\mathrm{tr}(g^{-1}\partial_{\mu}g~g^{-1}\partial_{\nu}g~g^{-1}\delta g).
\end{equation}
Let us assume that $G=SU(N)$ and the spacetime manifold is compactified as $S^2$. Then it is well-known that the second term on the right hand of the above equation is a variation of a Wess-Zumino term. Thus the free energy can be explicitly written as
\begin{equation}
W[g]=-\frac{1}{8\pi}\int d^2x~ \mathrm{tr}(g^{-1}\partial^\mu g~ g^{-1}\partial_\mu g)+\frac{1}{12\pi}\int d\tau d^2x~\epsilon^{\mu\nu\rho}\mathrm{tr}(\tilde{g}^{-1}\partial_{\mu}\tilde{g}~\tilde{g}^{-1}\partial_{\nu}\tilde{g}~\tilde{g}^{-1}\partial_\rho \tilde{g}),
\end{equation}
where $\tilde{g}(\tau,x)$ with $\tau\in[0,1]$ is a continuous extension of $g(x)$ with $\tilde{g}(0,x)=g(x)$ and $\tilde{g}(1,x)$ being constant.

Last, for the Dirac fields with both the spin and flavor as in Eq.(7) of the main text, the above derivation works but needs to be generalized with coupling to two non-Abelian gauge field $A_{\mu}$ and $B_{mu}$, resulting in the following Langrangian as,
\begin{equation}
\begin{split}
\mathcal{L} &=\bar{\psi}^{ia} (i\gamma^\mu \partial_\mu\delta^{ij}\delta^{ab}+A_\mu^{ij}\delta^{ab}+\delta^{ij}B_\mu^{ab})\psi^{jb}
   =\bar{\psi}(i\gamma^\mu \partial_\mu 1_n\otimes 1_m+A_\mu\otimes 1_m+1_n\otimes B_\mu)\psi,
\end{split}
\end{equation}
where accordingly we have$ A_+=ig_A^{-1}\partial_{+} g_A$, $A_-=ih_A\partial_{-} h_A$,
$B_+=ig_B^{-1}\partial_{+} g_B$,  and $B_-=ih_B\partial_{-} h_B$, such that
\begin{eqnarray}
A_+\otimes 1_m+1_n\otimes B_+=ig_A^{-1}\otimes g_B^{-1}\partial_+(g_A\otimes g_B)\\
A_-\otimes 1_m+1_n\otimes B_-=ih_A^{-1}\otimes h_B^{-1}\partial_+(h_A\otimes h_B)
\end{eqnarray}
For two arbitrary matrices $M$ and $N$, one has the following property $\mathrm{tr}(M\otimes N)=\mathrm{tr}(M)\mathrm{Tr}(N)$.
Moreover, for $g\in SU(n)$, $\mathrm{tr}(g\partial_\mu g^{-1})=0$ since the Lie algebra consists of $n\times n$ traceless Hermitian matrices. With the above two identities, it is straightforward to derive that the following WZW emerges:
\begin{equation}
W[g_A \otimes g_B]=MW[g_A]+NW[g_B].
\end{equation}
This is the decoupled WZW CFT, from which one can read of the fusion rules in order to obtain the Kondo fixed points, which are discussed in the main text for both cases, i.e., with and without the off-diagonal entries of the scattering $g_{ab}$ (Eq.(4) of the main text).

\section{perturbative RG calculation of $\beta$-functions of exchange couplings}
In order to determine the fixed points at the strong-coupling regime, we perform a perturbative RG calculation of the $\beta$-functions with respect to the derived effective 1D model, which reads as $H=H_0+H^{\prime}_p$, where $H_0$ is the rotated Hamiltonian in the valley space with respect to Eq.\eqref{eqs7}, which is of the form,
\begin{equation}\label{eqsn10}
H_0=\sum_{m=1,2}\sum_{\sigma}\int^{\infty}_{-\infty} d\epsilon(\epsilon-\mu)d^{\dagger}_{\epsilon,m,\sigma}d_{\epsilon,m,\sigma},
\end{equation}
and $H^{\prime}_p$ is approximated by using the density of states at the Fermi energy, $\rho_0$, leading to
\begin{equation}\label{eqsn11}
 H^{\prime}_p=g_1\int^{\infty}_{-\infty}d\epsilon d\epsilon^{\prime}d^{\dagger}_{\epsilon,1,\sigma}\boldsymbol{\tau}_{\sigma\sigma^{\prime}}\cdot\mathbf{S}_{imp}d_{\epsilon^{\prime},1,\sigma^{\prime}}
 +g_2\int^{\infty}_{-\infty}d\epsilon d\epsilon^{\prime}d^{\dagger}_{\epsilon,2,\sigma}\boldsymbol{\tau}_{\sigma\sigma^{\prime}}\cdot\mathbf{S}_{imp}d_{\epsilon^{\prime},2,\sigma^{\prime}},
\end{equation}
where a rotation in the valley space is performed, leading to the channel $m=1,2$, and $g_1=\rho_0(g_d+g_t)$, $g_2=\rho_0(g_d-g_t)$. The perturbative expansion over the two terms in $H^{\prime}_p$ can be constructed with Feynman diagrams to the two-loop order. Integrating out the fast mode momentum leads to the renormalization group flow as,
\begin{eqnarray}
  dg_1/dl &=& g^2_1-g_1(g^2_1+g^2_2)g_1/2, \\
  dg_2/dl&=&g^2_2-g_2(g^2_1+g^2_2)g_1/2.
\end{eqnarray}
where $dl=d\Lambda/\Lambda$ is the RG scaling parameter. The first term obtained from second order is relevant, showing the asymptotic free of the exchange coupling, and the second term from the third order contributes a suppression of the relevant flow, generating a channel-mixed fixed point with finite values of $g$'s, as shown by Fig.2(a) of the main text.

\section{The thermal conductivity from CS fermions}
We now list in this section details for calculation of the thermal conductivity, with respect to both of the two Kondo-generated fixed points. The starting point is Eq.\eqref{eqsn10} and $\eqref{eqsn11}$, which have been bosonized in the non-Abelian fashion before. The 1D model can be further mapped to half-infinite chain with left and right movers denoted by fields $d_{L/R,m,\sigma}$ \cite{Afflecka}. We consider the single-particle Green's function defined by $\langle d^{\dagger}_{L,m,\sigma}(z_1)d_{R,m,\sigma}(z_2)\rangle$, where $z$ lies in the complex plane representing for the 1+1D spacetime. We obtain $\langle d^{\dagger}_{L,m,\sigma}(z_1)d_{R,m,\sigma}(z_2)\rangle=0$ for the case with no boundary, and $\langle d^{\dagger}_{L,m,\sigma}(z_1)d_{R,m,\sigma}(z_2)\rangle_{Free}=1/(z_1-\overline{z}_2)$ for a trivial boundary (corresponding to the weak coupling regime). The correlation with respect to the Kondo fixed point can be calculated via the boundary state that in turn obtained by fusion \cite{Afflecka,Affleckb,Affleckc,Ludwig,Affleckd,Afflecke,Ludwigf}, leading to the scattering matrix $S$ connecting the correlations as $\langle d^{\dagger}_{L}(z_1)d_R(z_2)\rangle_{Kondo}=S\langle d^{\dagger}_{L,m,\sigma}(z_1)d_{R,m,\sigma}(z_2)\rangle_{Free}$ , where  $S=-1$ and $S=0$ for the FL and NFL fixed point.
Assuming a dilute impurity with density $n_{imp}$, the scattering time $\tau^{-1}_s=-2\mathrm{Im}\Sigma^R(\omega)$, where $\Sigma^R(\omega)$ is the retarded self-energy that is related to $S$ \cite{Afflecke}.

On the other hand, the thermal current can be readily derived from Eq.\eqref{eqsn10} and Eq.\eqref{eqsn11} as,
\begin{equation}\label{eqsn12}
\mathbf{j}_E=-v_F\hat{k}\sum_{m,\sigma}\int d\epsilon (\epsilon-\mu)d^{\dagger}_{\epsilon,m\sigma}d_{\epsilon,m,\sigma},
\end{equation}
where $\hat{k}=\mathbf{k}/k$ and $v_F$ is set to 1 in the following. The thermal conductivity can be evaluated through $T\sigma_E(\omega)=-\lim_{\omega\rightarrow0}\mathrm{Im}\pi^R(\omega)/\omega$, where $\pi^R$ is the retarded current-current correlation. We firstly calculate the current correlation function in Matsubara form as,
\begin{equation}\label{eqsn13}
\pi(i\omega_n)=-\frac{1}{3}\int^{\beta}_0 d\tau e^{i\omega_n\tau} \langle \hat{T}_{\tau}\mathbf{j}_E(\tau)\cdot\mathbf{j}_E(0)\rangle,
\end{equation}
where $\hat{T}$ is the imaginary time ordering. Inserting Eq.\eqref{eqsn12} into Eq.\eqref{eqsn13}, we obtain
\begin{equation}\label{eqsn14}
\pi(i\omega_n)=\frac{1}{3}\int d\epsilon\xi^2\frac{1}{\beta}\sum_{i\nu_n} g(\epsilon,i\omega_n+i\nu_n)g(\epsilon,i\nu_n),
\end{equation}
where $\xi=\epsilon-\mu$, and $g(\epsilon,i\omega_n)$ is the Matusbara Green's function of the d-fermions.  Performing the sum of Matusbara frequency $i\nu_n$ and after insertion of $\pi(i\omega_n)$ into $\sigma_E$, one obtains that
\begin{equation}\label{eqsn15}
T\sigma_E=\frac{1}{3}\int d\epsilon\xi^2\int d\nu\delta(\epsilon-\nu)\tau_s(-\frac{\partial}{\partial{\nu}}n_F(\nu)),
\end{equation}
where $n_F(\nu)=1/\mathrm{exp}[\beta(\nu-\mu)+1]$ is the Fermi distribution function.  We consider lowest order $T$ behavior of $T\sigma_E$. Since $\tau_s$ is dependent on $S$ which further relies on higher order $T$-terms, its $T$-dependence can be neglected for low $T$. After inserting the Fermi distribution function and completing the integrals, it is straightforward to find that the right-hand-side of Eq.\eqref{eqsn15} is proportional to $T^2$, leading to $\sigma_E(T)/T=\pi^2\tau_s/9=(\pi^3\rho_0)/[9(1-S)n_{imp}]$ for low $T$. With further taking into account the higher order $T$-dependence in $S$ \cite{Ludwigf,Afflecke}, one can obtain the different $T$-scalings of the thermal conductivity for the two fixed points, as shown by the main text.

\section{Details of numerical renormalization group calculations}
For the honeycomb lattice XY model discussed above, we derived an effective two-channel Kondo model with impurity spin $\mathbf{S}_{\text{imp}}=\frac{1}{2}$ and Kondo couplings $g_1$ and $g_2$.  This can be solved using numerical renormalization group. The density of states of the non-interacting bath takes a linear form,
\begin{equation}
\rho(\epsilon)=\frac{|\epsilon-\mu|}{2\pi},
\end{equation}
where, $\epsilon\in[-1, 1]$ and $\mu$ is the chemical potential. The bath is discretized logarithmically and mapped to a semi-infinite ``Wilson chain" with exponentially decaying hoppings, and the impurity coupled to the first chain site via Kondo constants $g_1$ and $g_2$. The chain is diagonalized iteratively while discarding high-energy states, thereby zooming in on low-energy properties: the finite-size level spacing of a chain ending at site $k$ is of order $\omega_k\propto \Lambda^{-k/2}$. Here $\Lambda>1$ is a discretization parameter, chosen to be 2 in this work. We use the full-density-matrix NRG~\cite{Andreasa} method to solve this model, exploiting its full U(1)$_{\text{charge}}$ $\times$ SU(2)$_{\text{spin}}$ $\times$ SU(2)$_{\text{channel}}$ symmetry when $g_1=g_2$ and U(1)$_{\text{charge}}$ $\times$ SU(2)$_{\text{spin}}$ when $g_1\ne g_2$ using the QSpace~\cite{Andreasb,Andreasc} tensor library.  We keep 4000 multiplets in the diagonalizations. The imaginary part of the impurity dynamical susceptibility, $\chi(\omega)=-\text{Im}\langle S||S\rangle_{\omega}$, was calculated at temperature $T=10^{-10}$. The FL [NFL] Kondo scale is determined as the energy at which $\chi(\omega)$ [$\frac{\text{d} (\log\chi(\omega))}{\text{d}(\log\omega)}$] has a maximum.


\section{Introduction}
The high-speed flight of an air-vehicle in supersonic and hypersonic flow regimes have the formation of shockwaves in front of the nose. The high pressure stagnation zone in the frontal shock wave can cause high wave drag, which leads to higher total drag of the vehicle. In high supersonic and hypersonic flow regimes, a major amount of flow kinetic energy will be converted to internal energy of the gas across the shock. It causes high aerodynamic heating in highly compressed stagnation zone\cite{erath}. Many active and passive drag control methods have been extensively studied by manipulation of the stagnation zone in supersonic and hypersonic flows such as mechanical spike, counterjet flow and energy deposition \cite{reviewp}. Among them, mechanical spike is a passive drag control method and counterjet flow and energy deposition are active drag control methods. Breathing blunt nose is another passive drag control methods studied in supersonic \cite{vashish} and hypersonic \cite{watanabe} flow regimes. The majority of the drag control methods work by manipulation of the frontal shock wave either by pushing away from the body as in the case of spike, counterjet flow and energy deposition or bringing closer the frontal shock and manipulating the base drag in breathing blunt nose. Among the passive flow controls, the main disadvantages of mechanical spike are unsteady flow-fields for certain length of spikes and it adds additional length to the vehicle which may produce additional side forces. The breathing blunt nose method was found to be promising but it may be required to further study the effectiveness in high enthalpy flows and heating control in internal flow. Among active control methods, the energy deposition method may require a complicated energy source device onboard either for electric discharge, laser or microwave based energy deposition. The counterjet flow have been found one of the effective technique to control the drag, which can be used to reduce the drag, when small mass flow is injected in counterjet flow manner or to increase drag (increase thrust, while re-entry as in retro-propulsion), when higher mass flow (mostly exhaust from engine) is ejected. However, most studies are performed with either air or non-reactive gases in drag reduction counterjet flow controls. Warren \cite{warren} studied helium and nitrogen injection and found that helium (lighter) gas was superior to reduce peak at pressure reattachment. A number of studies investigated hybrid drag control methods for example lateral jet injection from spike \cite{jiang} and electrical discharge from the spike tip \cite{kuo}. However, to author's knowledge, there are no hybrid studies involving energy addition (in the stagnation zone) and counterjet flow in order to achieve drag control. This study aims to investigate the feasibility of using hydrogen as counterjet flow gas in order to control the drag and compare the drag reduction with respect to counterjet air operating under the same conditions. The low ignition energy required to ignite hydrogen, its high diffusivity and being a lightest gas are the main reasons behind using hydrogen as a counterjet flow jet to control drag.
However, there may be concerns of high heating due to the high temperature flame in front of blunt nose. The current study is focused only on drag control by hydrogen injection in the stagnation zone. A complimentary heat transfer study will be carried out in future work.

To the best of author's knowledge, the current study is the first that targets drag control by direct hydrogen injection as counterjet flow in the stagnation zone. Srinivasan et al. \cite{srinivasan} have studied the hydrogen injection in a spiked tipped body and found 35-50 \% drag reduction. Zhao et al.\cite{zhao} have recently simulated non-equilibrium flow at hypersonic and shock induced combustion in hydrogen air mixture for spherical blunt nose as well as hydrogen injection (found flow behaviour close to air injection) from a numerical tool's capability point of view. Many studies have been performed by researchers for oblique detonation wave propagation in a hydrogen-air freestream premixed mixture \cite{qin} and traverse hot-jet induced detonation propagation\cite{cai} for a pressure gain propulsion system. It is well established that the counterjet flow may behave in short (SPM) and long penetration modes (LPM) \cite{hayashi}, which depends on stagnation pressure ratio of counterjet flow to freestream flow and nozzle design. Sonic (convergent) counterjet nozzles usually exhibit short penetration modes at most of operating pressures, while supersonic (convergent-divergent) counterjet nozzles exhibit long penetration mode (LPM) or short penetration modes (SPM), depends on operating stagnation pressure. In LPM mode operation, the shock will be moved at longer upstream distances, which can lead to higher drag reduction. In short penetration mode, the Mach disk can be formed near central region and a low pressure region can be formed near the blunt body, which can cause low drag on the body. The current study is based on axisymmeytric two-dimensional simulation of a counterjet flow injection from a hemispherical blunt nose using the open-source platform OpenFOAM \cite{weller} at hypersonic Mach number 6 and freestream pressure and temperature 10 kPa and 293 K, respectively. Additionally, high-pressure hydrogen injection near the center of hemispherical blunt nose is simulated using ddtFoam solver \cite{ettner}, which was developed and validated for deflagration to detonation transition in long tubes with blockages. To summarize, the objectives of the current numerical study are: 1) Investigate the effect of sonic hydrogen injection in stagnation zone of hemispherical blunt nose in hypersonic flow-field and compare with equivalent air injection. 2) Study the effects of supersonic injection of hydrogen in the stagnation zone of hemispherical blunt nose in hypersonic flow at different LPM or SPM modes.
\section{Numerical Method}
The numerical simulations are performed by solving unsteady compressible Navier-Strokes equations with $k-\omega$ SST turbulence closure model for an axisymmeytric domain. The ddtFoam \cite{ettner} solver at opensource OpenFOAM platform has been used for all the simulation of air injection as well as hydrogen ignition. The ddtFoam \cite{ettner} utilizes density based approximate Riemann solver with HLLC scheme for solving convective fluxes. The reaction mechanism of hydrogen air combustion is modelled as 9 species 19 reactions O’Conaire reaction mechanism \cite{conaire}. The deflagration and detonation combustion source terms in transport equation of reaction progress variable are evaluated using two look-up tables: reaction products table using mechanism in Cantera \cite{goodwin} and autoignition time delay using Cantera and EDL Toolbox \cite{brown}. The Weller model \cite{weller} has been used for the deflagration combustion modeling and detonation is modelled using auto ignition delay time. It is found that ddtFoam can predict various modes of DDT transition in obstacle laden macro shock-tube. As the ddtFoam solver is developed for modelling accidental fire in large domain of power plants, it utilizes a volume fraction method in each cell to determine the autoignition time-delay, which leads to adequate accurate solutions in coarser grid. On the other hand, to accurately capture the flame propagation, it is required to use 30-50 grid points in the flame thickness region, which may lead to higher computational cost. Hence, it can be said that ddtFoam solver can be utilized for various time dependent extreme combustion or detonation based problems and provide quick results with adequate accuracy for design and analysis purposes. This study further utilizes the ddtFoam solver to study the drag control by using high-pressure hydrogen injection in the stagnation zone. The details of computational domain, boundary conditions and solution procedure are outlines in the following subsections.
\subsection{Computational Domain:}
\begin{figure}[bht]
\centering
\includegraphics[width=0.8\textwidth]{Fig01.pdf}
\caption{\label{Fig01} (a) Computational Domain (b) Axisymmeytric Grid for preliminary study}
\end{figure}
The computational domain consists of a 50 mm diameter hemispherical blunt nose. The counterjet flow is injected through a 4 mm diameter hole at the center of blunt nose. Figure \ref{Fig01}a shows the overall computational domain. The computational domain extends 120 mm after the hemispherical blunt nose to establish the supersonic outlet. The incoming hypersonic flow have been modelled as air with an inlet boundary condition of Mach 6, freestream pressure of 10000 Pa and freestream temperature 293 K. Sonic and supersonic hydrogen or air injected from the center of the hole, which is modelled by imposing the inlet boundary conditions at a 2 mm arc near the center. The temperature of injection jet is maintained at 300 K for all the simulated cases and the stagnation pressure was varied in order to study their counterjet flow effects on drag control. The structured computational grids have been generated as shown in Fig. \ref{Fig01}b. To perform the axisymmetric simulations, the two-dimensional grid have been rotated for $5^{\circ}$ angle with symmetry in the X-Y plane. The wall is modelled as adiabatic wall and the outlet boundary is modelled as a supersonic outlet. The species at the outlet are modelled as a {\it zeroGradient} boundary condition. The numerical simulations for each case has been performed in two steps. Initially, a freestream hypersonic flow is established in front of hemispherical blunt nose by completing the simulation for 1 ms. The hole at the center is considered as adiabatic wall. In the second step, air or hydrogen is injected at sonic or supersonic speed with a temperature of 300 K and varied stagnation pressure. The second simulation is performed for an additional 2 ms physical time, which was sufficient to establish the flow-field in countejet flow.
\subsection{Grid Independence Study:}
\begin{figure}[htb]
\centering
\includegraphics[width=0.5\textwidth]{Fig02.pdf}
\caption{\label{Fig02} Air Injection at $PR$= 0.8, (a) Schlieren Visualization \cite{hayashi} (b) Mach contours from ddtFoam for medium and fine grid}
\end{figure}
In order to determine the effectiveness of the ddtFoam solver in counterjet flow application, one of the experimentally studied conditions by Hayashi et al.\cite{hayashi} was simulated for sonic counterjet flow with similar operating condition. The results and operating conditions are shown in Fig.\ref{Fig02}. As the ddtFoam solver was developed for atmospheric or high pressure detonation conditions, the very low pressure and temperature conditions may not be accurately modelled. Hence, in the simulated cases the freestream temperature was changed to 293 K and the sonic injection of air is simulated in place of the $N_2$ injection in experimental case. While the simulated case and experimental case slightly differ in their operating conditions, their flow features can be qualitatively compared. Also, a grid independence study has been performed for the same operating case for three grid sizes and the results are shown in Table 1. Figure\ref{Fig02}a and \ref{Fig02}b shows the qualitative comparison and captured flow-field related to sonic injection for the medium and fine grid configurations. The main flow features such as modified bow shock, Mach disk and the recirculation zones are well captured using the ddtFoam solver for air injection. As this study is mainly focused on drag control, the calculated drag coefficient without flow injection and with sonic injection for three different grid sizes are tabulated in Table 1. The mass flow rate has also been tabulated for all three cases. The drag coefficient for the no injection case shows less than 1\% deviation and with air injection shows less than 2\% deviation when comparing medium and fine grid cases. Hence, medium grid of $141 \times 101$ has been used for all the simulated cases.
\begin{center}
\begin{table}[!htb]
\centering
\caption{\label{jfonts}Grid Independence Study Results.}
\begin{tabular}{@{}l*{5}{l}}
\br
Grid Size & Drag Coefficient & Drag Coefficient & Mass Flow Rate\\
        & (No Air Injection) & Air Injection (PR = 0.8) & (gms /s)\\
\mr
$71 \times 51$ &0.9736 & 0.6753 & 29.17\\
$141 \times 101$ &0.9788 & 0.6981 & 29.17\\
$201 \times 151$ &0.9797 & 0.6869 & 29.17\\
\br
\end{tabular}
\end{table}
\end{center}
\section{Results and Discussions}
The effectiveness of counterjet flow control depends on several parameters, which have been modelled in many previous studies \cite{reviewp}. The three main parameters for governing performance of counterjet flow drag controls are injection mass flow rate ($\dot{m_j}$), stagnation pressure ratio of injecting jet and freestream flow ($PR$) and momentum ratio ($R_{MA}$) for injecting flow and freestream flow. These parameters can be defined as follows in terms of free stream or injection density, velocity and cross-section area:
\begin{equation}
\dot{m_j} = \rho_{j}V_{j}A_{j}
\end{equation}
where, $\rho_j$, $V_j$ and $A_j$ represent the counterjet exit density, velocity and cross section area, respectively.
\begin{equation}
 PR = \frac{P_{0,j}}{P_{0,\infty}}
\end{equation}
where, $P_{0,j}$ and $P_{0,\infty}$ represent the counterjet stagnation pressure and freestream flow stagnation pressure, respectively.
\begin{equation}
 R_{MA} = \frac{\rho_{j}V_{j}^2 A_{j}}{\rho_{\infty}V_{\infty}^2A_{\infty}}
\end{equation}
additionally, here $\rho_{\infty}$, $V_{\infty}$ and $A_{\infty}$ represent the freestream density, freestream velocity and cross section area of blunt nose, respectively. \\
The total drag in each simulations have been calculated by using the following equation:
\begin{equation}
 C_{D,Total} = C_{D,A} + \dot{m_j}V_j + (P_j-P_{\infty})A_j
\end{equation}
here, $C_{D,Total}$ is total drag force experienced by blunt nose, $C_{D,A}$ is drag force experienced by blunt nose (pressure + viscous), and the next two terms represent drag due to thrust generated by counterjet flow.
In this study, the freestream conditions for air have been kept fixed as hypersonic Mach number $M_{\infty} = 6$, freestream pressure $P_{\infty} = 10000$ Pa and temperature $T_{\infty} = 293$ K. The cross-sectional area of hemispherical blunt nose $A_{\infty}$ and cross-sectional area of counterjet exit $A_{j}$ are also fixed. The counterjet exit temperature, $T_{j} = 300$ K has also been kept fixed. It should be noted that the drag reduction can be achieved by counterjet flow, when the momentum ratio is significantly less than 1 in order to generate the least thrust from the counterjet. If the momentum ratio is higher than 1, it will produce high thrust and work as retro-propulsion. The stagnation pressure ratio $PR$ signifies the strength of a counterjet flow, it also governs the mass flow rate of a counterjet. In this study, when the comparison of the effectiveness of a counterjet flow between air and hydrogen was made, $PR$ and $R_{MA}$ were kept constant, which lead to a smaller mass flow rate for hydrogen injection than air and higher exit velocity for hydrogen than air. The results are discussed for sonic ($M_j = 1$) injection of hydrogen and air, and supersonic injection ($M_j = 2.94$)of hydrogen and air. The sonic counterjet injection cases for air and hydrogen injections have been compared for a single $PR = 0.12$, which shows the flow-features as short penetration mode. Further more supersonic counterjet scenarios have been compared for air and hydrogen injection for different $PR$ values in order to simulate short and long penetration modes.

\subsection{Sonic Air and hydrogen Injection:}
Figure 3 shows Mach contours for a blunt nose without any counterjet injection in the freestream air with a Mach number 6. The computed drag coefficient for the hemispherical blunt nose is 0.94. Figure \ref{Fig03}b and \ref{Fig03}c show the flow features for air and hydrogen counterjet flows exiting from the sonic nozzle at $PR$= 0.12. For both cases, the momentum ratio remains constant as 0.0178. For this stagnation pressure ratio, the air and hydrogen counterjet flows show similar flow features in their Mach contours. However, the bow shock at the stagnation line moves slightly away from body for hydrogen injection in comparison to air injection. However, there is no significant effect on drag reduction, as both the gas injection scenarios lead to approximately a 21 \% drag reduction. However, the required mass flow rate for hydrogen is almost 3.76 times less than the required mass flow rate of air to achieve the same performance.
\begin{figure}[htb]
\centering
\includegraphics[width=0.8\textwidth]{Fig03.pdf}
\caption{\label{Fig03} Mach contours for (a) Blunt Nose (b) Sonic Air Injection (c) Sonic Hydrogen Injection}
\end{figure}
\begin{figure}[htb]
\centering
\includegraphics[width=0.9\textwidth]{Fig04.pdf}
\caption{\label{Fig04} (a) Temperature, (b) Hydrogen Mixture Fraction (c)$H_{2}O$ Mass Fraction, (d) OH Mass Fraction for sonic hydrogen injection case}
\end{figure}
Fig. \ref{Fig04} shows the contour plots for temperature, hydrogen mixture fraction, $H_{2}O$ mass fraction and $OH$ mass fraction for hydrogen injection by the sonic nozzle at $PR$= 0.12. The mixture fraction was modelled as unburned fuel mass fraction if the reaction is not present \cite{ettner}. The temperature plot shows a high temperature zone away from the blunt body, which can be clearly seen as a standing flame in the $OH$ mass fraction plot. The $H_{2}O$ mass fraction shows that the reaction product can completely submerge the blunt body. The main finding is that the counterjet of a lighter or heavier gas may provide similar levels of drag control in SPM modes for the same momentum ratio and PR, but the lighter gas requires a lower mass flow rate resulting in reduced storage requirements.
\subsection{Supersonic Air and hydrogen Injection:}
\begin{figure}[htb]
\centering
\includegraphics[width=0.9\textwidth]{Fig05.pdf}
\caption{\label{Fig05} (a) Mass Flow Rates (b) Thrust and Drag Coefficient variations with change in $PR$ for supersonic ($M_j = 2.94$)air or hydrogen injection}
\end{figure}
\vspace{-1cm}
\begin{center}
\begin{table}[h]
\centering
\caption{\label{table2} Observed behaviour of Supersonic Counterjet flow.}
\begin{tabular}{@{}l*{6}{l}}
\br
PR & 0.20 & 0.30 & 0.40 & 0.50 & 1.0 & 2.0\\
\mr
Air & LPM & LPM & LPM & LPM & SPM & SPM \\
Hydrogen & Unstable & Oscillatory & Oscillatory & Oscillatory & Intermittent & Intermittent \\
         & Solution & LPM & LPM & LPM & SPM-LPM & SPM-LPM \\
\br
\end{tabular}
\end{table}
\end{center}
The sonic counterjets generally behave in a short penetration mode. In order to achieve a long penetration mode with a wide range of operating pressures, the use supersonic nozzle is required. The supersonic injection of air or hydrogen (at $M_j = 2.94$) have been modelled as inlet conditions at the center hole location of a blunt nose with same diameter as 4 mm. The simulations were performed by changing the stagnation pressure ratio from $PR$= 0.25 to 2.0 as shown in Fig. \ref{Fig05}. The momentum ratio $R_{MA}$ for each $PR$ have been constant for air and hydrogen injection, which results in a 3.76 times lower mass flow rate for hydrogen injection in comparison to air injection (Fig. \ref{Fig05}a). The hydrogen injection for $PR$ less than 0.3 were not stable during the simulation, hence were not included in the results. The time-averaged aerodynamic drag (including pressure and viscous drag on the blunt nose), drag due to thrust by counterjet flow and total drag in each case of air and hydrogen injection are plotted in Fig. \ref{Fig05}b. Table \ref{table2} summarises the behaviour of counterjet for both air and hydrogen injection. It can be seen that the drag due to thrust from the counterjet increases linearly as $PR$ increases for both hydrogen and air, because the $R_{MA}$ remains constant for each $PR$. In the case of air injection, the aerodynamic drag decreases with an increase in PR and the counterjet changes its behaviour from LPM to SPM. However, at higher $PR$, the thrust also increases, which leads to an over-all increase in the total drag at higher $PR$. For hydrogen injection, $PR$ = 0.3, 0.4 and 0.5 exhibit a long penetration mode and a drag reduction up to 60 \%. The flow-field is oscillatory in nature. At higher $PR$ = 1 and 2 the aerodynamic drag increases as flow-field changes to SPM for a longer time and intermittently change to LPM mode because of the high pressure burning. It can be concluded that at low $PR$, operating at a smaller mass flow rate, hydrogen injection provides higher drag reduction by pushing the shock wave far away from blunt body. A future heat transfer study will provide a clear idea of the overall effectiveness of hydrogen injection in LPM mode for drag as well as aerodynamic heating control.

\section{Conclusions}
The effectiveness of a hybrid drag control method for a sonic and supersonic hydrogen injection in the stagnation zone of a blunt nose at hypersonic Mach number 6 have been numerically investigated. The sonic hydrogen injection lead to a SPM flow-field, similar to air injection in front of the blunt nose and provide same drag reduction performance at 26 \% lower mass flow rate. The supersonic hydrogen injection behaves differently to the supersonic air injection, leading to a higher drag reduction up to 60 \% at lower $PR$ (= 0.3 to 0.5) and mass-flow rate as low as 5.6 - 10 g/s. At higher $PR$ values, air injection performs better than hydrogen injection, but may require very high mass flow rates up to 140 g/s. Future work will focus on a heat transfer study to analyse the effects of hydrogen injection on drag as well as aerodynamic heating at LPM modes.

\section*{References}
\medskip


\section{Motivation and significance}
\label{sec:intro}

Gravitational waves were originally predicted by Einstein in
1916~\cite{Einstein:1916cc} as a consequence of general relativity, which
describes gravity as the warping of space and time caused by mass and
energy~\cite{Einstein:1915ca}.  Two extremely massive objects orbiting one
another e.g., black holes or neutron stars, warp space dynamically and send
ripples across the universe that can be observed here on Earth. As they pass
by, gravitational waves stretch and squeeze the space around Earth by less than the
width of an atom compared to the Earth's diameter.  Due to their tiny effect on
scientific instruments, gravitational waves were not observed until 100 years
after their initial prediction. Technological advances in laser interferometry
led to the discovery of gravitational waves from a merging binary black hole in
2015~\cite{Abbott:2016blz}. This watershed moment was made possible by the
Advanced Laser Interferometer Gravitational-wave
Observatory (LIGO) ~\cite{TheLIGOScientific:2014jea} and the scientists of the LIGO and
Virgo Collaborations.

\begin{figure}[h!]
\includegraphics[width=\textwidth]{infographic.pdf}
\caption{\label{fig:info} Gravitational wave infographic.  Gravitational wave
data is time series, audio frequency data that is noise dominated.  GstLAL
identifies signals consistent with the predictions of general relativity as
measured by multiple gravitational wave detectors and assesses the probability
that these signals come from merging neutron stars and/or black holes in near
real-time.}
\end{figure}

Advanced LIGO and Advanced Virgo~\cite{TheLIGOScientific:2014jea,
TheVirgo:2014hva} are the currently operating worldwide network of
kilometer-scale laser interferometric gravitational wave observatories which
have measured gravitational wave signals. These detectors provide a new way to
observe our Universe and enable a vast amount of new science. LIGO's
observations have already deepened our understanding of the populations of
compact objects such as neutron stars and black
holes~\cite{LIGOScientific:2018mvr, LIGOScientific:2018jsj}, and they have also
offered new tests of fundamental physics~\cite{Abbott:2017xzu, Abbott:2018lct,
LIGOScientific:2019fpa, Abbott:2019yzh}.  The strong gravity regime probed by
compact binary mergers is a laboratory for novel tests of general relativity,
and a joint observation of gravitational waves along with electromagnetic
waves~\cite{GBM:2017lvd} has taught us how matter behaves in the most extreme
conditions~\cite{Abbott:2018exr, Nicholl:2017ahq}.

The science made possible by LIGO and Virgo is reliant on measuring
miniscule changes in the arm lengths of the interferometers known
as strain. The perturbations caused by incident gravitational waves manifest
themselves as variations in the intensity of laser light output.
\emph{Detector calibration} aims to accurately map the intensity of the output
to differential changes in the arm length through real-time signal-processing.
The calibrated strain data contains the encoded properties of the astrophysical
systems that produce gravitational waves.  The analysis of this data is
complicated by the presence of a vast array of transient noise sources.
\emph{Detector characterization} aims to quantify departures from stationary
noise to identify times where instrumental issues are so severe that the data
should not be analyzed or there may be a coupling between environmental sensors
(such as seismometers) and the gravitational wave strain data. Once data is
calibrated and assessed for quality, it is analyzed by a host of detection
algorithms to identify potential gravitational wave signals. In many cases the
signals are invisible to the naked-eye in raw data and discovering them
requires sophisticated techniques that may involve checking millions of models
against each segment of data.  All three of these activities require
substantial cyberinfrastructure.  The GstLAL software framework~\cite{gstlal}
was initially designed to support low-latency compact binary searches to
facilitate multi-messenger astronomy, but since its conception it has grown to
be a key component of the software used to produce accurately calibrated strain
data~\cite{Viets:2017yvy}, and recently it has contributed to detector
characterization
efforts~\cite{idq-methods, idq-in-gstlal}. The GstLAL framework is now
contributing key cyberinfrastructure to all three of these key aspects of
gravitational wave data analysis.  This paper will describe how the GstLAL
software is used in gravitational-wave
searches~\cite{Messick:2016aqy,Sachdev:2019vvd}, provide examples, and describe
the history and impact of GstLAL on gravitational wave discovery.

\section{Software description}
\label{sec:description}

Gravitational wave strain data quantifies how the distance between two points
will change as a gravitational wave passes. The current gravitational wave
observatories are sensitive to changes in strain and measure the stretching and
squeezing of space as a function of time.  The Advanced
LIGO~\cite{TheLIGOScientific:2014jea} and Advanced
Virgo~\cite{TheVirgo:2014hva} gravitational wave detectors are most sensitive
to strain frequencies between 10Hz--10kHz, which is remarkably close to the
frequency range of the human ear~\cite{olson1967music}. For this reason, there
is a close connection between the analysis of gravitational wave data and the
analysis of audio data.  Indeed, techniques such as low pass filtering, high
pass filtering, channel mixing, and gating apply equally well to both audio
processing and gravitational wave data processing, which provides the
motivation for basing GstLAL on Gstreamer~\cite{gstreamer}.

Gstreamer~\cite{gstreamer} is an open-source, cross-platform multimedia
processing framework designed to execute audio and video processing graphs
organized into three basic elements: sources, filters and sinks, which are
provided by dynamically loaded plugins.  A valid Gstreamer graph,
called a pipeline, connects elements together ensuring that the capabilities of
each element are satisfied.  The data are passed along in \emph{buffers} that
store both the memory location of the raw data as well as rich metadata.
Pipelines can be used to construct complex workflows and scale to thousands of
elements. GstLAL combines standard Gstreamer signal processing elements with
custom elements to analyze LIGO strain data.
\begin{figure}
\includegraphics[width=\textwidth]{gst.pdf}
\caption{\label{fig:gst}A basic Gstreamer graph.  Data starts at a source,
``src", e.g., a file on disk or a network socket, and then is passed through a
filter element, ``filter 1", which transforms the data, e.g., by performing a
low pass filter.  A second data stream starts from ``control src" and the
output of filter 1 is moderated by a gate controlled by the state of ``control
src".  The output of the gate is filtered through ``filter 2" and sent to a
sink which could be another file on disk or a network socket.}
\end{figure}

The GstLAL software began development in 2008 through the exploration of novel
techniques for filtering gravitational wave data~\cite{Cannon:2010qh}.  It
derives its name from ``Gstreamer wrappings of the LIGO Algorithm Library\footnote{\url{https://git.ligo.org/lscsoft/lalsuite}}
(LAL)". GstLAL began to take on its modern form by 2009 and has been since
actively developed as open source software. GstLAL currently resides in the
LIGO Scientific Collaboration hosted GitLab instance at
\url{https://git.ligo.org/lscsoft/gstlal}~\cite{gstlal}.

GstLAL is primarily a mix of Python and C with contributions from $75$ authors
distributed across North America, Europe, Asia and Australia.  The master
branch currently has over 13,000 commits and 250,000 lines of code.  GstLAL is
released under the GPLv2 license with 44 distinct releases since
2011~\cite{ligo-software}.  In 2012, code solely used for gravitational wave
searches for compact binaries was split off into its own package:
GstLAL-Inspiral. GstLAL-Inspiral has had 45 distinct releases since then.  In
2014, code used for gravitational wave burst detection was split into its own
package, GstLAL-Burst, with nine distinct releases. And finally, in 2014 code
used primarily for LIGO strain data calibration was split into its own package,
GstLAL-Calibration, with 60 releases. In addition to tar-ball releases, RedHat
and Debian compatible packages were produced for the LIGO Data Grid reference
platforms~\cite{ligo-packages}.

At present, Docker containers with the full GstLAL/LALSuite software stack are
built and distributed using the LIGO Container Registry~\cite{ligo-containers}.
The containers are built on top of the CentOS-based Scientific Linux 7, which
currently serves as the reference operating system on the LIGO Data Grid.
Binary executables are linked against Intel's high-performance Math Kernel
Library (MKL), and compiled to leverage Advanced Vector Extensions. Optimized
versions of the GstLAL software stack tuned at the compiler-level to best
leverage the native features of the local computing environment are often
custom-built by users. Previous studies have demonstrated a $\gtrsim 2$ times
increase in overall code throughput as a result of software tuning.  GstLAL is
also available through CondaForge~\cite{conda-forge}.

\subsection{Software Architecture}
\label{ssec:architecture}

Prior to 2015, gravitational waves had not been directly
observed~\cite{Abbott:2016blz}.  Although analysis techniques had been studied
for decades~\cite{Finn:1992xs, Allen:2005fk}, the character of gravitational
wave data evolved with the interferometers detectors during the initial
operation of LIGO and Virgo from 2002--2010~\cite{Abbott:2003pj, Abbott:2005pe,
Abbott:2007xi, Abbott:2009tt, Abbott:2009qj, Colaboration:2011np}.  Therefore,
we aimed to make the GstLAL software modular and easy to adapt to the
challenges of Advanced LIGO and Virgo data.

The joint observation of gravitational waves and electromagnetic signals, known
as multi-messenger astronomy, was a significant goal for Advanced LIGO and
Advanced Virgo.  In this scenario, gravitational-wave observations were
expected to be followed up by observations with telescopes across the
electromagnetic spectrum hoping to catch a short-lived transient light source.
Discovering gravitational waves quickly is critical because electromagnetic
counterparts may quickly fade.  GstLAL was designed to offer analysts an
extremely short time-to-solution to help ensure that electromagnetic
counterparts could be observed quickly.

The key design principles of GstLAL are:
\begin{description}
\item{\textbf{Plugin-based:}} Libraries within GstLAL provide Gstreamer
plugins to perform gravitational wave specific signal processing tasks. These
are mixed together with stock Gstreamer plugins to produce gravitational wave
analysis pipelines.  Plugins provide elements, which are the building blocks of
signal processing workflows.  These elements can be ordered in multiple ways
with minimal coding effort which allows for quick exploratory work and
development of new methods.
\item{\textbf{Streaming:}}  The goal of Gstreamer is to provide ultra low
latency signal processing suitable for audio and video playback and editing.
GstLAL pipelines typically work with streaming data buffers $\lesssim 1$ s in
duration.
\item{\textbf{Event-driven:}} GstLAL analysis pipelines are designed to run
continuously as data is collected.  Each application runs an event loop which
controls both application level operations as well as settings within a given
plugin.  This allows for the control of program behavior to be altered while
the application is running.  Dynamic program control is facilitated through
embedding the microservices framework Bottle~\cite{bottle}.  Simple http APIs
push information to the program, retrieve information or alter the program's
behavior.
\item{\textbf{Scalable:}} The GstLAL framework is designed both to scale to
dozens of cores on a single computer using the multi-threading provided by
Gstreamer and to scale to thousands of cores across a compute cluster
leveraging HTCondor directed acyclic graph (DAG) scheduling~\cite{condor}.
\end{description}

The GstLAL project is currently comprised of five distinct packages. 1)
\texttt{gstlal}, 2) \texttt{gstlal-ugly}, 3) \texttt{gstlal-inspiral}, 4)
\texttt{gstlal-calibration} and 5) \texttt{gstlal-burst}, all of which are
described below.

\subsection{{\upshape \texttt{gstlal}} package}
\label{ssec:gstlal}

The \texttt{gstlal} package curates a set of core plugins, functions and
applications that are used by nearly every analysis workflow developed within
GstLAL. The \texttt{gstlal} package is a dependency for all of the remaining
packages which are described in subsequent sections.  The \texttt{gstlal}
package provides Gstreamer elements for Finite-Impulse-Response filtering
(\texttt{lal\_firbank}), $N\to M$-channel matrix operations
(\texttt{lal\_matrixmixer}), data whitening (\texttt{lal\_whiten}), and data
gating (\texttt{lal\_gate}).  The \texttt{gstlal} package also provides basic
python APIs for building Gstreamer pipelines in the module \texttt{pipeparts},
basic data access routines in the module \texttt{datasource}, and a base class
for event handling in the module \texttt{simplehandler}.

\subsection{{\upshape \texttt{gstlal-ugly}} package}
\label{ssec:gstlal-ugly}

The \texttt{gstlal-ugly} package is an incubator package for software that is
in development. Eventually all \texttt{gstlal-ugly} software is migrated to one
of the other packages.

\subsection{{\upshape \texttt{gstlal-inspiral}} package}
\label{ssec:gstlal_inspiral}

The primary purpose of the \texttt{gstlal-inspiral} package is to house the
GstLAL-based search for compact binaries~\cite{Messick:2016aqy,
Sachdev:2019vvd}, which centers around the application
\texttt{gstlal\_inspiral}.  The GstLAL-based compact binary pipeline was
created to make near real-time gravitational-wave detections and aimed to one
day detect electromagnetically bright systems before
coalescence~\cite{Cannon:2011vi,Sachdev:2020lfd}.

The GstLAL-based compact binary search is a matched-filter search that
incorporates efficient time-domain filtering~\cite{Cannon:2010qh,
Cannon:2011tb, Cannon:2011xk, Cannon:2011rj, Cannon:2012gq, Smith:2012du,
Cannon:2011vi, Tsukada:2017cuf} of a set of template waveforms that match the
gravitational wave signals of merging black holes and neutron
stars~\cite{Mukherjee:2018yra}.  LIGO and Virgo detectors are prone to bursts
of nonstationary noise called glitches~\cite{TheLIGOScientific:2016zmo} and
determining the difference between gravitational waves and glitches is well
suited for many classification algorithms. GstLAL Inspiral implements a
classification scheme that is a hybrid of hypothesis testing techniques with
some elements coming from machine learning approaches.  The classifier is an
approximate likelihood ratio comprised of many terms which began as a custom
implementation of Naive Bayes classification~\cite{zhang2004optimality} applied
to gravitational wave searches~\cite{Cannon:2008zz}.  It was realized early in
the project that two things were apparent.  First, it wasn't practical to treat
the classifier as entirely data driven relying purely on training sets.
Training sets to adequately classify the full parameter space were too
expensive to produce.  Second, correlations between some parameters had to be
tracked in order to classify well.  The first point was addressed by developing
semi-analytic models to describe parts of the data~\cite{Hanna:2019ezx} and the
second point was addressed by factoring the multi-dimensional likelihood ratios
into groups of lower dimensional, but not one dimensional,
distributions~\cite{Cannon:2015gha, Cannon:2012zt}.

The GstLAL-based compact binary search has two modes.  The first is a
near-real-time, ``low-latency" mode that discovers and reports compact binaries
within tens of seconds of the signal arriving at Earth.  The second is an
``offline" mode that efficiently processes data in batch jobs where
time-to-solution is not as important as computational efficiency and
reproducibility.  Although both modes share $\gtrsim 95\%$ of the same code,
their behavior and design are very different in order to address the differing
concerns of real-time vs. batch processing.

The low-latency mode is a collection of typically $\mathcal{O}(1000)$
microservices that communicate (modestly) with one another asynchronously
through http using python Bottle, through an Apache Kafka queue, and through a
shared file system.  Each one of these microservices processes a portion of the
nearly 2 million models used in the current low-latency compact binary search.
The low-latency workflow is designed to be fault tolerant. If a job dies,
another is restarted to take its place.  Since information is exchanged
asynchronously and there is no guarantee of job success, the behavior in this
mode is non-deterministic.  In contrast, the offline mode has a fully
deterministic execution which can be reproduced to floating point precision.
The determinism is imposed by organizing each job in a directed acyclic graph
(DAG) using HTCondor.

\subsection{{\upshape \texttt{gstlal-burst}} package}
\label{ssec:gstlal_burst}

The GstLAL-based burst package is a collection of utilities intended to
search for gravitational-wave sources other than compact binaries as well
as non-astrophysical noise transients. One of the
recent developments is a pipeline searching for cosmic strings, which are
hypothetical objects considered to have formed in the early universe. The
pipeline uses time-domain stream-based signal processing algorithms, along
with a classification scheme using parameters specific to the search.
The algorithms used are mostly in common with the \texttt{gstlal-inspiral}
package, but simplified due to the smaller number of templates required for
matched filtering.

In addition, \texttt{gstlal-burst} provides utilities to identify and extract
features from non-Gaussian noise transients in near real-time ($O(5 s)$) via the
Stream-based Noise Acquisition and eXtraction, or SNAX toolkit. The SNAX toolkit
also leverages time-domain signal processing, but instead utilizes a sine-Gaussian
basis to identify the presence of and extract features from many types of
non-Gaussian noise in strain and auxiliary data. Its main data product is
multivariate time-series data containing the extracted features, including
SNR and phase information as well as the waveform parameters of interest.

\subsection{{\upshape \texttt{gstlal-calibration}} package}
\label{ssec:gstlal_calibration}

The \texttt{gstlal-calibration} package houses the unique software used for calibration of the LIGO strain data.  Software in the \texttt{gstlal-calibration} package produces the only official LIGO strain data product used in all subsequent analysis.  Calibration of LIGO strain data involves standard signal processing and digital filtering techniques in order to derive the differential arm motion observed in the LIGO detectors from the detector's digital readouts \cite{Abbott:2016jsd, Viets:2017yvy, Tuyenbayev2016}.  Many of the signal processing and digital filtering plugins used by the calibration pipeline \texttt{gstlal\_compute\_strain} are housed in the \texttt{gstlal} or \texttt{gstlal-ugly} packages.  A few plugins unique to the calibration process as well as calibration-specific python APIs are housed in the \texttt{gstlal-calibration} package.

Much like the GstLAL-based compact binary pipeline, the LIGO calibration pipeline is built to operate in two modes: a ``low-latency" mode and an ``offline" mode.  The low-latency LIGO calibration pipeline operates on hardware located physically at the two LIGO detector sites, LIGO Hanford  in Hanford, WA and LIGO Livingston in Livingston, LA.  This pipeline produces calibrated LIGO data and a bit-wise state-vector that indicates the fidelity of the calibrated data within $O(5 s)$ for each detector at the respective detector sites.  The low-latency calibration process involves using a combination of digital filtering performed in the LIGO front-end computers, which are directly connected to the LIGO detectors and employ the CDS Real-time Code Generator (RCG) core software \cite{RCGReleaseNotes}, and further digital filtering and processing performed by the GstLAL calibration software running on non-front-end hardware located at the LIGO detector sites.  This two-step process takes advantage of the access the front-end computing system has to the installed detector filters and models and the advanced stream-based filtering techniques housed in the GstLAL software packages \cite{Viets:2017yvy}.

There is often a need to re-calibrate the strain data after the initial low-latency data calibration in order to improve calibration accuracy based on more sophisticated modeling or to remove systematic errors present in the low-latency calibrated strain data \cite{Viets:2017yvy, Cahillane2017, Sun2019}.  The re-calibration of LIGO data is performed using the offline mode of the \texttt{gstlal\_compute\_strain} pipeline.  In this mode, the entire calibration process is performed by software housed in the GstLAL software packages.  The offline calibrated data is processed in batch jobs using HTCondor in order to optimize computational efficiency and is completely reproducible to floating-point precision.  All analyses derived from LIGO strain data use the calibrated data produced either by the low-latency GstLAL calibration pipeline or the offline GstLAL calibration pipeline.

\section{Illustrative Examples}
\label{sec:examples}

\subsection{Example Gstreamer pipeline with GstLAL}

The following example was run on a newly instantiated CentOS 7 64 bit virtual
machine with miniconda~\cite{miniconda} installed by doing:
\begin{lstlisting}[language=Bash]

wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
conda create -n myenv python=2.7
conda activate myenv
conda install -c conda-forge gstlal-inspiral=1.7.3

\end{lstlisting}
To verify that it works, try:
\begin{lstlisting}[language=Bash]

$ gst-inspect-1.0 gstlalinspiral | head -n4
Plugin Details:
  Name                     gstlalinspiral
  Description              Various bits of the LIGO Algorithm Library wrapped in gstreamer elements
  Filename                 <your miniconda path>/lib/gstreamer-1.0/libgstgstlalinspiral.so

\end{lstlisting}
Next we will try a very simple Gstreamer pipeline that uses two GstLAL
elements: \texttt{lal\_peak} and \texttt{lal\_nxydump}, along with additional
Gstreamer elements to construct a pipeline that generates 10 Hz Gaussian, white
noise, finds the peak sample every second and streams the result to the
terminal screen as ASCII text.  It is possible to construct simple pipelines
such as this without any code using the Gstreamer tool
\texttt{gst-launch}~\cite{gst-launch}:
\begin{lstlisting}[language=Bash]

$ gst-launch-1.0 audiotestsrc wave=9 ! capsfilter caps=audio/x-raw,rate=10 ! lal_peak n=10 ! lal_nxydump ! filesink location=/dev/stdout

\end{lstlisting}
The first element, \texttt{audiotestsrc} is a Gstreamer element that can
provide many test signals.  The \texttt{wave=9} property sets it to be unit
variance white noise. The second element, \texttt{capsfilter} specifies that we
want the format of the output to be floating point audio data with a sample
rate of 10 Hz. Next, \texttt{lal\_peak} is the first GstLAL element. In this
example it is configured to find the largest absolute value of the signal every
10 sample points.  \texttt{lal\_nxydump} is the second GstLAL element which
converts the time-series data to two column ASCII text.  Finally,
\texttt{filesink} dumps the ASCII output to standard out.  You should see the
following output (with variations caused by the fact that the data is random):
\begin{lstlisting}[language=Bash]{GstLaunchOut.sh}

$ gst-launch-1.0 audiotestsrc wave=9 ! capsfilter caps=audio/x-raw,rate=10 ! lal_peak n=10 ! lal_nxydump ! filesink location=/dev/stdout | head -n 15
Setting pipeline to PAUSED ...
Pipeline is PREROLLING ...
Pipeline is PREROLLED ...
Setting pipeline to PLAYING ...
New clock: GstSystemClock
0.000000000     0
0.100000000     0
0.200000000     0
0.300000000     0
0.400000000     0
0.500000000     0
0.600000000     0.95523328
0.700000000     0
0.800000000     0
0.900000000     0

\end{lstlisting}
you can see that the maximum sample point was chosen in the 10 sample interval
on line 13. Other values are set to 0.

\texttt{gst-launch} is useful tool for quickly testing a simple pipeline, or
debugging, however it is not suitable for writing large applications with many
elements or situations where program control is exposed dynamically to the
user.  For building applications, GstLAL relies on the Python bindings for
Gstreamer and adds a substantial amount of gravitational wave specific
application code written in python.  An example
of the pipeline above written in the style of GstLAL applications is below.
\begin{lstlisting}[language=Python]{demo.py}

# boiler plate Gstreamer imports
import gi
gi.require_version('Gst', '1.0')
from gi.repository import GObject, Gst
GObject.threads_init()
Gst.init(None)

# Gstlal imports
from gstlal import datasource
from gstlal import pipeparts
from gstlal import simplehandler

# initialize an event loop, a pipeline and an event handler
mainloop = GObject.MainLoop()
pipeline = Gst.Pipeline("softwarex_demo")
handler = simplehandler.Handler(mainloop, pipeline)

src = pipeparts.mkaudiotestsrc(pipeline, wave = 9)
src = pipeparts.mkcapsfilter(pipeline, src, caps = "audio/x-raw, rate=10")
src = pipeparts.mkpeak(pipeline, src, n = 10)
src = pipeparts.mknxydumpsink(pipeline, src, "/dev/stdout")

if pipeline.set_state(Gst.State.PLAYING) == Gst.StateChangeReturn.FAILURE:
        raise RuntimeError("pipeline failed to enter PLAYING state")
mainloop.run()

\end{lstlisting}
We have found that, using python to procedurally build Gstreamer graphs, we can
construct enormous pipelines containing tens of thousands of distinct elements.
A prime example of this is our workhorse signal processing pipeline used for
discovering compact binary mergers as described in the next section.

\subsection{Compact binary searches}
\label{ssec:cbc_searches}

\texttt{Makefile.softwarex\_test} provides an example of the general workflow
involved in offline gravitational wave analyses, which primarily rely on the
\texttt{gstlal} and \texttt{gstlal-inspiral} packages. The miniconda
installation of GstLAL will run this example in $\sim$30 minutes on a single
machine. Production level compact binary searches analyzing data from the three advanced interferometers, on the other hand, take $\sim$1 week when distributed over $\mathcal{O}$(1000) core computing clusters with optimized software builds. The test Makefile is
entirely self-contained; the target and dependency relationships and brief
comments within describe the workflow. Although the structure presented in this
test is linear in nature, full scale searches for compact binaries are heavily
parallelized to take advantage of DAG scheduling. The requisite commands to run
this test analysis are shown below.
\begin{lstlisting}[language=Bash]

mkdir workflow-test && cd workflow-test
export LAL_PATH=${CONDA_PREFIX} GSTLAL_FIR_WHITEN=0 TMPDIR=/tmp
wget https://dcc.ligo.org/public/0168/P2000195/004/Makefile.softwarex_test
make -f Makefile.softwarex_test

\end{lstlisting}

\section{Impact}
\label{sec:impact}

GstLAL has played an integral role in the history of gravitational wave
detections, and has participated in gravitational wave searches since
S5~\cite{wade2015}.  Although GstLAL was designed specifically for
near-real-time applications, the low-latency pipeline was prevented from
searching for binary black holes at the start of Advanced LIGO's first
observing run (O1). The pipeline was only allowed to use a template bank sensitive
to electromagnetically binary neutron stars and neutron star - black hole
binaries. The software was fully capable of detecting binary black holes in near-real-time. The LIGO
collaboration desired blind binary black hole (BBH) analyses, and since BBH
systems are not expected to produce electromagnetic radiation there was no
perceived need to detect them in near-real-time. The restriction on the allowed
template bank rendered GstLAL unable to detect GW150914 in low-latency, though
it was one of two matched filter pipelines used to analyze archival data and
validate the event~\cite{Abbott:2016blz}. GW150914 was initially detected in
low-latency by the weakly-modeled burst pipeline CWB~\cite{Klimenko:2015ypf},
which demonstrated that there could be no truly blind analysis while
low-latency burst pipelines produced alerts. As a result, in late 2015 the
GstLAL pipeline was approved to include BBHs with total mass $\lesssim 100
M_\odot$ in its low-latency configuration. GstLAL quickly demonstrated its
ability to recover BBHs in low-latency; it became the first matched-filter
pipeline using waveforms based on general relativity to make a near real-time
detection of a compact binary with the discovery of
GW151226~\cite{Abbott:2016nmj}.

Development work between Advanced LIGO's first and second observing runs
focused on enabling single detector discoveries and
incorporating data from Virgo in the analysis.
These efforts were rewarded by August 2017 as GstLAL became the first three-detector
matched-filter search in the Advanced LIGO era and, more notably, the first
(and to date, the only) gravitational-wave detection pipeline to observe a
binary neutron star merger in low-latency~\cite{GCN24168,G298048}.  Although
both Advanced LIGO interferometers and the Advanced Virgo interferometer were
operating at the time of GW170817, it was initially observed in a single
interferometer. This marked the first single detector observation of a gravitational
wave, and the autonomous identification of the candidate enabled rapid offline
follow-up of the candidate within the LIGO-Virgo collaboration.

Advanced LIGO's third observing run (O3) marked the beginning of open
public alerts (OPA)\cite{GCN24045}. For the first time, candidates with
false-alarm-rates below 1 per month\footnote{after applying a trials factor to
account for the number of concurrently running searches.} were made public at
the time of discovery. Although the first public alert was distributed in error
by the collaboration~\cite{GCN24109}, the identification of binary black hole
candidate S90408an~\cite{GCN24069} marked the successful start of the era of
automated public alerts. At first only candidates appearing in two or more
detectors were approved for automated release, but GW170817 had already demonstrated
the importance of single detector searches. Indeed, two weeks into the
third observing run GstLAL was the only pipeline to detect GW190425 in
near-real-time~\cite{GCN24168,Abbott:2020uma}, further highlighting the
necessity of single detector searches.

Two months into the observing
run, GstLAL became the only pipeline approved to release single detector
candidates as OPAs. This was a high risk, high reward endeavor. Matched filter
searches have traditionally been able to suppress the background by demanding
coincidence across interferometers; single detector candidates do not benefit
from this effect and can therefore be more susceptible to short term noise
transients. Unfortunately, this resulted in several retractions throughout O3
as the GstLAL team worked on ways to mitigate the effects of noise transients
in single detectors. By the end of the second half of the observing run,
pipeline tuning had reduced the rate of retractions.

GstLAL has contributed to all gravitational-wave
discoveries published by the LIGO Scientific
Collaboration~\cite{Abbott:2020uma,LIGOScientific:2018mvr,TheLIGOScientific:2017qsa,Abbott:2017oio,Abbott:2017gyy,Abbott:2017vtc,Abbott:2016nmj,Abbott:2016blz},
but it has also contributed to searches for as yet undetected sources.
Sub-solar mass and intermediate mass black holes both pose problems for
conventional models of stellar evolution, and GstLAL has directly contributed
to searches of
both~\cite{Abbott:2018oah,Salemi:2019ovz,Authors:2019qbw,Abbott:2019prv}.
Although these searches have not yet yielded any detections, the null results
have been able to place strict limits on the abundance of such objects and have
also provided the tightest limit to date on a primordial black hole model of
the dark matter.

\section{Conclusions}
\label{sec:conclusions}

The GstLAL library has significantly impacted the progress of
gravitational-wave astrophysics, not only via compact binary searches, but also
through contributions to detector calibration and characterization efforts. The
low-latency GstLAL based inspiral pipeline was instrumental in the first
multi-messenger discovery with gravitational waves, and strives to lead the
march towards more remarkable observations with ground-based interferometers.

\section{Conflict of interest}

No conflict of interest exists: we wish to confirm that there are no known
conflicts of interest associated with this publication and there has been no
significant financial support for this work that could have influenced its
outcome.

\section*{Acknowledgments}
\label{}

Funding for this work was provided by the National Science Foundation through
awards: PHY-1454389, OAC-1642391, PHY-1700765, OAC-1841480, PHY-1607178, and
PHY-1847350.  Funding for this project was provided by the Charles E. Kaufman
Foundation of The Pittsburgh Foundation. Computations for this research were
performed on the Pennsylvania State University’s Institute for Computational
and Data Sciences Advanced CyberInfrastructure (ICDS-ACI) and VM hosting.  We
are grateful for computational resources provided by the Leonard E Parker
Center for Gravitation, Cosmology and Astrophysics at the University of
Wisconsin-Milwaukee.  Computing support was provided by the LIGO Laboratory
through National Science Foundation grant PHY-1764464.  GstLAL relies on many
other open source software libraries; we gratefully acknowledge the development
and support of NumPy~\cite{numpy}, SciPy~\cite{scipy}, PyGTK~\cite{pygtk},
PyGST~\cite{pygst}, Bottle~\cite{bottle}, Kafka~\cite{kafka},
Fftw3F~\cite{fftw3f}, Intel MKL~\cite{intelmkl}, GLib2~\cite{glib2}, GNU
Scientific Library~\cite{gsl}, and GWpy~\cite{gwpy}.  The authors gratefully
acknowledge the LIGO-Virgo-Kagra collaboration for support, review, and
valuable critiques throughout various stages of development of the GstLAL
library. We are especially thankful for collaborations within the Compact
Binary Coalescence working group.

\bibliographystyle{elsarticle-num}