\section{Introduction} \label{sec:intro}

In the past two years evidence has accumulated pointing to discrepancies between theoretical predictions and experimental results in angular observables in the rare decay $B\to K^* \mu^+\mu^-$~\cite{Aaij:2013qta,LHCb-CONF-2015-002}. It was noted early on~\cite{Descotes-Genon:2013wba,Altmannshofer:2013foa,Beaujean:2013soa,Hurth:2013ssa} that new physics (NP) contributing to the operator,
\begin{equation}
\label{eqn:Hamiltonian}
 \mathscr{H}_\text{eff} = \Delta \mathscr{C}_9^\mu (\bar s \gamma_\alpha P_L b)(\bar \mu \gamma^\alpha \mu)~,
\end{equation}
can explain the discrepancies. This operator consists of a left-handed quark current and a vector current for the muons.
The theoretical predictions of the Standard Model (SM) involve non-perturbative QCD effects that could mimic a NP contribution to exactly such an operator. It is not easy to reliably estimate some parts of the non-perturbative QCD effects~\cite{Khodjamirian:2010vf,Jager:2012uw,Lyon:2014hpa,Descotes-Genon:2014uoa,Jager:2014rwa}. Thus, the aforementioned discrepancies might still be the result of underestimated SM contributions. Nevertheless, a variety of NP models have been proposed in the literature as the origin of the operator in Eq.~(\ref{eqn:Hamiltonian})~\cite{Descotes-Genon:2013wba,Altmannshofer:2013foa,Gauld:2013qba,Buras:2013qja,Gauld:2013qja,Datta:2013kja,Buras:2013dea,Altmannshofer:2014cfa,Biancofiore:2014wpa,Buras:2014fpa,Altmannshofer:2014rta,Gripaios:2014tna,Crivellin:2015mga,Crivellin:2015lwa,Sierra:2015fma,Celis:2015ara,Belanger:2015nma,Niehoff:2015iaa}.

In Ref.~\cite{Altmannshofer:2014cfa} we presented a model with an extra vector-boson associated with the gauging (and spontaneous breaking) of muon-number minus tau-number, $L_\mu - L_\tau$. Gauging of $L_\mu - L_\tau$ guarantees a vectorial coupling of the vector-boson to muons. The new vector-boson can couple to the different quark flavours indirectly through its coupling to very heavy colored fermions that mix with the quarks as in Ref.~\cite{Fox:2011qd}. If the mass of the new vector-boson is sufficiently heavy compared to the $B$-meson mass, then it can be integrated out of the low-energy theory to yield the operator in Eq.~(\ref{eqn:Hamiltonian}).

One clean prediction of this model, first made in Ref.~\cite{Altmannshofer:2014cfa}, is a $\sim 20\%$ suppression of the rate of the inclusive $B \to X_s \mu^+\mu^-$ decay relative to the electron mode $B \to X_s e^+e^-$ which remains SM-like. The structure of $L_\mu - L_\tau$ also predicts a corresponding $\sim 20\%$ enhancement of the decay rate of the tauonic mode with respect to its SM prediction. The predicted suppression in the muonic mode has subsequently received experimental support with the recent hints for violation of lepton flavor universality (LFU) in the exclusive decay $B \to K \ell^+\ell^-$ by LHCb~\cite{Aaij:2014ora}
\begin{eqnarray} \label{eq:RK}
 R_K &=& \frac{\text{BR}(B \to K \mu^+\mu^-)}{\text{BR}(B \to K e^+e^-)}  \\\nonumber &=& 0.745^{+0.090}_{-0.074} ~\text{(stat)}~\pm 0.036 ~\text{(syst)} ~,
\end{eqnarray}
where the branching ratios refer to a dimuon invariant mass squared, $q^2$, between 1~GeV$^2 < q^2 < 6$~GeV$^2$.
A number of analyses have appeared since, interpreting the hint of lepton flavor universality violation~(\ref{eq:RK}) in the context of new physics models~\cite{Alonso:2014csa,Hiller:2014yaa,Ghosh:2014awa,Biswas:2014gga,Hurth:2014vma,Glashow:2014iga,Hiller:2014ula,Bhattacharya:2014wla,Niehoff:2015bfa,Becirevic:2015asa,Alonso:2015sja,Greljo:2015mma,Calibbi:2015kma}.

In this short paper we wish to sharpen the prediction for lepton flavor universality violation made in Ref.~\cite{Altmannshofer:2014cfa}, including the theoretical uncertainties associated with it. Furthermore, we provide a variety of new predictions for lepton flavour universality violation in other, independent, observables.

\section{The Model} \label{sec:model}

We briefly review the defining components of the model and refer the reader to Ref.~\cite{Altmannshofer:2014cfa} for more details. We extend the SM by a new $U(1)^\prime$ gauge group with a new gauge-boson, $Z^\prime$, corresponding to muon-number minus-tau number, $L_\mu - L_\tau$. The gauge group is anomaly free with the SM particle content. The $L_\mu - L_\tau$ gauge symmetry is spontaneously broken by the vacuum expectation value $v_\Phi$ of a new heavy scalar $\Phi$ charged under $L_\mu - L_\tau$. The gauge-boson obtains a mass $m_{\zprime} = g^\prime v_\Phi$, with the $U(1)^\prime$ gauge coupling $g^\prime$. As we will see below, the $B$-meson anomalies are not very sensitive to either of these parameters, $m_{\zprime}$, $g^\prime$, or $v_\Phi$. Other observations~\cite{Altmannshofer:2014cfa, Altmannshofer:2014pba} constrain the model, but leave ample room in the parameter space defined by $m_{\zprime}$ and $g^\prime$ with $O(10) {~\rm GeV} \lesssim m_{\zprime} \lesssim \text{few} {~\rm TeV}$ with a corresponding increasing range in the allowed gauge
coupling $g^\prime$.

\begin{figure}
\includegraphics[scale=0.65]{effective_coupling_1.pdf}
\caption{Example diagram in the high energy theory that leads to flavor-changing effective couplings of the $Z^\prime$ to left-handed SM quarks. Similar diagrams involving other heavy colored fermions can lead to an effective coupling of the $Z^\prime$ to right-handed SM quarks.}
\label{fig:effective_coupling}
\end{figure}

The coupling of the $Z^\prime$ to quarks at low-energies is derived from an anomaly-free theory involving new heavy colored fermions, which are vector-like charged under the $U(1)^\prime$ and the SM gauge group. These heavy colored fermions couple to the SM quarks through Yukawa couplings with the scalar $\Phi$. The spontaneous breaking of the $U(1)^\prime$ leads to mass mixing between the SM quarks and the heavy fermions. In turn, the mass mixing results in an effective coupling of the $Z^\prime$ to SM quarks, as illustrated in Fig.~\ref{fig:effective_coupling}. The low-energy effective coupling is therefore determined by the type of heavy fermions present (electroweak singlets vs.\ doublets), their masses, and their Yukawa couplings to the SM quarks and $\Phi$. At the even lower energies scales relevant for $B$-meson decays, after the $Z^\prime$ has been integrated out, we are left with contributions to the Wilson coefficients of the semileptonic $b \to s \ell\ell$ operators with leptonic vector currents
\begin{equation} \label{eq:Heff}
\mathscr{H}_\text{eff} = \sum_{\ell = e, \mu, \tau} \Big[ \Delta \mathscr{C}_9^\ell (\bar s \gamma_\alpha P_L b) +  \Delta \mathscr{C}_9^{\prime \ell} (\bar s \gamma_\alpha P_R b) \Big] (\bar \ell \gamma^\alpha \ell) ~,
\end{equation}
with
\begin{eqnarray}
\label{eqn:def_C_9}
\Delta \mathscr{C}_9^e = 0 ~&,&~~ \Delta \mathscr{C}_9^{\mu} = - \Delta \mathscr{C}_9^{\tau} = \frac{Y_{Qb} Y_{Qs}^*}{2m_Q^2}~, \\
\label{eqn:def_Cp_9}
\Delta \mathscr{C}_9^{\prime e} = 0 ~&,&~~ \Delta \mathscr{C}_9^{\prime \mu} = - \Delta \mathscr{C}_9^{\prime \tau} = -\frac{Y_{Db} Y_{Ds}^*}{2m_D^2} ~,
\end{eqnarray}
where $m_Q$ and $m_D$ are the masses of heavy electroweak doublet quarks and heavy electroweak singlet down type quarks, respectively, and  $Y_{Qi}$, $Y_{Di}$ are their corresponding Yukawa couplings with $\Phi$, see~\cite{Altmannshofer:2014cfa}. In addition to the operator in~(\ref{eqn:Hamiltonian}), also its right-handed counterpart is generated in general.
It is common to use dimensionless Wilson coefficients $\Delta C_9^{(\prime)}$ defined as
\begin{equation}
 - \frac{4G_F}{\sqrt{2}} V_{tb} V_{ts}^* \frac{e^2}{16\pi^2} \Delta C_9^{(\prime)} = \Delta \mathscr{C}_9^{(\prime)} ~.
\end{equation}
As alluded to earlier, interestingly the Wilson coefficients do not depend on $m_{\zprime}$, $g^\prime$, or $v_\Phi$, but only on the properties of the heavy fermions and their Yukawa couplings to SM quarks. This is so because the $m_{\zprime}^2$ that would  appear in the denominator from integrating out the $Z^\prime$, is exactly canceled by the product of $(g^\prime)^2$ (coming from each vertex) and $v_\Phi^2$ (coming from the mass mixing) in the numerator. In addition, the $L_\mu - L_\tau$ gauge symmetry forces $\Delta C_9^e = \Delta C_9^{\prime e}= 0$  and opposite sign between the muonic and tauonic contributions. We are only free to adjust $\Delta C_9^\mu$ and $\Delta C_9^{\prime \mu}$ in Eqs.~(\ref{eqn:def_C_9}) and~(\ref{eqn:def_Cp_9}) through the choice of Yukawa coupling and heavy fermion masses. This tight structure is the origin of the sharp prediction for the presence of LFU violations in $B$-meson decay: once $\Delta C_9^\mu$ and $\Delta C_9^{\prime \mu}$ are fixed to explain the angular anomalies in $B\to K^* \mu^+\
mu^-$, they immediately effect the branching ratio of $B\to K \mu^+\mu^-$ and various other observables in rare $B$-meson decays based on the $b \to s \mu\mu$ transition. All electron modes, $b \to s ee$, remain SM-like. In the following we will drop the lepton flavor superscript on the Wilson coefficients and denote the muon specific coefficients with $\Delta C_9$ and $\Delta C_9^\prime$.

\section{Probes of Lepton Flavor Universality}
\label{sec:LFU}

We consider observables in various decays based on the quark level $b \to s \ell\ell$ transition that are affected by the operators in Eq.~(\ref{eq:Heff}): the exclusive decays $B \to K \ell^+\ell^-$, $B \to K^* \ell^+\ell^-$, and $B_s \to \phi \ell^+\ell^-$, as well as the inclusive $B \to X_s \ell^+\ell^-$ decay.

\subsection{Exclusive Decays}
\label{sec:exc}

We consider the following ratios of branching ratios that are tests of lepton flavor universality~\cite{Hiller:2003js,Bobeth:2007dw,Hiller:2014ula}
\begin{eqnarray}
R_K &=& \frac{\text{BR}(B \to K \mu^+\mu^-)}{\text{BR}(B \to K e^+e^-)}  ~, \\
R_{K^*} &=& \frac{\text{BR}(B \to K^* \mu^+\mu^-)}{\text{BR}(B \to K^* e^+e^-)}  ~, \\
R_\phi &=& \frac{\text{BR}(B \to \phi \mu^+\mu^-)}{\text{BR}(B \to \phi e^+e^-)}  ~.
\end{eqnarray}
In the SM, these ratios are expected to be unity to a very high accuracy.\footnote{Radiation of collinear photons off the leptons can lead to logarithmically enhanced QED corrections that can be different in the electron and muon modes. In the $R_K$ measurement performed by LHCb~\cite{Aaij:2014ora}, such effects are included in the event simulation.}
In the following, superscripts on the LFU observables indicate the bin of di-lepton invariant mass in which the observable is measured. For example, $R_K^{[1,6]}$ refers to the $q^2$ bin $1\,$GeV$^2 < q^2 < 6\,$GeV$^2$, etc.

In the $L_\mu - L_\tau$ framework, the above LFU ratios can be expressed in terms of the muon specific Wilson coefficients $\Delta C_9$ and $\Delta C_9^\prime$. For all three decays we consider one bin at low and high $q^2$ each. We provide explicit expressions in Appendix~\ref{sec:appendix1}.

The ratios $R_\phi$ and $R_{K^*}$ depend on the left-handed and right-handed Wilson coefficients $\Delta C_9$ and $\Delta C_9^\prime$ in a very similar way. Differences only arise from $SU(3)$ flavor breaking and the rather sizable life time difference in the $B_s$-meson system. However, given the expected experimental sensitivities in the near future, one has to an excellent approximation $R_\phi \simeq R_{K^*}$. Comparing $R_\phi$ and $R_{K^*}$ with $R_K$ allows in principle to disentangle contributions from left-handed and right-handed quark currents~\cite{Hiller:2014ula}.

In addition to the branching ratio ratios, we also consider lepton flavor universality tests using angular observables in the $B \to K^* \ell^+\ell^-$ decay. As some angular observables have zero crossings we do not use lepton flavor ratios but rather construct lepton flavor differences
\begin{eqnarray}
D_{A_\text{FB}} &=& A_\text{FB}(B \to K^* \mu\mu) - A_\text{FB}(B \to K^* ee) ~, \\
D_{F_L} &=& F_L(B \to K^* \mu\mu) - F_L(B \to K^* ee) ~, \\
D_{S_5} &=& S_5(B \to K^* \mu\mu) - S_5(B \to K^* ee) ~.
\end{eqnarray}
Here, $A_\text{FB}$ is the lepton forward-backward asymmetry, $F_L$ is the longitudinal polarization fraction of the $K^*$ vector meson, and the angular observable $S_5$ is related to the ``$B \to K^* \mu^+\mu^-$ anomaly'' in the observable $P_5^\prime = S_5/\sqrt{F_L (1-F_L)}$~\cite{Aaij:2013qta}.
Other probes of lepton flavor universality have been considered in the recent literature. For instance: lepton flavor ratios of angular observables in $B \to K^* \ell^+ \ell^-$~\cite{Altmannshofer:2014rta,Jager:2014rwa} and lepton flavor specific shifts of the zero crossings of the angular observables $A_\text{FB}$ and $S_5$~\cite{Jager:2014rwa}.

In contrast to $B \to K^* \ell^+\ell^-$, the $B_s \to \phi \ell^+\ell^-$ decay is not self-tagging. Therefore, the observables $D_{A_\text{FB}}$ and $D_{S_5}$ cannot be easily accessed in the $B_s \to \phi \ell^+\ell^-$ decay.
We therefore concentrate on $B \to K^* \ell^+\ell^-$.

At high $q^2$, the angular observables $A_\text{FB}$, $F_L$, and $S_5$ do not differ significantly from their SM predictions in regions of NP parameter space that provide a good fit of the anomalies~\cite{Altmannshofer:2014rta}. Thus, we consider the corresponding LFU differences only in the low $q^2$ region.
Expressions for the new $B \to K^* \ell^+\ell^-$ LFU differences in terms of the Wilson coefficients $\Delta C_9$ and $\Delta C_9^\prime$ are given in Appendix~\ref{sec:appendix1}.

The new observables provide additional means to test lepton flavor universality and, in the case that evidence for violation of lepton flavor universality can be established, can help to distinguish between left-handed and right-handed currents.

\subsection{Inclusive Decay}
\label{sec:inc}

The inclusive decay $B \to X_s \ell^+\ell^-$ gives access to three independent observables, $H_L$, $H_T$, and $H_A$~\cite{Lee:2006gs}. We will consider the following combinations of these observables:
(a) the branching ratio BR$(B \to X_s \ell^+\ell^-) = H_L + H_T$;
(b) the normalized lepton forward-backward asymmetry $A_\text{FB} = \frac{3}{4} H_A/(H_L + H_T)$;
(c) the observable $F_L = H_L/(H_L + H_T)$, in analogy to the longitudinal polarization fraction of the $K^*$ in the $B \to K^* \ell^+ \ell^-$ decay.
In slight abuse of notation, we denote these observables in the same way in both the $B \to K^* \ell^+ \ell^-$ decay and the $B \to X_s \ell^+ \ell^-$ decay. Analogous to the exclusive decays we will consider a LFU ratio for the branching ratios, and LFU differences for the angular observables $A_\text{FB}$ and $F_L$
\begin{eqnarray}
R_{X_s} &=& \frac{\text{BR}(B \to X_s \mu\mu)}{\text{BR}(B \to X_s ee)} ~, \\
D_{F_L^{X_s}} &=& F_L(B \to X_s \mu\mu) - F_L(B \to X_s ee) ~, \\
D_{A_\text{FB}^{X_s}} &=& A_\text{FB}(B \to X_s \mu\mu) - A_\text{FB}(B \to X_s ee)  ~.
\end{eqnarray}
Expressions for these observables in terms of the Wilson coefficients $\Delta C_9$ and $\Delta C_9^\prime$ are collected in Appendix~\ref{sec:appendix1}.
In contrast to the exclusive decays, the right-handed coefficient $\Delta C_9^\prime$ enters these expressions only at the quadratic level.
The LFU observables in the inclusive decay are therefore complementary to those in the exclusive decays discussed above.

\section{Predictions for Lepton Flavor Universality Violation}
\label{sec:predictions}

\begin{figure*}[th]
\centering
\includegraphics[width=0.5\textwidth]{bands_gen.pdf}\includegraphics[width=0.5\textwidth]{bands_MFV.pdf}
\caption{Predictions for lepton flavor universality ratios and differences in models with gauged $L_\mu - L_\tau$. Left: generic quark couplings. Right: MFV quark couplings. The solid lines show the predicted central values, the $1\sigma$ and $2\sigma$ ranges are shown by the shaded bands. The SM predictions are indicated by the dotted line. The experimental measurement of $R_K^{[1,6]}$ is also given for comparison.}
\label{fig:LFU}
\end{figure*}

\renewcommand{\arraystretch}{2.4}
\setlength\tabcolsep{0pt}
\begin{table*}[th]
\begin{center}
\begin{tabularx}{\textwidth}{ccccccc}
\hline\hline
 & \multicolumn{2}{c}{(i) $L_\mu - L_\tau$ generic} & \multicolumn{2}{c}{(ii) $L_\mu - L_\tau$ MFV} & \multicolumn{2}{c}{(iii) LH quarks \& muons}  \\ \hline
 \,~~~~~~~~~~~~~~~~~~~~~~~&~~~~~~~~~~ $1\sigma$ ~~~~~~~~~~&~~~~~~~~~~ $2\sigma$ ~~~~~~~~~~&~~~~~~~~~~ $1\sigma$ ~~~~~~~~~~&~~~~~~~~~~ $2\sigma$ ~~~~~~~~~~&~~~~~~~~~~ $1\sigma$ ~~~~~~~~~~&~~~~~~~~~~ $2\sigma$ \\
\hline\hline
\rowcolor[gray]{.9} $R_K^{[1,6]}$ & $0.85^{+0.09}_{-0.08}$ & $[0.66,1.02]$ & $0.77^{+0.08}_{-0.06}$ & $[0.63,0.92]$ & $0.76^{+0.10}_{-0.09}$ & $[0.55,0.93]$   \\
$R_K^{[15,22]}$ & $0.86^{+0.09}_{-0.08}$ & $[0.67,1.02]$ & $0.78^{+0.08}_{-0.06}$ & $[0.63,0.91]$ & $0.76^{+0.10}_{-0.09}$ & $[0.55,0.93]$   \\
\rowcolor[gray]{.9} $R_{K^*}^{[1,6]}$ & $0.77^{+0.07}_{-0.06}$ & $[0.64,0.91]$ & $0.83^{+0.06}_{-0.04}$ & $[0.72, 0.94]$ & $0.77^{+0.10}_{-0.08}$ & $[0.58,0.93]$   \\
$R_{K^*}^{[15,19]}$ & $0.72^{+0.09}_{-0.07}$ & $[0.57,0.89]$ & $0.78^{+0.08}_{-0.06}$ & $[0.64,0.92]$ & $0.75^{+0.11}_{-0.09}$ & $[0.53,0.92]$ \\
\rowcolor[gray]{.9} $R_{\phi}^{[1,6]}$ & $0.77^{+0.08}_{-0.06}$ & $[0.64,0.92]$ & $0.83^{+0.06}_{-0.05}$ & $[0.73, 0.94]$ & $0.77^{+0.10}_{-0.08}$ & $[0.57,0.93]$   \\
$R_{\phi}^{[16,19]}$ & $0.72^{+0.09}_{-0.07}$ & $[0.55,0.89]$ & $0.78^{+0.08}_{-0.06}$ & $[0.64,0.91]$ & $0.75^{+0.10}_{-0.09}$ & $[0.54,0.93]$   \\
\rowcolor[gray]{.9} $R_{X_s}^{[1,6]}$ & $0.82^{+0.04}_{-0.03}$ & $[0.75,0.89]$ & $0.82^{+0.04}_{-0.03}$ & $[0.75,0.89]$ & $0.73^{+0.08}_{-0.08}$ & $[0.56,0.88]$   \\
$R_{X_s}^{>14.4}$ & $0.90^{+0.06}_{-0.05}$ & $[0.79,1.02]$ & $0.90^{+0.06}_{-0.05}$ & $[0.78,1.02]$ & $0.84^{+0.11}_{-0.11}$ & $[0.62,1.04]$    \\
\hline\hline
\rowcolor[gray]{.9} $D_{P_5^\prime}^{[1,6]}$ & $+0.27^{+0.07}_{-0.13}$ & $[+0.07,+0.47]$ & $+0.27^{+0.07}_{-0.12}$ & $[+0.08,+0.45]$ & $+0.112^{+0.049}_{-0.093}$ & $[-0.014,+0.281]$    \\
$D_{S_5}^{[1,6]}$ & $+0.109^{+0.030}_{-0.050}$ & $[+0.029,+0.189]$ & $+0.110^{+0.026}_{-0.050}$ & $[+0.032,+0.184]$ & $+0.045^{+0.016}_{-0.026}$ & $[+0.008,+0.093]$  \\
\rowcolor[gray]{.9} $D_{A_\text{FB}}^{[1,6]}$ & $-0.082^{+0.031}_{-0.025}$ & $[-0.139,-0.029]$ & $-0.074^{+0.031}_{-0.019}$ & $[-0.125,-0.025]$ & $-0.034^{+0.016}_{-0.015}$ & $[-0.071,-0.007]$   \\
$D_{A_\text{FB}^{X_s}}^{[1,3.5]}$ & $-0.070^{+0.017}_{-0.014}$ & $[-0.102,-0.040]$ & $-0.070^{+0.018}_{-0.015}$ & $[-0.103,-0.037]$ & $-0.043^{+0.014}_{-0.012}$ & $[-0.073,-0.018]$   \\
\rowcolor[gray]{.9} $D_{A_\text{FB}^{X_s}}^{[3.5,6]}$ & $-0.081^{+0.023}_{-0.020}$ & $[-0.124,-0.039]$ & $-0.079^{+0.026}_{-0.020}$ & $[-0.124,-0.033]$ & $-0.035^{+0.013}_{-0.012}$ & $[-0.066,-0.013]$    \\
$D_{F_L}^{[1,6]}$ & $-0.075^{+0.030}_{-0.022}$ & $[-0.128,-0.024]$ & $-0.054^{+0.022}_{-0.015}$ & $[-0.091,-0.018]$ & $-0.017^{+0.010}_{-0.007}$ & $[-0.038,-0.003]$  \\
\rowcolor[gray]{.9} $D_{F_L^{X_s}}^{[1,3.5]}$ & $+0.014^{+0.020}_{-0.017}$ & $[-0.024,+0.051]$ & $+0.014^{+0.022}_{-0.018}$ & $[-0.024,+0.052]$ & $+0.041^{+0.020}_{-0.015}$ & $[-0.001,+0.070]$    \\
$D_{F_L^{X_s}}^{[3.5,6]}$ & $+0.033^{+0.010}_{-0.009}$ & $[+0.014,+0.051]$ & $+0.035^{+0.011}_{-0.008}$ & $[+0.015,+0.053]$ & $+0.065^{+0.003}_{-0.004}$ & $[+0.056,+0.073]$   \\
\hline\hline
\end{tabularx}
\end{center}
\caption{\small Predictions for lepton flavor universality ratios and differences in models with gauged $L_\mu - L_\tau$ and in models with purely left-handed quark and muon currents.}
\label{tab:summary}
\end{table*}

Having defined the observables that probe lepton flavor universality, we go on and provide new physics predictions for these observables. In the context of $L_\mu - L_\tau$, we consider two cases:
\begin{itemize}
 \item[(i)] the generic case where the mixing Yukawas are arbitrary and both $\Delta C_9$ and $\Delta C_9^\prime$ can be affected significantly by new physics;
 \item[(ii)] the physically motivated case of minimal flavor violation (MFV), where the flavor structure of the left-handed mixing Yukawas are determined by CKM angles, $Y_{Qb} Y_{Qs}^* \propto V_{tb} V_{ts}^*$, and the right-handed mixing Yukawas are negligible. In that case only $C_9$ receives new physics contributions, while $\Delta C_9^\prime \sim 0$.
\end{itemize}

In~\cite{Altmannshofer:2014rta} the following best fit values for $\Delta C_9^{(\prime)}$ and improvements in the $\chi^2$ were found in a global fit to all relevant $ b \to s \mu^+\mu^- $ data
\begin{align} \label{eq:fit}
 \text{(i)}~ & \Delta C_9 = -1.10 \,,~ \Delta C_9^\prime = 0.45 \,, & \Delta \chi^2 = 15.6 \,, \\
 \text{(ii)}~ & \Delta C_9 = -1.07  \,, & \Delta \chi^2 = 13.7 \,.
\end{align}
We use the corresponding best fit ranges for $\Delta C_9^{(\prime)}$ determined in~\cite{Altmannshofer:2014rta} to make predictions for the LFU observables defined in Sec.~\ref{sec:LFU}.

The results are summarized in Fig.~\ref{fig:LFU} and Tab.~\ref{tab:summary}.
Shown are the predicted central values together with the $1\sigma$ and $2\sigma$ ranges. In Fig.~\ref{fig:LFU} also the SM predictions and the experimental result for $R_K^{[1,6]}$ are shown for comparison.
For completeness we also include predictions for the lepton flavor difference of the angular observable $P_5^\prime$
\begin{equation}
 D_{P_5^\prime} = P_5^\prime(B \to K^* \mu\mu) - P_5^\prime(B \to K^* ee) ~,
\end{equation}
with $P_5^\prime = S_5 /\sqrt{F_L(1-F_L)}$.

We observe that both the LFU ratios, $R_i$, and the LFU differences, $D_i$, differ significantly from the respective SM predictions.
All branching ratio ratios that we consider are below 1, with central values ranging approximately between $\sim 0.75$ and  $\sim 0.9$.
The predictions for $D_{P_5^\prime}$ and $D_{S_5}$ are above their SM values, while the various differences of $A_\text{FB}$ and $F_L$ are predicted below their SM values.
We expect that within the next few years of LHCb data taking the experimental sensitivity should be high enough to either confirm or exclude the predicted non-standard values for $R_K$ and $R_{K^*}$.
Given that $B_s$-meson production is suppressed compared to $B$-meson production at the LHC, establishing a deviation of $R_\phi$ from 1 will be more challenging.
At Belle II, the predicted non-standard values for $R_{X_s}$ should be well within reach.
Concerning the LFU differences, it has to be seen if experimental sensitivities at LHCb and Belle II will be sufficient to observe also the predicted deviations in $D_{P_5^\prime}$, $D_{S_5}$, $D_{A_\text{FB}}$, and $D_{F_L}$ in the near future.

We now comment on the differences between the generic case (i) and the MFV case (ii).
Notable differences do exist in the predictions for the ratios $R_K$, $R_{K^*}$ and $R_\phi$.
Because of the parity of the final state mesons, the ratios are corrected by the combination $C_9+C_9^\prime$ ($R_K$) or mainly by $C_9-C_9^\prime$ ($R_{K^*}$ and $R_\phi$). As the global fit results~(\ref{eq:fit}) show preference for a negative value for $C_9$ and a slightly positive value for $C_9^\prime$, this leads to the MFV prediction to be consistently higher (lower) than the generic prediction for $ R_{K^*}$ and $R_{\phi}$ ($R_K$).
As the inclusive decay $B \to X_s \ell^+\ell^-$ depends on the $\Delta C_9^\prime$ only at the quadratic level, the predictions for the corresponding LFU observables are virtually identical in the case (i) and case (ii).
A $\sim 20\%$ suppression of $R_{X_s}$ compared to the SM prediction is therefore a very robust prediction of the $L_\mu - L_\tau$ scenario.
We observe that also the predictions for the $B \to K^* \ell^+\ell^-$ angular differences $D_{P_5^\prime}$, $D_{S_5}$, and $D_{A_\text{FB}}$, are very similar in the cases (i) and (ii).

\begin{figure}[th]
\centering
\includegraphics[width=0.49\textwidth]{bands_LH.pdf}
\caption{Predictions for lepton flavor universality ratios and differences in models with purely left-handed quark and muon currents. The solid lines show the predicted central values, the $1\sigma$ and $2\sigma$ ranges are shown by the shaded bands. The SM predictions are indicated by the dotted line. The experimental measurement of $R_K^{[1,6]}$ is also given for comparison.}
\label{fig:LFU_LH}
\end{figure}

Aside from comparison with the SM predictions, the lepton flavor universality observables can also be useful in distinguishing between alternative models of new physics.
Here, we investigate to what extent the $L_\mu - L_\tau$ setup can be distinguished from a scenario where the anomalies in $B \to K^* \mu^+\mu^-$ are resolved by NP in an operator that consists of left-handed quark and muon currents
\begin{equation}
 \mathscr{H}_\text{eff} =  - \frac{4G_F}{\sqrt{2}} V_{tb} V_{ts}^* \frac{e^2}{16\pi^2} \Delta C_L (\bar s \gamma_\alpha P_L b) (\bar \mu \gamma^\alpha P_L \mu) ~.
\end{equation}
Various NP scenarios of this type have been discussed recently~\cite{Glashow:2014iga,Bhattacharya:2014wla,Niehoff:2015bfa,Greljo:2015mma}.
Indeed, compared to scenarios with NP in $\Delta C_9$, NP in $\Delta C_L$ gives only slightly worse fits to the data~\cite{Altmannshofer:2014rta}
\begin{align} \label{eq:fitLH}
 \Delta C_L = -1.06  ~,~~ \Delta \chi^2 = 9.8 ~.
\end{align}
Expressions for the LFU observables in such a scenario are given in Appendix~\ref{sec:appendix2}.

In Fig.~\ref{fig:LFU_LH} we summarize the predicted ranges for the LFU ratios and LFU differences in this setup.
The results are also collected in the last columns of Tab.~\ref{tab:summary}.
The ranges for the LFU ratios, $R_i$, almost entirely include the ranges found in the cases (i) and (ii) discussed above. This is due to the larger range for $\Delta C_L$ that is allowed by the global fit~\cite{Altmannshofer:2014rta}.
We observe that the non-standard effects in the LFU differences, $D_i$, are considerably smaller compared to the cases (i) and (ii). Precision measurements of these observables might allow to distinguish the $L_\mu - L_\tau$ scenarios from NP purely in left-handed currents.

\section{Summary} \label{sec:sum}

In this short paper we provided explicit expressions and numerical predictions for a variety of lepton flavor universality observables associated with the model of Ref.~\cite{Altmannshofer:2014cfa} based on gauged $L_\mu - L_\tau$. We considered known LFU ratios of branching ratios as well as LFU differences of angular observables.

Generically, a robust prediction is that all branching ratio ratios are suppressed compared to the SM prediction by $O(20\%)$, both in the case of generic quark couplings and in the case of quark couplings that follow the principle of minimal flavor violation. The predictions for the $B \to K \ell^+\ell^-$ decay, $R_K =0.85^{+0.09}_{-0.08}$ (generic quark couplings) and $R_K = 0.77^{+0.08}_{-0.06}$ (MFV quark couplings), have recently received experimental support from the measurement of this ratio by LHCb~\cite{Aaij:2014ora}. The LFU ratios represent a very clear target of opportunity for future measurements in the search for new physics.

The LFU differences of angular observables provide additional means to probe new physics in rare $B$ decays. We predict non-standard effects at the level of $O(10\%)$ in the context of gauged $L_\mu - L_\tau$.
These observables might also serve as discriminants between alternative new physics models if precision measurements will become available in the future.

\begin{acknowledgments}
Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research and Innovation. IY is supported in part by funds from the NSERC of Canada.
\end{acknowledgments}

\begin{appendix}\onecolumngrid
\section{Numerical Expressions in the \texorpdfstring{\boldmath $L_\mu - L_\tau$}{Lmu - Ltau} Scenario} \label{sec:appendix1}

For the branching ratio ratios of the exlcusive decays we find
\begin{eqnarray} \label{eq:LFU1}
 R_K^{[1,6]} &\simeq& 1.00 + 0.243 (\Delta C_9 + \Delta C_9^\prime) + 0.030 (\Delta C_9 + \Delta C_9^\prime)^2 ~, \\
 R_K^{[15,22]} &\simeq& 1.00 + 0.241 (\Delta C_9 + \Delta C_9^\prime) + 0.030 (\Delta C_9 + \Delta C_9^\prime)^2 ~, \\[10pt]
 R_{K^*}^{[1,6]} &\simeq& 1.00 + 0.192 \Delta C_9 - 0.198 \Delta C_9^\prime + 0.034 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right) - 0.052 \Delta C_9 \Delta C_9^\prime~, \\
 R_{K^*}^{[15,19]} &\simeq& 1.00 + 0.241 \Delta C_9 - 0.182 \Delta C_9^\prime + 0.033 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right) - 0.049 \Delta C_9 \Delta C_9^\prime~, \\[10pt]
 R_{\phi}^{[1,6]} &\simeq& 1.00 + 0.193 \Delta C_9 - 0.200 \Delta C_9^\prime + 0.034 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right) - 0.052 \Delta C_9 \Delta C_9^\prime~, \\
 R_{\phi}^{[16,19]} &\simeq& 1.00 + 0.241 \Delta C_9 - 0.201 \Delta C_9^\prime +0.032 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right) - 0.053 \Delta C_9 \Delta C_9^\prime ~, \label{eq:LFU6}
\end{eqnarray}
in reasonable agreement with~\cite{Hiller:2014ula}.
The expressions hold for real Wilson coefficients $\Delta C_9$ and $\Delta C_9^\prime$. A study of CP violating scenarios is beyond the scope of the present work.
The theoretical uncertainty on the leading piece (unity) is expected to be at the sub-percent level, whereas we estimate that the uncertainty on the remaining numerical coefficients in Eqs.~(\ref{eq:LFU1})-(\ref{eq:LFU6}) is at the $10\%$ level. In our numerical analysis in Sec.~\ref{sec:predictions} we neglect the uncertainty of the leading piece and use uncorrelated $20\%$ relative uncertainties for the remaining numerical coefficients.
These uncertainties have only minor impact on our results, as the dominant uncertainty arises from the allowed ranges of  $\Delta C_9$ and $\Delta C_9^\prime$.

For the differences of angular observables in $B \to K^* \ell^+\ell^-$ we obtain
\begin{eqnarray} \label{eq:LFU7}
 D_{A_\text{FB}}^{[1,6]} &\simeq& \frac{0.00 + 0.057 \Delta C_9 }{ 1 + 0.192 \Delta C_9 - 0.198 \Delta C_9^\prime + 0.034 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right) - 0.052 \Delta C_9 \Delta C_9^\prime}  ~, \\[10pt]
 D_{F_L}^{[1,6]} &\simeq& \frac{0.00 + 0.039  \Delta C_9 - 0.035 \Delta C_9^\prime - 0.002 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right) - 0.008 \Delta C_9 \Delta C_9^\prime}{ 1 + 0.192 \Delta C_9 - 0.198 \Delta C_9^\prime + 0.034 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right) - 0.052 \Delta C_9 \Delta C_9^\prime} ~, \\[10pt]
 D_{S_5}^{[1,6]} &\simeq& \frac{0.00 - 0.079 \Delta C_9 - 0.033 \Delta C_9^\prime + 0.006 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right)-0.009 \Delta C_9 \Delta C_9^\prime}{ 1 + 0.192 \Delta C_9 - 0.198 \Delta C_9^\prime + 0.034 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right) - 0.052 \Delta C_9 \Delta C_9^\prime} ~, \label{eq:LFU9}
\end{eqnarray}
valid for real Wilson coefficients.
Similar to the case of the LFU ratios in~(\ref{eq:LFU1})-(\ref{eq:LFU6}), the leading piece (zero) can be predicted with very high accuracy. The uncertainty on the remaining numerical coefficients is expected to be at the $10\%$ level. In the numerical analysis, we use uncorrelated $20\%$ relative uncertainties for the corresponding numerical coefficients in~(\ref{eq:LFU7})-(\ref{eq:LFU9}).

For the LFU observables in the inclusive decay we make use of the expressions in~\cite{Huber:2015sra,Huber:2005ig} and arrive at
\begin{eqnarray} \label{eq:LFU10}
 R_{X_s}^{[1,6]} &\simeq& (R_{X_s}^{[1,6]})_\text{SM} + 0.181 \Delta C_9 + 0.036 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right)~, \\
 R_{X_s}^{>14.4} &\simeq& (R_{X_s}^{>14.4})_\text{SM} + 0.281 \Delta C_9 + 0.042 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right) ~, \\[10pt]
 D_{F_L^{X_s}}^{[1,3.5]} &\simeq& \frac{(D_{F_L}^{[1,3.5]})_\text{SM} + 0.065 \Delta C_9 + 0.003 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right)}{1 + 0.167 \Delta C_9 + 0.035 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right)} ~, \\
 D_{F_L^{X_s}}^{[3.5,6]} &\simeq& \frac{(D_{F_L}^{[3.5,6]})_\text{SM} + 0.031 \Delta C_9 - 0.001 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right)}{1 + 0.212 \Delta C_9 + 0.039 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right)} ~, \\[10pt]
 D_{A_\text{FB}^{X_s}}^{[1,3.5]} &\simeq& \frac{(D_{A_\text{FB}}^{[1,3.5]})_\text{SM} + 0.050 \Delta C_9 + 0.003 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right) }{1 + 0.167 \Delta C_9 + 0.035 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right)}~, \\
 D_{A_\text{FB}^{X_s}}^{[3.5,6]} &\simeq& \frac{(D_{A_\text{FB}}^{[3.5,6]})_\text{SM} + 0.052 \Delta C_9 - 0.003 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right) }{1 + 0.212 \Delta C_9 + 0.039 \left( (\Delta C_9)^2 + (\Delta C_9^\prime)^2 \right)}~, \label{eq:LFU16}
\end{eqnarray}
valid for real Wilson coefficients.
For the SM predictions we find from~\cite{Huber:2015sra}
\begin{eqnarray}\label{eq:SM1}
(R_{X_s}^{[1,6]})_\text{SM} = 0.970 ~&,&~ (R_{X_s}^{>14.4})_\text{SM} = 1.151~,\\
(D_{F_L^{X_s}}^{[1,3.5]})_\text{SM} = 0.079 ~&,&~ (D_{F_L^{X_s}}^{[3.5,6]})_\text{SM} = 0.062 ~, \\
(D_{A_\text{FB}^{X_s}}^{[1,3.5]})_\text{SM} = -0.010 ~&,&~ (D_{A_\text{FB}^{X_s}}^{[3.5,6]})_\text{SM} = -0.005 ~.\label{eq:SM3}
\end{eqnarray}
We expect uncertainties of $O(0.01)$ for these SM predictions and we neglect them in our numerical analysis.
The small differences in the SM predictions of $R_{X_s}$ from unity and of $D_{F_L}$ and $D_{A_\text{FB}}$ from zero, arise from logarithmically enhanced QED corrections.
Given the good theoretical control in the prediction of the inclusive $B \to X_s \ell^+\ell^-$ decay~\cite{Huber:2015sra}, the uncertainties on the remaining numerical coefficients in Eqs.~(\ref{eq:LFU10})-(\ref{eq:LFU16}) are likely below $10\%$. We use uncorrelated $10\%$ relative uncertainties in our numerical analysis. This choice of uncertainties has only very minor impact on our results.

\section{Numerical Expressions in the Scenario with Left-Handed Currents} \label{sec:appendix2}

In the scenario with NP in left-handed currents only, we find for the LFU observables in exclusive decays
\begin{eqnarray}
 R_K^{[1,6]} &\simeq& 1.00 + 0.245 \Delta C_L + 0.015 (\Delta C_L)^2 ~, \\
 R_K^{[15,22]} &\simeq& 1.00 + 0.244 \Delta C_L + 0.015 (\Delta C_L)^2 ~, \\[10pt]
 R_{K^*}^{[1,6]} &\simeq& 1.00 + 0.235 \Delta C_L + 0.017 (\Delta C_L)^2 ~, \\
 R_{K^*}^{[15,19]} &\simeq& 1.00 + 0.254 \Delta C_L + 0.016 (\Delta C_L)^2~, \\[10pt]
 R_{\phi}^{[1,6]} &\simeq& 1.00 + 0.234 \Delta C_L + 0.017 (\Delta C_L)^2~, \\
 R_{\phi}^{[16,19]} &\simeq& 1.00 + 0.253 \Delta C_L + 0.016 (\Delta C_L)^2~,
\end{eqnarray}
\begin{eqnarray}
 D_{A_\text{FB}}^{[1,6]} &\simeq& \frac{0.00 + 0.029 \Delta C_L + 0.003 (\Delta C_L)^2}{ 1 + 0.235 \Delta C_L  + 0.017 (\Delta C_L)^2}  ~, \\[10pt]
 D_{F_L}^{[1,6]} &\simeq& \frac{0.00 + 0.011 \Delta C_L - 0.001 (\Delta C_L)^2}{ 1 + 0.235 \Delta C_L  + 0.017 (\Delta C_L)^2}  ~, \\[10pt]
 D_{S_5}^{[1,6]} &\simeq& \frac{0.00 - 0.037 \Delta C_L -0.004 (\Delta C_L)^2}{ 1 + 0.235 \Delta C_L  + 0.017 (\Delta C_L)^2}  ~.
\end{eqnarray}
Our treatment of uncertainties in the numerical analysis is equivalent to the $L_\mu - L_\tau$ scenario discussed above.

For the LFU observables in the inclusive decay we find from~\cite{Huber:2015sra,Huber:2005ig}
\begin{eqnarray}
 R_{X_s}^{[1,6]} &\simeq& (R_{X_s}^{[1,6]})_\text{SM} + 0.243 \Delta C_L + 0.018 (\Delta C_L)^2 ~, \\
 R_{X_s}^{>14.4} &\simeq& (R_{X_s}^{>14.4})_\text{SM} + 0.319 \Delta C_L + 0.022 (\Delta C_L)^2 ~, \\[10pt]
 D_{F_L^{X_s}}^{[1,3.5]} &\simeq& \frac{(D_{F_L}^{[1,3.5]})_\text{SM} + 0.046 \Delta C_L + 0.002 (\Delta C_L)^2}{1 + 0.236 \Delta C_L + 0.018 (\Delta C_L)^2} ~, \\
 D_{F_L^{X_s}}^{[3.5,6]} &\simeq& \frac{(D_{F_L}^{[3.5,6]})_\text{SM} + 0.014 \Delta C_L - 0.0002 (\Delta C_L)^2}{1 + 0.272 \Delta C_L + 0.020 (\Delta C_L)^2} ~, \\[10pt]
 D_{A_\text{FB}^{X_s}}^{[1,3.5]} &\simeq& \frac{(D_{A_\text{FB}}^{[1,3.5]})_\text{SM} + 0.026 \Delta C_L + 0.004 (\Delta C_L)^2}{1 + 0.236 \Delta C_L + 0.018 (\Delta C_L)^2}~, \\
 D_{A_\text{FB}^{X_s}}^{[3.5,6]} &\simeq& \frac{(D_{A_\text{FB}}^{[3.5,6]})_\text{SM} + 0.022 \Delta C_L + 0.003 (\Delta C_L)^2 }{1 + 0.272 \Delta C_L + 0.020 (\Delta C_L)^2}~.
\end{eqnarray}
The SM predictions have already been given above in Eqs.~(\ref{eq:SM1}) - (\ref{eq:SM3}), and our treatment of uncertainties in the numerical analysis is equivalent to the $L_\mu - L_\tau$ scenario.

\end{appendix}
\twocolumngrid


\section{Introduction} \label{sec:intro}

For nearly a century, absorption features known as the diffuse interstellar bands (DIBs) have been
observed superimposed on the spectra of stars, galaxies, and quasars along numerous lines of sight through the Milky Way
and other galaxies \citep[e.g.,][]{Heger_1922_stationarylines,Cox_2011_DIBssummary}.
The current strongest consensus is that most (if not all) of the distinct features are produced by
carbon-dominated molecules, such as polycyclic aromatic hydrocarbons (PAHs) or fullerenes
\citep[e.g.,][]{Herbig_1995_dibs,Sarre_2006_DIBssummary},
but so far, definitive association of any specific species with a specific DIB has proven elusive.

The vast majority of, perhaps all, DIB detections have occurred in the
diffuse interstellar medium (ISM), unassociated with stars \citep[][]{Seab_1995_DCBs}.
However, circumstellar environments are particularly interesting potential hosts of DIB carriers,
because often different types of information concerning their temperature, density, and chemistry are available than for
a particular patch of diffuse interstellar cloud at a generally unknown distance.
Furthermore, as circumstellar environments are typically significantly warmer and denser than diffuse ISM clouds, they
provide additional, complementary regimes in which to characterize DIB behavior and thus their carriers.
Searches for DIB features in circumstellar environments have been
either unsuccessful or inconclusive in a variety of sources, from protostars to planetary nebulae
\cite[PNe; e.g.,][]{LeBertre_1993_DCBs,Oudmaijer_1997_YSO-DIBs,GarciaHernandez_2013_fullerenedibs}.
Slightly more promising candidates have been found by
\citet{DiazLuis_2015_fullereneDCBs}, who noted absorption features towards the PN Tc~1
at wavelengths consistent with DIBs blue-shifted within the range of the nebular emission velocity;
\citet{Dahlstrom_2013_herschel36dibs} and \citet{Oka_2013_Herschel36DIBs},
interpreted some unusually broad and asymmetric DIBs towards the infrared (IR) source Herschel~36
as most likely due to excitation of the DIB carrier molecules by the local radiation field, but see also \citet{Bernstein_2015_6614DIB}.
Broadening and asymmetry may also be caused by the superposition of DIB features from multiple clouds
along the line of sight,
and many DIBs have been determined to have intrinsic substructure and asymmetry
\citep[e.g., 6203\AA\, and 13175\AA;][]{Porceddu_1991_6200angDIBs,Rawlings_2014_JbandDIBs}.

Here we present the discovery of a
strongly asymmetric instance of
the (normally symmetric) 1.5272~$\mu$m IR DIB \citep[vacuum wavelength;][]{Geballe_2011_IRdibs,Zasowski_2015_dibs}
along the line of sight towards the Red Square Nebula \citep[RSN, also MWC~922;][]{Merrill_1949_mwc}.
The RSN,
a luminous IR source variously classified
as a proto-planetary nebula and a B[e] star \citep{Lamers_1998_Be-classification,Meixner_1999_PPNcandidates},
has strong hydrogen
and metallic emission lines that provide
constraints on the geometrical and kinematical structure of the nebula
\citep{Andrillat_1976_emissionline-IRexcess,Tuthill_2007_redsquare,Whelan_2014_mwc922}.
The circumstellar environment hosts a substantial amount of dust, atomic, and molecular material
\citep{Andrillat_1976_emissionline-IRexcess,vanDiedenhoven_2004_PAHprofiles},
and we demonstrate here that it is also responsible for approximately half of the 1.5272~$\mu$m DIB absorption
along this line of sight, based on evidence from foreground reddening and the absorption feature shape and velocity.

In Section~\ref{sec:data},
we describe the data used
and the extraction of the RSN's 1.5272~$\mu$m DIB profile.
In Sections~\ref{sec:analysis}-\ref{sec:other_interpretations},
we present the properties of this unusual detection, compare it to the typical profiles observed in the ISM,
briefly examine other IR and optical DIBs along this line of sight, and explore plausible interpretations.
In Section~\ref{sec:conclusions},
we summarize our conclusions.

\section{Data} \label{sec:data}

\subsection{$H$-Band Spectra} \label{sec:apogee_spectra}

The primary data used are high-resolution, high signal-to-noise
$H$-band spectra of the RSN observed by the Apache Point Observatory Galactic Evolution Experiment (APOGEE)
survey \citep[][2015 in prep]{Majewski_2012_apogee}.
APOGEE is a high resolution ($R \sim 22,500$), $H$-band ($\lambda = 1.51-1.70~\mu$m)
spectroscopic survey studying the chemodynamical evolution
of the Milky Way (MW) with a large sample of predominantly red giant stars
\citep{Zasowski_2013_apogeetargeting}.
The survey uses a custom-built 300-fiber spectrograph \citep{Wilson_2012_apogee} on the 2.5-m Sloan telescope \citep{Gunn_2006_sloantelescope};
public release of the data began in 2013 with the Sloan Digital Sky Survey III \citep[SDSS-III;][]{Eisenstein_11_sdss3overview}
Data Release 10 \citep[DR10;][]{Ahn_2014_dr10}
and culminated in 2015 with DR12 \citep{Alam_2015_SDSSDR12}.

Because of their usefulness in identifying terrestrial atmospheric (telluric) absorption features,
several thousand early type stars were also observed by the survey.
The RSN was observed as part of this sample on 4, 23, and 25 May 2013.
Figure~\ref{fig:RSN_image} is a Keck NIRC2 $H$-band image
of the RSN, along with the size and position of the APOGEE fiber field of view (red) and some slits (blue and orange)
used for supplementary information on this sightline's optical DIBs.

\begin{figure}[!hptb]
\begin{center}
   \includegraphics[trim=0.1in 4.6in 3in 1.5in, clip, angle=0, width=0.5\textwidth]{f1.pdf}
\end{center}
\caption{
Keck NIRC2 $H$-band image of the RSN (grayscale),
with the apertures of the 2$^{\prime\prime}$ APOGEE fiber (red circle),
a $1.2^{\prime\prime} \times 8^{\prime\prime}$ optical slit (outer blue rectangle),
and a $1.0^{\prime\prime} \times 8^{\prime\prime}$ optical slit (inner orange rectangle).
}
\label{fig:RSN_image}
\end{figure}

All spectra taken with APOGEE are reduced using a custom reduction pipeline \citep{Nidever_2015_apogeereduction},
which produces an ``apVisit'' data file for each stellar observation.
This pipeline calculates the radial velocity (RV) for each visit, but as it is not optimized for emission line spectra,
we calculate our own RVs with a more robust treatment of the spectral features
and use those to stack the apVisit data into a single source spectrum (with ${\rm S/N} > 460$).

\subsection{Extraction of 1.5272~$\mu$m DIB Profile} \label{sec:spectra_fit}

The 1.5272~$\mu$m DIB
overlaps the Brackett-19 hydrogen emission line (Br19; $\lambda_0 = 15\,264.725$~\AA).
Because the DIB appears very close to where the Br19 line drops to the continuum,
careful treatment of the intrinsic Br19 profile shape is required.
Here, we model the Br19 profile empirically using all of the unblended
hydrogen lines available in the APOGEE spectrum (i.e., Br11, Br13, Br16, Br17, and Br20).
In the rest velocity frame, we calculate a quadratic relationship between the observed flux
$F_\lambda(\Delta v)$ and the $\lambda_0$ at each $\Delta v$ of the Brackett lines.
Then, we interpolate this relationship at the Br19 rest wavelength to construct our empirical Br19 profile.
Figure~\ref{fig:RSN_fits}a shows these profiles (colored lines) and the interpolated Br19 model (black dashed line),
along with the observed Br19+DIB profile (thick blue line).
We measure a heliocentric RV of 37.8 km~s$^{-1}$ from this line,
corresponding to a velocity relative to the Local Standard of Rest ($v_{\rm LSR}$) of 49.8~km~s$^{-1}$.

\begin{figure*}[!hptb]
\begin{center}
  \includegraphics[trim=0.5in 5.8in 0.9in 1.5in, clip, angle=0, width=\textwidth]{f2.pdf}
\end{center}
\caption{
Measurement of the overlappping Br19 and 1.5272~$\mu$m DIB profiles.
{\it (a):} The unblended Brackett profiles (thin colored lines) are used to interpolate the intrinsic Br19 profile (dashed black line).
The observed Br19+DIB profile is shown in blue.
{\it (b):} The residual DIB absorption, versus velocity $v_{\rm LSR}$,
with single (top, red) and dual (bottom, green) Gaussian fits.
The vertical lines indicate the centers of the Gaussian components, and the gray arrow is the $v_{\rm LSR}$ of the RSN's Br19 line.
}
\label{fig:RSN_fits}
\end{figure*}

\section{Analysis} \label{sec:analysis}

Figure~\ref{fig:RSN_fits}b shows the strong residual absorption ($R_\lambda$) remaining after normalization
by this empirical Br19 profile model
(with two fitting scenarios that are discussed in Section~\ref{sec:dib_fits}).
To our knowledge, there are no other atomic or molecular transitions at this wavelength that could
be strong and broad enough in this source to produce the absorption observed.

We calculate an integrated equivalent width $W_{\rm DIB} = 557$~m\AA\,
by summing $(1-R_\lambda)$ in the window $-150 \le v_{\rm LSR} \le 250$~km~s$^{-1}$;
random uncertainties based on the ``flat'' span of $R_\lambda$ between $-300 \le v_{\rm LSR} \le -150$~km~s$^{-1}$ are $\lesssim$1~m\AA.
Assuming the mean ratio $W_{\rm DIB}/A_V=102$~m\AA~mag$^{-1}$ for this DIB \citep{Zasowski_2015_dibs},
and a typical $R_V=3.1$ extinction law, this equivalent width corresponds to a line-of-sight dust reddening of $E(B-V) \approx 1.75$.
At high reddening, the dispersion in the $W_{\rm DIB}/A_V$ ratio is $\lesssim$40~m\AA~mag$^{-1}$.

Of course, for any circumstellar material, the extinction law (and $W_{\rm DIB}/A_V$ ratio) may be atypical.
However, we can make some independent estimates of the interstellar reddening.
The optical DIBs also measured as part of this analysis
suggest a total reddening of $E(B-V) \sim 1.0-1.8$~mag, with a dispersion of about 0.5~mag,
depending on the optical spectrum used (see Section~\ref{sec:optical_dibs}).
The foreground-only reddening, as measured from PanSTARRS+2MASS photometry by
\citet{Green_2014_PS1dustmap,Green_2015_PS1dustmap_aas}
for the $7^\prime \times 7^\prime$ region around the RSN's position,
is $E(B-V) \approx 0.57$$^{+0.03}_{-0.08}$ and $\approx 0.87$$^{+0.05}_{-0.05}$ mag
at heliocentric distances of $\sim$2 and $\sim$3~kpc, respectively.
These lower values are consistent with the presence of {\it additional} dust along the line of sight being traced by the IR DIB,
and perhaps some of the optical ones.

\subsection{Asymmetry of Profile}

To quantitatively compare the asymmetry in the RSN's 1.5272~$\mu$m DIB
to that of the typical 1.5272~$\mu$m DIB features observed in the ISM,
we use a simplified version of the line bisector estimator
\cite[e.g.,][]{Toner_1988_starlinemeasurement},
which compares the midpoints of a profile near its top and bottom (15\% and 85\% of the amplitude; left panel of Figure~\ref{fig:rsn_bis}).
A difference in these midpoints (here, $\delta v$) indicates asymmetry in the profile.
When compared using this metric to a set of several thousand 1.5272~$\mu$m DIB features
that sample the ISM towards cool giant stars\footnote{This comparison
set comprises all sightlines from \citet{Zasowski_2015_dibs} with $\sigma(R_\lambda) \le 0.03$,
$N_{\rm visits} = 1$ or $RV_{\rm scatter} \le 0.5$~km~s$^{-1}$,
$\sigma(R_\lambda)/\sigma(F_\lambda) < 0.55$, and $N_{\rm pix} \ge 1350$; see Section~2.2 of that paper for details.}
at identical spectral resolution,
the RSN sightline appears more asymmetric than all but $\sim$2\% of the larger sample (right panel of Figure~\ref{fig:rsn_bis}).
Note that this comparison set does contain noisy sources, because restricting it to only those sightlines with ``good''
DIB fits would explicitly remove any asymmetric features.  Visual inspection of those $\sim$2\% asymmetric sources
confirms that the vast majority are contaminated by residual stellar features,
though we cannot rule out a few potentially intrinsically
lopsided DIBs that will be good candidates for follow-up observations.
None of the sightlines with potentially asymmetric 1.5272~$\mu$m DIB features lie closer than 2.4$^\circ$ to the RSN,
and none of the DIBs in the $\sim$60 sightlines in the dataset at closer projected distances
show any signs of asymmetry.

\begin{figure*}[!hptb]
\begin{center}
   \includegraphics[trim=5.5in 1.2in 1in 4.4in, clip, angle=180, width=0.45\textwidth]{f3a.pdf}
   \includegraphics[trim=0.75in 4.5in 5.7in 1.1in, clip, angle=180, width=0.45\textwidth]{f3b.pdf}
\end{center}
\caption{
{\it Left:} The measurement of the asymmetry estimate $\delta v$ for the RSN's DIB.
{\it Right:} The distribution of ISM $\delta v$ values, compared to that from the RSN sightline.
}
\label{fig:rsn_bis}
\end{figure*}

\subsection{DIB Fit Components} \label{sec:dib_fits}

We turn now to different options for fitting this feature.
The upper $R_\lambda$ profile in Figure~\ref{fig:RSN_fits}b is overplotted in red with the
Gaussian that best reproduces the blue (lower velocity) side of the feature;
the extended red excess is obvious.  (Fitting the entire feature with a single Gaussian yields significantly noisier residuals.)
The bulk of the feature is well-represented by this Gaussian with a FWHM of 3.5~\AA.
If this DIB feature has a single origin, the red excess may be due to changes in the DIB carrier's electronic properties (Section~\ref{sec:other_interpretations}).
The fact that the deepest part of the feature appears at {$v_{\rm LSR} \sim 12$~km~s$^{-1}$,
nearly 37~km~s$^{-1}$ offset from the nebular Brackett emission, is a strike against its origin in the nebula.
The mismatch does not completely rule out a circumstellar origin, as a
range of velocities (of tens of km~s$^{-1}$) is observed in other
absorption and emission lines definitively associated with the nebula,
arising from
the superposition of multiple nebular components
\citep[e.g.,][]{Whelan_2014_mwc922}.
However, a simpler solution may be sought.

The preferred scenario is shown with the lower $R_\lambda$ profile of Figure~\ref{fig:RSN_fits}b,
overplotted in green with the best simultaneous fit dual-Gaussian components,
which together reproduce the entire feature very well.
Each component corresponds to an approximately equal $W_{\rm DIB}$ contribution of 277~m\AA,
and thus an $E(B-V) \approx 0.87$~mag.  This is consistent with the {\it foreground} reddening predicted by the
PanSTARRS+2MASS map in the distance range of the RSN, implying an additional reddening of 0.87 mag due to {\it circumstellar} material.
Using instead a dense environment extinction curve of $R_V = 5$ with one component,
this additional reddening is closer to $E(B-V) \sim 0.54$~mag,
which brings the total reddening to 1.4~mag.
Both of these are within the (albeit broad) range of reddening spanned by the optical DIBs,
suggesting that some of the optical DIB absorption may arise from the circumstellar medium as well.

The use of reddening arguments {\it alone} as evidence for circumstellar DIBs is notoriously tricky,
because of the large uncertainties and intrinsic dispersions in the DIB/$E(B-V)$ relations \citep[e.g.,][]{Luna_2008_DIBs-in-pAGBs}.  Thus,
further critical evidence for this scenario lies in the line-of-sight velocities of the two components.
The narrower feature (${\rm FWHM} = 3.1$~\AA) has $v_{\rm LSR} = 10.3$~km~s$^{-1}$,
well within the range spanned by the optical DIBs and the interstellar Na~D1 and K~I lines.
In contrast, the wider feature has $v_{\rm LSR} = 49.6$~km~s$^{-1}$,
nearly identical to that measured for the Br emission lines of the RSN itself.
Under the interpretation of a circumstellar origin, the larger width (${\rm FWHM} = 8.7$~\AA)
may be seen as simply velocity broadening due to turbulence or rotation within the nebula.
If the foreground component approximately represents the intrinsic feature width,
this implies $\sim$47~km~s$^{-1}$ of broadening within the nebula, consistent with the range of velocities
observed in different nebular components \citep[e.g.,][]{Whelan_2014_mwc922}.
Thus, if the two postulated components indeed represent different physical origins,
then they serve as strong evidence for DIB carriers in the {\it circumstellar} environment of the RSN.

\subsection{Other $H$-Band DIBs} \label{sec:hband_dibs}

Ten known DIBs lie within APOGEE's IR wavelength coverage \citep{Geballe_2011_IRdibs,Cox_2014_xshooterDIBs},
including the 1.5272~$\mu$m feature.
Based on the small number of sightlines in which the others have been observed,
their strengths are expected to be at most $\sim$30\% of the 1.5272~$\mu$m DIB, if detected at all.
Despite the strong 1.5272~$\mu$m absorption in the RSN, however,
the difficulty of measuring these weaker bands is greatly exacerbated
(indeed, dominated) by the complex spectrum of the RSN ---
the numerous strong, blended, and sometimes unidentified emission
lines make determination of the underlying continuum extraordinarily tricky.
(The 1.5272~$\mu$m DIB lies in an exceptionally ``smooth'' region.)
Figure~\ref{fig:other_dibs} demonstrates this with six of the DIBs.
These features are indeed quite weak, but we do not see strong evidence for any feature asymmetries.

\begin{figure*}[!hptb]
\begin{center}
   \includegraphics[trim=0.7in 3in 2in 1.8in, clip, width=0.45\textwidth]{f4a.pdf}
   \includegraphics[trim=0.7in 3in 2in 1.8in, clip, width=0.45\textwidth]{f4b.pdf}
\end{center}
\caption{
Examples of other DIBs in APOGEE's wavelength range.
{\it Top:} Observed flux (black) and polynomial continuum (green).
{\it Bottom:} Continuum normalized flux, with DIB central wavelengths and FWHMs
\citep[black brackets; average FWHMs from][]{Geballe_2011_IRdibs,Cox_2014_xshooterDIBs}.
Pixels flagged as BAD or with high errors in any of the apVisit flux arrays have been masked.
}
\label{fig:other_dibs}
\end{figure*}

\subsection{Optical DIBs} \label{sec:optical_dibs}

We also studied two high-resolution optical spectra, to assess the properties of optical DIBs and interstellar lines
in the context of our proposed scenario.
One spectrum (hereafter, ``MDO'') was observed with the Harlan~J.~Smith 2.7~m telescope at McDonald Observatory
on 6/7 September 2009, using the
$1.2^{\prime\prime}$ slit shown in blue in Figure~\ref{fig:RSN_image}.
The $R=60,000$ spectrum, which spans 4787--8861~\AA\, and
has been smoothed with a moving boxcar filter (with a kernel of 11 points),
has a S/N of $\sim$20 near 4860~\AA\, and $\sim$90 near 6563~\AA.
The other (hereafter, ``UVES'') was observed with the ESO VLT+UVES on 27 August 2009, using the
$1^{\prime\prime}$ slit shown in orange in Figure~\ref{fig:RSN_image}.
The UVES spectrum spans 5655--9464\AA\, with $R=42,000$ and a median S/N of 47.2.
Of the DIBs that are in regions relatively free of contaminating lines,
we could reliably fit seven features using a simple Gaussian, shown in Figure~\ref{fig:optical_dibs} (MDO on the left, UVES on the right).

In comparing the two spectra, the largest differences are in the NaD doublet, which is saturated in the UVES
spectrum but not in the MDO one, and in the strength of the 5780\AA\, and 5797\AA\, DIB features (bottom curves in
Figure~\ref{fig:optical_dibs}).  The UVES 5780\AA\, and 5797\AA\, DIB equivalent widths are approximately twice as high
as their counterparts in the MDO spectrum, and the other DIB equivalent widths are between 1.3--1.9$\times$ as high.
The slit widths used differ only slightly (1.2$^{\prime\prime}$ for MDO and 1.0$^{\prime\prime}$ for UVES),
and the slits were centered on the same position to within the uncertainty of the telescope pointing reported in the file headers.
To assess the possible impact of small angular offsets sampling small-scale variations in the ISM across the nebula,
we also visually examined the NaD and K~I lines from three optical VLT+X-Shooter spectra\footnote{These data, spanning 534--1020\AA,
were obtained from the ESO archive.  One spectrum were observed on 02 July 2010 with median ${\rm S/N}=123$, and the
other two on 07 July 2012, with median ${\rm S/N}=145$ and 144.} separated from each other by $1-2^{\prime\prime}$ on the sky.
The spectral resolution is lower, but we see no difference in the shape or strength of these interstellar features.
More detailed mapping of any ISM structure or DIB variability across the face of the nebula is reserved for future work;
here, we take the consistency of the X-Shooter spectra as evidence that the discrepancies between the MDO and UVES spectra
are not due to strong variations in the ISM density on the scale of arcseconds.
For the remainder of this analysis, we will quote values derived from the two spectra separately.

Qualitatively, the spectra are very similar.
The critical point is that none of the DIB features, in either optical spectrum,
display the strong asymmetry seen in the 1.5272~$\mu$m DIB
\citep[except possibly 6203\AA, which has been previously noted as intrinsically asymmetric;][]{Porceddu_1991_6200angDIBs}.
We attempted to fit two simultaneous Gaussians to these features,
emulating the dual components seen in the 1.5272~$\mu$m DIB,
but the resulting fits were {\it quantitatively} no better than those with a single Gaussian (including for 6203\AA).
We note that at these wavelengths, the shift expected for a 37.8~km~s$^{-1}$ offset is $\sim$0.7\AA,
which is discernible even in these broad (and potentially substructured) features;
however, there are degeneracies among the fit parameters (e.g., amplitude ratio vs.\, velocity broadening),
and thus we cannot conclusively rule out the presence of red excess absorptions
with strength ratios or velocity properties slightly different from those of the IR DIB components.

\begin{figure*}[!hptb]
\begin{center}
   \includegraphics[trim=1in 4.9in 1.1in 1.8in, clip, width=\textwidth]{f5.pdf}
\end{center}
\caption{
Single Gaussian fits (blue lines) to seven optical DIBs.  The normalized spectra are offset by increments of 0.5 for visual clarity.
{\it Left:} Smoothed $R=60,000$ spectrum from McDonald Observatory.
{\it Right:} Original (gray) and smoothed (black) $R=42,000$ spectrum from VLT+UVES.
}
\label{fig:optical_dibs}
\end{figure*}

The total equivalent widths of the optical features ($W_{\rm opt}$), combined with the
$W_{\rm opt}/W_{1.5272\mu m}$ relations calculated by \citet{Cox_2014_xshooterDIBs},
predict an absorption of the 1.5272~$\mu$m DIB that is $\sim$$30-50\pm10$\% what is observed
(from the MDO and UVES data, respectively);
we note that these relations are based on a small number of sightlines,
which probe interstellar material with potentially different chemistry
or radiation exposure than that towards or in the RSN.

Using the mean $W_{\rm opt}/E(B-V)$ relations from \citet{Jenniskens_1994_dibs}, we calculate
a reddening of $E(B-V) \sim 1.0\pm0.2$~mag (MDO) and $E(B-V) \sim 1.8\pm0.5$~mag (UVES).
These are larger than the foreground reddening from the PanSTARRS+2MASS map at 2--3~kpc,
approaching that of the total estimated from the 1.5272~$\mu$m DIB.
This is consistent with a circumstellar contribution to the DIB absorption, on the order
of 15--50\% of the total.
This result on its own is merely suggestive, and we emphasize the large dispersion in the $W_{\rm opt}/E(B-V)$ relations
(both intrinsic scatter and that due to measurement uncertainties), as noted in Section~\ref{sec:dib_fits}.
However, when combined with the velocity comparisons below,
it is worth noting the consistency of the reddening with excess DIB absorption as measured in several DIBs.

Adopting the central wavelength $\lambda_c$ values of \citet{Hobbs_2008_DIBcatalog},
we find DIB $v_{\rm LSR}$ values spanning $2-20$~km~s$^{-1}$,
centered on the value measured for the ``foreground'' 1.5272~$\mu$m DIB ($v_{\rm LSR} = 10.3$~km~s$^{-1}$),
the interstellar 5895~\AA\, Na~D1 line in the MDO spectrum ($v_{\rm LSR} = 8.5$~km~s$^{-1}$),
and the interstellar 7701~\AA\, K~I line in the UVES spectrum ($<$$v_{\rm LSR}$$> = 15.2$~km~s$^{-1}$).
The optical DIB feature widths are also consistent with the widths from the interstellar cloud studies of
\citet{Jenniskens_1994_dibs} and \citet{Hobbs_2008_DIBcatalog}, with some scatter of $\lesssim$8\%;
as noted above, though a dual-component fit does not produce a quantitatively better match to these features,
we cannot rule out the possibility that a small fraction of circumstellar material may be contributing to the observed widths.

\section{Other Interpretations} \label{sec:other_interpretations}

There exist other potential explanations for the ``additional'' long-wavelength 1.5272~$\mu$m DIB absorption,
but all face serious problems.
One interpretation is that the ISM's velocity distribution towards the RSN has a peculiar structure,
such that the ``high'' $v_{\rm LSR}$ feature is due to an unassociated ISM cloud.
In a two-cloud scenario, this would require the second one to have nearly precisely the same radial velocity as the RSN
and an internal velocity dispersion of $\gtrsim$47~km~s$^{-1}$,
neither of which is readily apparent in the observed profiles of the optical DIBs.

To assess this in greater detail, we examined the velocity structure of the 7701~\AA\, K~I line from the UVES spectrum (Section~\ref{sec:optical_dibs}).
This is shown in the top panel of Figure~\ref{fig:ISM_velocity_model}, along with the best three-Gaussian fit
representing the three resolved subpeaks.
In the bottom panel, we ``model'' the 1.5272~$\mu$m DIB with these three components,
fixing the central velocities, velocity dispersions, and relative amplitudes, but scaling the
shortest-wavelength amplitude peak to match the deepest DIB absorption.
These components (and their summed total) are shown in red in the bottom panel of Figure~\ref{fig:ISM_velocity_model},
compared with the dual-component best fit presented in Figure~\ref{fig:RSN_fits}b and Section~\ref{sec:dib_fits}.
It is clear that while the K~I profile does reflect ISM components at higher line-of-sight velocities,
these are not sufficient to explain the long-wavelength 1.5272~$\mu$m DIB absorption,
under the assumption that the DIB carrier and K~I are mixed in the ISM \citep[like carriers of other DIBs are; e.g,][]{Galazutdinov_2005_ISMlines}.

\begin{figure}[!hptb]
\begin{center}
   \includegraphics[trim=1.0in 5.0in 4.1in 1.8in, clip, width=0.5\textwidth]{f6.pdf}
\end{center}
\caption{
Extrapolating a multi-component 1.5272~$\mu$m DIB from the K~I ISM line.
{\it Top:} The UVES spectrum towards the RSN, showing the clearly structured K~I line at 7701\AA.
The blue lines represent the best-fitting Gaussian components (dotted) and summed total (solid).
{\it Bottom:} The observed 1.5272~$\mu$m DIB (black),
with the best-fit components presented in Figure~\ref{fig:RSN_fits}b and Section~\ref{sec:dib_fits} (green).
The red dashed lines are simulated DIB components with the same central velocities, velocity dispersions,
and relative amplitudes as the K~I components above.  The solid red line is their sum.
}
\label{fig:ISM_velocity_model}
\end{figure}

Another possible explanation is the presence of an additional,
hitherto undetected DIB blended with the 1.5272~$\mu$m feature
\citep[e.g., as postulated by][for the Her~36 asymmetric DIBs]{Bernstein_2015_6614DIB},
or intrinsic asymmetry in this feature.
However, this is not supported by any other observation of this feature to date,
including in heavily reddened, early-type stars with strong 1.5272~$\mu$m features that are symmetric
at a high level of confidence.
Some examples are shown in the right-hand panels of Figure~\ref{fig:comps_interp} (\ref{fig:comps_interp}b and \ref{fig:comps_interp}d).
These spectra are also APOGEE data and have a resolution identical to that shown in Figure~\ref{fig:RSN_fits}.
The red lines are single Gaussian fits.
If the long wavelength extension of the RSN's 1.5272~$\mu$m DIB were a result of an intrinsic asymmetry or a distinct DIB feature,
it should be apparent in these spectra as well, but all sightlines show a clear symmetry.

We confirmed that the asymmetry is not an artifact intrinsic to the APOGEE data by comparing the Br19+DIB profile
to that in an X-Shooter spectrum\footnote{Obtained from the ESO archive, spanning 994--2479\AA,
observed on 02 July 2010 with a median ${\rm S/N}=538$.}; we find that they share the same shape.
Figure~\ref{fig:comps_interp} in full presents two additional examples of extracting 1.5272~$\mu$m DIB features
from the APOGEE spectra of stars with Br19 emission, using the methodology described in Section~\ref{sec:spectra_fit}.
In the same style as Figure~\ref{fig:RSN_fits}a,
Figure~\ref{fig:comps_interp}a shows the weak, double-peaked Br emission profiles of ABE-022 \citep[][]{Chojnowski_2015_BeStars},
with the DIB clearly visible.  Figure~\ref{fig:comps_interp}c shows the same for ABE-137, a star with single-peaked emission
profiles more like the RSN.  The right-hand panel of each row contains the residual absorption after removal of the interpolated
Br19 profile as black lines, along with a single Gaussian fit (red lines).
For comparison, we also reproduce the RSN's 1.5272~$\mu$m DIB in Figure~\ref{fig:comps_interp}b in blue, and
in Figure~\ref{fig:comps_interp}d, we overplot in green the normalized spectrum of ABE-Q09 \citep[also][]{Chojnowski_2015_BeStars},
a very early type star with no emission lines to complicate the DIB measurement.
The isolated DIB absorptions are highly symmetric, with shapes that are nearly perfectly Gaussian and that are extremely similar to the
DIB profile in the emission-free spectrum.
This demonstrates clearly that the asymmetry in the RSN's 1.5272~$\mu$m DIB is not an artifact of our deblending process.

\begin{figure}
\begin{center}
   \includegraphics[trim=0.5in 2.9in 0.3in 1.4in, clip, width=0.5\textwidth]{f7.pdf}
\end{center}
\caption{
Extraction of the 1.5272~$\mu$m DIB feature from two additional emission-line stars \citep[ABE-022 and ABE-137;][]{Chojnowski_2015_BeStars}.
Panels~({\it a}) and ({\it c}) are identical in form to Figure~\ref{fig:RSN_fits}a.
Panels~({\it b}) and ({\it d}), like Figure~\ref{fig:RSN_fits}b, show the residual absorption in these spectra after the Br19 feature has been removed (black lines),
along with a single Gaussian fit (red lines).
In panel~({\it b}), we also show the RSN's DIB profile (blue line) for comparison.  In panel~({\it d}),
we scale and overplot the spectrum of the early type star ABE-Q09 \citep[green line][]{Chojnowski_2015_BeStars}.
See text for details.
}
\label{fig:comps_interp}
\end{figure}

We cannot fully rule out the possibility that the RSN's 1.5272~$\mu$m DIB feature arises entirely from circumstellar carriers strongly impacted
by their local environment.
\citet{Dahlstrom_2013_herschel36dibs} reported the discovery of optical DIBs with pronounced red wings along sightlines towards Herschel~36,
a system of O stars illuminating the M8 HII region and in close proximity to the Her~36~SE strong IR source.
Qualitatively, these wings are very similar to the one observed here in the RSN (see, e.g., their Figure 6).
Chemical models calculated by \citet{Oka_2013_Herschel36DIBs} demonstrate that these wings can be produced
when a strong radiation field drives some electrons in the DIB carrier molecules into excited states.
However, the need to invoke an outflow of nearly exactly the same velocity as the RSN's motion relative to the Sun,
such that the main feature is centered at the same velocity as the foreground ISM,
renders this a less likely interpretation.

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

We have presented the detection of a peculiar, highly asymmetric instance of the
1.5272~$\mu$m DIB along the line of sight towards the Red Square Nebula,
using high resolution $H$-band spectra from the APOGEE survey.
Strongly asymmetric profiles have never yet been observed in this particular DIB feature, or any at longer wavelengths.

The likeliest explanation is that the feature towards the RSN is a blend of absorption arising from the foreground ISM
and from the nebular material itself.  This interpretation is supported by the radial velocities of the two blended components,
which align nearly perfectly with the foreground ISM (e.g., as traced by interstellar Na or K)
and with the RSN's hydrogen emission lines, respectively.
In addition, the amount of interstellar reddening ``predicted'' by the interstellar DIB component nearly exactly matches
the {\it foreground} reddening traced by high-quality three-dimensional dust maps.
The {\it total} reddening predicted by the whole IR DIB feature (interstellar plus circumstellar) is
roughly comparable than the reddening predicted by the optical DIBs,
suggesting that at least some of the optical DIB carriers reside in the nebula itself
or that the relationship between the 1.5272~$\mu$m DIB and the optical DIBs, currently only coarsely constrained,
is substantially different in the circumstellar material.

The asymmetric 1.5272~$\mu$m absorption feature observed towards the Red Square Nebula
is likely the most confident detection to date of a diffuse circumstellar band
\citep[DCB; see also][]{DiazLuis_2015_fullereneDCBs},
and if genuine, is also the first detection of an IR DCB,
cleanly separated from the ISM along the same line of sight in both velocity and dust reddening.
Because the circumstellar environment (its radiative properties, temperature profile, geometry, etc.) can be understood and disentangled,
this source will become a powerful laboratory for constraining the response of the DIB carriers to these environmental factors,
and thus constraining the nature of the carriers themselves.

\begin{acknowledgments}

We thank the referee for helpful comments that improved the clarity and strength of the paper.

GZ is supported by an NSF Astronomy \& Astrophysics Postdoctoral Fellowship under Award No.\ AST-1203017.

Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is \url{http://www.sdss3.org/}.

SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, University of Cambridge, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, The Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

This paper includes data taken at The McDonald Observatory of The University of Texas at Austin, and
data obtained from the ESO Science Archive Facility under programme IDs 083.C-0676(A), 385.C-0720(A), and 089.C-0874(A).

\end{acknowledgments}

\bibliographystyle{aa}


\section{Acknowledgements}
We thank Luming Duan, Zhe-Xuan Gong, Tarun Grover, Chris Laumann, Shengtao Wang, Norman Yao, and Jiehang Zhang for helpful discussions. This work is supported by the ARO Atomic and Molecular Physics Program, the AFOSR MURI on Quantum Measurement and Verification, the IARPA MQCO program, and the NSF Physics Frontier Center at JQI. PH and MH acknowledge  the Deutsche Akademie der Naturforscher Leopoldina (grant No.~LPDS 2013-07), the EU IP SIQS, the SFB FoQuS (FWF Project No.~F4016-N23) and the ERC synergy grant UQUAM.