# Tulipa Energy Model: Mathematical Formulation

Diego A. Tejada-Arangoa, Germán Morales-Españaa, Lauren Clisbya, Ni Wanga,
Abel S. Siqueirab, Ali Subayu b, Laurent Soucasse b, and Zhi Gao c

a TNO - Energy Transition Studies
b eScience Center
c Utrecht University

## Abstract

Tulipa aims to optimise the investment and operation of the electricity
market, considering its coupling with other sectors, such as hydrogen and
heat, that can also be electrified. The problem is analysed from the
perspective of a central planner who determines the expansion plan that is
most beneficial for the system as a whole, either by maximising social welfare
or by minimising total costs. The formulation provides a general description
of the objective function and constraints in the optimisation model based on
the concept of energy assets representing any element in the model. The model
uses subsets and specific methods to determine the constraints that apply to a
particular technology or network, allowing more flexibility in the code to
consider new technologies and constraints with different levels of detail in
the future.

##  1 Tulipa Energy Model: An electricity market and sector-coupling model for
investment and operation

In general terms, this document presents a classical problem of joint
generation and transmission expansion planning in an electric energy system
[1]. However, Tulipa’s formulation is more generalised and includes other
sectors coupled with the electricity sector. In addition, the formulation is
based on Graph Theory [2], which provides a more flexible framework to model
energy assets in the system as vertices and flows between energy assets as
edges. Figure (1) shows an example of an energy system network coupling
electricity, methane (i.e., gas), and hydrogen ($H_{2}$).

Figure 1: Test system in Tulipa Energy Model

The following sections describe the notation, constraints, and objective
function in the formulation used in the model.

##  2 Notation

###  2.1 Sets for assets

The formulation in Tulipa relies on the definitions of energy assets (i.e.,
vertices of the network), which have been taken from the Energy System
Description Language (ESDL) [3].

Name |  Description | Elements
---|---|---
$\mathcal{A}$ |  Energy assets | $a\in\mathcal{A}$
$\mathcal{A}_{cv}$ |  Conversion assets | $\mathcal{A}_{cv}\subset\mathcal{A}$
$\mathcal{A}_{p}$ |  Production assets | $\mathcal{A}_{p}\subset\mathcal{A}$
$\mathcal{A}_{c}$ |  Consumption assets | $\mathcal{A}_{c}\subset\mathcal{A}$
$\mathcal{A}_{t}$ |  Transport assets | $\mathcal{A}_{t}\subset\mathcal{A}$
$\mathcal{A}_{s}$ |  Storage assets | $\mathcal{A}_{s}\subset\mathcal{A}$

As shown in the previous Table, energy assets, denoted by $a$, can belong to
different subsets based on the definitions, such as conversion, production,
consumption, transport, and storage assets. Here there are some examples:

  * •

Conversion ($\mathcal{A}_{cv}$): Power plants, fuel cells, heat pumps, etc.

  * •

Production ($\mathcal{A}_{p}$): Solar PV panels, wind farms, etc.

  * •

Consumption ($\mathcal{A}_{c}$): Electricity demand, heat consumption of a
city, etc.

  * •

Transport ($\mathcal{A}_{t}$): Electricity nodes and lines, gas networks,
district heating, etc.

  * •

Storage ($\mathcal{A}_{s}$): Battery, pumped storage, heat buffer, hydrogen
storage, etc.

Figure 2: Energy assets in ESDL

General constraints can be created at the energy asset level ($a$), such as
maximum capacity, and specific constraints can be applied to these subsets or
groups, such as unit commitment constraints for the conversion units
($\mathcal{A}_{cv}$). This structure allows the model to be flexible in
representation and levels of detail.

###  2.2 Extra sets

The formulation also uses more set definitions, such as:  Name Description
Elements $\mathcal{A}_{b}$ Assets with balance constraint method
$\mathcal{A}_{b}\subseteq\mathcal{A}$ $\mathcal{A}^{in}(a)$ Assets that are
inputs for asset $a$ $\mathcal{A}^{in}(a)\subseteq\mathcal{A}$
$\mathcal{A}^{out}(a)$ Assets that are outputs for asset $a$
$\mathcal{A}^{out}(a)\subseteq\mathcal{A}$ $\mathcal{Y}$ Years
$y\in\mathcal{Y}$ $\mathcal{Y}_{m}$ Milestone years
$\mathcal{Y}_{m}\subseteq\mathcal{Y}$ $\mathcal{K}$ Representative periods
$k\in\mathcal{K}$ $\mathcal{T}$ Time steps for the lowest time resolution in
the system $t\in\mathcal{T}$ $\tau(a)$ Time steps defined for each asset $a$
$\tau(a)\subseteq\mathcal{T}$ $\mathcal{P}$ Auxiliary set for constraints with
several time steps definitions $p\in\mathcal{P}\subseteq\mathcal{T}$
$\mathcal{R}$ Types of system reserves $r\in\mathcal{R}$ $\mathcal{R}_{a}^{+}$
Upwards reserves provided by asset $a$
$\mathcal{R}_{a}^{+}\subseteq\mathcal{R}$ $\mathcal{R}_{a}^{-}$ Downwards
reserves provided by asset $a$ $\mathcal{R}_{a}^{-}\subseteq\mathcal{R}$
$\Omega_{ap\tau}$ Set to map between two temporal definitions $p$ and $\tau$

The model uses sets to define methods, which are used to determine constraints
and detail levels for energy assets. For instance, depending on whether it
belongs to $\mathcal{R}_{a}^{+}$ or $\mathcal{R}_{a}^{-}$, an energy asset may
have a reserve requirement method. The balance method determines if the asset
uses the balance constraint when it belongs to $\mathcal{A}_{b}$. The concept
of methods has been previously defined in [4] and [5].

###  2.3 Parameters

The model’s parameters are listed in the table below. Unit measures and scales
may vary based on input data.

Name |  Description | Units
---|---|---
$\overline{F}_{a}$ |  Maximum asset capacity | [MW]
$\overline{F}^{max}_{ya}$ |  Maximum potential capacity | [MW]
$\tilde{F}^{max}_{yakt}$ |  Maximum energy asset production profile | [p.u.]
$\tilde{F}^{min}_{yakt}$ |  Minimum energy asset production profile | [p.u.]
$\overline{U}_{ya}$ |  Initial installed assets | [-]
$S_{yakt}^{I}$ |  Energy asset inflows per representative | [MWh]
$\eta_{a\alpha}$ |  Efficiency when converting energy from $a$ to $\alpha$ | [p.u.]
$R_{yrkt}$ |  Reserve requirement | [MW]
$D_{yakt}$ |  Demand of asset $a$ | [MWh]
$P_{yakt}$ |  Production of asset $a$ | [MWh]
$LT_{a}$ |  Lifetime of the asset $a$ | [years]
$W_{y}^{m}$ |  Weight of milestone year $y$ | [years]
$W_{yk}^{op}$ |  Weight of representative period $k$ for operational cost | [hours]
$M_{\alpha p\tau}^{op}$ |  Mapping matrix between two temporal definitions $p$ and $\tau$ | [-]
$C^{I}$ |  Total investment cost of the energy system | [€ ]
$C^{O}$ |  Total operational cost of the energy system | [€ ]
$C_{ya}^{op}$ |  Variable operational cost of asset $a$ at year $y$ | [€ /MWh]
$C_{ya}^{T}$ |  Total investment cost of the energy asset $a$ at year $y$ | [€ /MW]
$SV_{ya}$ |  Salvage cost of the energy asset $a$ at year $y$ | [€ /MW]
$IR$ |  Interest rate | [p.u.]

###  2.4 Continuous Variables

Name |  Description | Units
---|---|---
$f_{ya\alpha kt}$ |  Flow from asset $a$ to asset $\alpha$ | [MWh]
$r_{yarkt}$ |  Reserve provision of asset $a$ to reserve requirement $r$ | [MW]
$s_{yakt}$ |  Storage level of asset $a$ | [MWh]

###  2.5 Binary and Integer Variables

The following variables are, by definition, integers; however, these
definitions could be relaxed to a continuous variable for simplicity during
the solution process.  Name Description Units $u_{yakt}$ Number of units on
[-] $\overline{u}_{ya}$ Invested number of units [-]

##  3 Constraints and basic methods for energy assets

This section includes the constraints and basic methods for all energy assets
$a$, i.e., producers, storage units, conversion units, and transport elements.

###  3.1 Balance constraint for assets (One constraint to rule them all!)

This constraint is the foundation of the modelling framework in Tulipa.
Equation (1) provides a general form, but simpler versions of this constraint
exist depending on the asset method. It also includes the possibility of
having flexible temporal resolution on each asset [6], allowing the model to
have sectors with different temporal resolutions (e.g., hourly decisions for
the electricity sector and six-hour decisions for the gas sector). In the
following sections, we will explain these versions with examples to enhance
comprehension.

$\displaystyle\underbrace{\sum_{\alpha\in\mathcal{A}^{in}(a)}\quad\sum_{\tau(\alpha)\in\Omega_{\alpha
p\tau}}M_{\alpha p\tau}f_{y\alpha ak\tau}}_{\textit{flows into the asset}}$
(1)
$\displaystyle+\underbrace{\sum_{\tau(a)\in\Omega_{ap\tau}}M_{ap\tau}P_{yak\tau}}_{\textit{for
producer assets}}$
$\displaystyle+\underbrace{\sum_{\tau(a)\in\Omega_{ap\tau}}M_{ap\tau}\cdot\left(s_{yak\tau}-s_{yak,\tau-1}+S_{yak\tau}^{I}\right)}_{\textit{for
storage assets}}$ $\displaystyle\left\\{\begin{array}[]{l}=\\\ \geqslant\\\
\leqslant\end{array}\right\\}$
$\displaystyle\underbrace{\sum_{\alpha\in\mathcal{A}^{out}(a)}\quad\sum_{\tau(\alpha)\in\Omega_{\alpha
p\tau}}M_{\alpha p\tau}\frac{f_{ya\alpha
k\tau}}{\eta_{a\alpha}}}_{\textit{flows outgoing the asset}}$
$\displaystyle+\underbrace{\sum_{\tau(a)\in\Omega_{ap\tau}}M_{ap\tau}D_{yak\tau}}_{\textit{for
consumption assets}}$
$\displaystyle\qquad\forall{y}\in\mathcal{Y}_{m},\forall{a}\in\mathcal{A},\forall{k}\in\mathcal{K},\forall{p}\in\mathcal{P}$

In equation (1), the left-hand side (LHS) shows all the flows entering asset
$a$ from other connected assets $\alpha$. The first term on the right-hand
side (RHS) represents the outflow from the asset. This constraint has
additional terms for producers, consumption, and storage assets, which will be
explained later. Users can choose the constraint sense ($=$, $\leq$, or
$\geq$) based on the asset’s balance behavior. The constraint assumes a linear
transfer function between input and output flow, which is the default method
for the first version. However, we will also include an alternative method
that uses a convex piece-wise linear function.

####  3.1.1 Balance method for conversion assets

In this case, equation (1) applies to conversion assets $\mathcal{A}_{cv}$
using the balance asset method ($\mathcal{A}_{b}$). It only considers the
inflows and outflows of the asset, and uses efficiency ($\eta$) as a
conversion factor to calculate the outgoing flow for each incoming flow.

$\displaystyle\sum_{\alpha\in\mathcal{A}^{in}(a)}\quad\sum_{\tau(\alpha)\in\Omega_{\alpha
p\tau}}M_{\alpha p\tau}f_{y\alpha
ak\tau}=\sum_{\alpha\in\mathcal{A}^{out}(a)}\quad\sum_{\tau(\alpha)\in\Omega_{\alpha
p\tau}}M_{\alpha p\tau}\frac{f_{ya\alpha k\tau}}{\eta_{a\alpha}}$
$\displaystyle\qquad\forall{y}\in\mathcal{Y}_{m},\forall{a}\in\mathcal{A}_{cv}\cap\mathcal{A}_{b},\forall{k}\in\mathcal{K},\forall{p}\in\mathcal{P}$

Let’s take a look at an example of a fuel cell $a$ for a twelve-hourly
representative period ($k$); see Figure (3):

Figure 3: Conversion asset example - Fuel cell

  * •

Input: Hydrogen asset ($H_{2}$) with a resolution of 3h, i.e.,
$\tau(\alpha=H_{2})=\\{1,2,3,4\\}$

  * •

Output: Electricity asset ($E$) with a resolution of 1h, i.e.,
$\tau(\alpha=E)=\\{1,2,...12\\}$, and $\eta_{a,E}=40\%$

  * •

Output: Heat asset ($H$) with a resolution of 4h, i.e.,
$\tau(\alpha=H)=\\{1,2,3\\}$, and $\eta_{a,H}=20\%$

  * •

$\mathcal{A}^{in}(a)=\\{H_{2}\\}$ and $\mathcal{A}^{out}(a)=\\{E,H_{2}\\}$

The resolution determines how often we have variables for each sector (i.e.,
hydrogen, electricity, and heat). For instance, the hydrogen asset flow is
measured by four variables ($f_{y,H_{2},a,k,1}$, $f_{y,H_{2},a,k,2}$,
$f_{y,H_{2},a,k,3}$, and $f_{y,H_{2},a,k,4}$), while the total electricity
generated per hour is also taken into account ($f_{y,a,E,k,1}$,
$f_{y,a,E,k,2}$ … $f_{y,a,E,k,12}$). Similarly, the fuel cell’s flow to the
heat asset is measured using three variables ($f_{y,a,H,k,1}$,
$f_{y,a,H,k,2}$, and $f_{y,a,H,k,3}$). Please note that the resolutions are
not multiples of each other and that the higher the resolution, the fewer
variables. The mapping of variables to the representative period ($k$) is
shown in Table 1. Furthermore, we have:

$p$ | Periods in | $H_{2}$ | $E$ | $H$
---|---|---|---|---
| the $k$ | variables | variables | variables
1 | 1 | $1/3f_{y,H_{2},a,k,1}$ | $f_{y,a,E,k,1}$ | $1/4f_{y,a,H,k,1}$
| 2 | $1/3f_{y,H_{2},a,k,1}$ | $f_{y,a,E,k,2}$ | $1/4f_{y,a,H,k,1}$
| 3 | $1/3f_{y,H_{2},a,k,1}$ | $f_{y,a,E,k,3}$ | $1/4f_{y,a,H,k,1}$
| 4 | $1/3f_{y,H_{2},a,k,2}$ | $f_{y,a,E,k,4}$ | $1/4f_{y,a,H,k,1}$
2 | 5 | $1/3f_{y,H_{2},a,k,2}$ | $f_{y,a,E,k,5}$ | $1/4f_{y,a,H,k,2}$
| 6 | $1/3f_{y,H_{2},a,k,2}$ | $f_{y,a,E,k,6}$ | $1/4f_{y,a,H,k,2}$
| 7 | $1/3f_{y,H_{2},a,k,3}$ | $f_{y,a,E,k,7}$ | $1/4f_{y,a,H,k,2}$
| 8 | $1/3f_{y,H_{2},a,k,3}$ | $f_{y,a,E,k,8}$ | $1/4f_{y,a,H,k,2}$
3 | 9 | $1/3f_{y,H_{2},a,k,3}$ | $f_{y,a,E,k,9}$ | $1/4f_{y,a,H,k,3}$
| 10 | $1/3f_{y,H_{2},a,k,4}$ | $f_{y,a,E,k,10}$ | $1/4f_{y,a,H,k,3}$
| 11 | $1/3f_{y,H_{2},a,k,4}$ | $f_{y,a,E,k,11}$ | $1/4f_{y,a,H,k,3}$
| 12 | $1/3f_{y,H_{2},a,k,4}$ | $f_{y,a,E,k,12}$ | $1/4f_{y,a,H,k,3}$
Table 1: Mapping between variables and periods for the fuel-cell $a$ example

  * •

The index $p$ represents the maximum resolution of connections to the asset.
In this case, $p=\\{1,2,3\\}$ and each value of $p$ will span 4h. Table 1 uses
three different colours to highlight the periods in the $k$ that are
considered for each value of $p$.

  * •

The matrices $M_{\alpha p\tau}$ can be obtained from Table 1. These matrices
will map the values from the auxiliary index $p$ to each time resolution.

$M_{H_{2},p,\tau(H_{2})}=\begin{array}[]{l}{\color[rgb]{.75,.75,.75}p=1}\\\
{\color[rgb]{.75,.75,.75}p=2}\\\
{\color[rgb]{.75,.75,.75}p=3}\end{array}\left(\begin{array}[]{cccc}{\color[rgb]{0,0,1}1}&{\color[rgb]{0,0,1}1/3}&0&0\\\
0&{\color[rgb]{0,1,0}2/3}&{\color[rgb]{0,1,0}2/3}&0\\\
0&0&{\color[rgb]{0.5,0,1}1/3}&{\color[rgb]{0.5,0,1}1}\\\ \end{array}\right)$

$M_{E,p,\tau(E)}=\begin{array}[]{l}{\color[rgb]{.75,.75,.75}p=1}\\\
{\color[rgb]{.75,.75,.75}p=2}\\\
{\color[rgb]{.75,.75,.75}p=3}\end{array}\left(\begin{array}[]{cccccccccccc}{\color[rgb]{0,0,1}1}&{\color[rgb]{0,0,1}1}&{\color[rgb]{0,0,1}1}&{\color[rgb]{0,0,1}1}&0&0&0&0&0&0&0&0\\\
0&0&0&0&{\color[rgb]{0,1,0}1}&{\color[rgb]{0,1,0}1}&{\color[rgb]{0,1,0}1}&{\color[rgb]{0,1,0}1}&0&0&0&0\\\
0&0&0&0&0&0&0&0&{\color[rgb]{0.5,0,1}1}&{\color[rgb]{0.5,0,1}1}&{\color[rgb]{0.5,0,1}1}&{\color[rgb]{0.5,0,1}1}\\\
\end{array}\right)$

$M_{H,p,\tau(H)}=\begin{array}[]{l}{\color[rgb]{.75,.75,.75}p=1}\\\
{\color[rgb]{.75,.75,.75}p=2}\\\
{\color[rgb]{.75,.75,.75}p=3}\end{array}\left(\begin{array}[]{ccc}{\color[rgb]{0,0,1}1}&0&0\\\
0&{\color[rgb]{0,1,0}1}&0\\\ 0&0&{\color[rgb]{0.5,0,1}1}\\\
\end{array}\right)$

The coefficients in the matrices determine the fraction of each variable taken
into account for each constraint. For example, when $p=1$, the complete value
of $f_{y,H_{2},a,k,1}$ is considered, but only $1/3$ of $f_{y,H_{2},a,k,2}$ is
taken into account since $p$ represents a four-hour interval (refer to Table
1). Similarly, when $p=2$, $2/3$ of $f_{y,H_{2},k,2}$ and $2/3$ of
$f_{y,H_{2},a,k,3}$ are considered, and so on.

$p=1$:

$f_{y,H_{2},a,k,1}+\frac{1}{3}f_{y,H_{2},a,k,2}=\frac{f_{y,a,E,k,1}}{0.4}+\frac{f_{y,a,E,k,2}}{0.4}+\frac{f_{y,a,E,k,3}}{0.4}+\frac{f_{y,a,E,k,4}}{0.4}+\frac{f_{y,a,H,k,1}}{0.2}$

$p=2$:

$\frac{2}{3}f_{y,H_{2},a,k,2}+\frac{2}{3}f_{y,H_{2},a,k,3}=\frac{f_{y,a,E,k,5}}{0.4}+\frac{f_{y,a,E,k,6}}{0.4}+\frac{f_{y,a,E,k,7}}{0.4}+\frac{f_{y,a,E,k,8}}{0.4}+\frac{f_{y,a,H,k,2}}{0.2}$

$p=3$:

$\frac{1}{3}f_{y,H_{3},a,k,1}+f_{y,H_{2},a,k,4}=\frac{f_{y,a,E,k,9}}{0.4}+\frac{f_{y,a,E,k,10}}{0.4}+\frac{f_{y,a,E,k,11}}{0.4}+\frac{f_{y,a,E,k,12}}{0.4}+\frac{f_{y,a,H,k,3}}{0.2}$

To define the sets $\Omega_{\alpha p\tau}$, we use a tuple of three elements
derived from the matrices $M_{\alpha p\tau}$ so that the positive values
compound the elements in the set.

$\Omega_{\alpha p\tau}=\\{(\alpha,p,\tau)|M_{\alpha p\tau}>0\\}$

####  3.1.2 Balance method for production assets

Production assets are elements that produce energy from a resource and
distribute it to one or more energy assets in the network. These assets can
include solar, wind, nuclear, and energy imports. When the method is applied
to production assets, equation (1) can be simplified into the following
expression.

$\displaystyle\sum_{\tau(a)\in\Omega_{ap\tau}}M_{ap\tau}P_{yak\tau}\left\\{\begin{array}[]{l}=\\\
\geqslant\\\
\leqslant\end{array}\right\\}\sum_{\alpha\in\mathcal{A}^{out}(a)}\quad\sum_{\tau(\alpha)\in\Omega_{\alpha
p\tau}}M_{\alpha p\tau}f_{ya\alpha k\tau}$
$\displaystyle\qquad\forall{y}\in\mathcal{Y}_{m},\forall{a}\in\mathcal{A}_{p}\cap\mathcal{A}_{b},\forall{k}\in\mathcal{K},\forall{p}\in\mathcal{P}$

On the LHS, we have the representation of the production resource
availability, while on the RHS, we have the flows that are directed towards
other assets. It is important to note that the efficiency is equal to one
($\eta_{\alpha a}=1$) due to the absence of any energy conversion in the
production assets. This approach can prove helpful for various production
assets, but it’s particularly relevant for energy imports. As an example, take
a look at Figure 4, which showcases a methane source that can provide fuel to
a gas-fired plant and a $H_{2}$ generator like a steam-methane reformation
(SMR), having a twelve-hourly representative period ($k$) again:

  * •

Production: Methane ($CH_{4}$) with a resolution of 6h, i.e.,
$\tau(\alpha=CH_{4})=\\{1,2\\}$

  * •

Output: Gas turbine ($GT$) with a resolution of 3h, i.e.,
$\tau(\alpha=GT)=\\{1,2,3,4\\}$

  * •

Output: Steam-methane reformation ($SMR$) with a resolution of 4h, i.e.,
$\tau(\alpha=SMR)=\\{1,2,3\\}$

  * •

$\mathcal{A}^{out}(CH_{4})=\\{GT,SMR\\}$

Figure 4: Production asset example - Methane producer

As shown in the previous section, we can map variables to the representative
period ($k$) using Table 2. The index $p$ represents the maximum resolution of
connections to the asset, with $p$ being either 1 or 2, each representing a
6-hour span. Table 2 uses two colours to highlight the periods in the $p$
considered for each $p$ value. The $M_{\alpha p\tau}$ matrices can be derived
from the table values.

$p$ | Periods in | $CH_{4}$ | $GT$ | $SMR$
---|---|---|---|---
| the $k$ | profile values | variables | variables
1 | 1 | $1/6P_{y,CH_{4},k,1}$ | $1/3f_{y,CH_{4},GT,k,1}$ | $1/4f_{y,CH_{4},SMR,k,1}$
| 2 | $1/6P_{y,CH_{4},k,1}$ | $1/3f_{y,CH_{4},GT,k,1}$ | $1/4f_{y,CH_{4},SMR,k,1}$
| 3 | $1/6P_{y,CH_{4},k,1}$ | $1/3f_{y,CH_{4},GT,k,1}$ | $1/4f_{y,CH_{4},SMR,k,1}$
| 4 | $1/6P_{y,CH_{4},k,1}$ | $1/3f_{y,CH_{4},GT,k,2}$ | $1/4f_{y,CH_{4},SMR,k,1}$
| 5 | $1/6P_{y,CH_{4},k,1}$ | $1/3f_{y,CH_{4},GT,k,2}$ | $1/4f_{y,CH_{4},SMR,k,2}$
| 6 | $1/6P_{y,CH_{4},k,1}$ | $1/3f_{y,CH_{4},GT,k,2}$ | $1/4f_{y,CH_{4},SMR,k,2}$
2 | 7 | $1/6P_{y,CH_{4},k,2}$ | $1/3f_{y,CH_{4},GT,k,3}$ | $1/4f_{y,CH_{4},SMR,k,2}$
| 8 | $1/6P_{y,CH_{4},k,2}$ | $1/3f_{y,CH_{4},GT,k,3}$ | $1/4f_{y,CH_{4},SMR,k,2}$
| 9 | $1/6P_{y,CH_{4},k,2}$ | $1/3f_{y,CH_{4},GT,k,3}$ | $1/4f_{y,CH_{4},SMR,k,3}$
| 10 | $1/6P_{y,CH_{4},k,2}$ | $1/3f_{y,CH_{4},GT,k,4}$ | $1/4f_{y,CH_{4},SMR,k,3}$
| 11 | $1/6P_{y,CH_{4},k,2}$ | $1/3f_{y,CH_{4},GT,k,4}$ | $1/4f_{y,CH_{4},SMR,k,3}$
| 12 | $1/6P_{y,CH_{4},k,2}$ | $1/3f_{y,CH_{4},GT,k,4}$ | $1/4f_{y,CH_{4},SMR,k,3}$
Table 2: Mapping between variables and periods for the methane imports
$CH_{4}$ example

$M_{CH_{4},k,\tau(CH_{4})}=\begin{array}[]{l}{\color[rgb]{.75,.75,.75}p=1}\\\
{\color[rgb]{.75,.75,.75}p=2}\end{array}\left(\begin{array}[]{cc}{\color[rgb]{0,0,1}1}&0\\\
0&{\color[rgb]{0,1,0}1}\\\ \end{array}\right)$

$M_{GT,k,\tau(GT)}=\begin{array}[]{l}{\color[rgb]{.75,.75,.75}p=1}\\\
{\color[rgb]{.75,.75,.75}p=2}\end{array}\left(\begin{array}[]{cccc}{\color[rgb]{0,0,1}1}&{\color[rgb]{0,0,1}1}&0&0\\\
0&0&{\color[rgb]{0,1,0}1}&{\color[rgb]{0,1,0}1}\\\ \end{array}\right)$

$M_{SMR,k,\tau(SMR)}=\begin{array}[]{l}{\color[rgb]{.75,.75,.75}p=1}\\\
{\color[rgb]{.75,.75,.75}p=2}\end{array}\left(\begin{array}[]{ccc}{\color[rgb]{0,0,1}1}&{\color[rgb]{0,0,1}1/2}&0\\\
0&{\color[rgb]{0,1,0}1/2}&{\color[rgb]{0,1,0}1}\\\ \end{array}\right)$

In this example, the equation (1) can be used with different options to
indicate minimum production, maximum production, or fixed production. We will
use the maximum production limit as an example in this case.

$p=6$:

$\displaystyle P_{y,CH_{4},k,1}\geq$ $\displaystyle
f_{y,CH_{4},GT,k,1}+f_{y,CH_{4},GT,k,2}+f_{y,CH_{4},SMR,k,1}+\frac{1}{2}f_{y,CH_{4},SMR,k,2}$

$p=12$:

$\displaystyle P_{y,CH_{4},k,2}\geq$ $\displaystyle
f_{y,CH_{4},GT,k,3}+f_{y,CH_{4},GT,k,4}+\frac{1}{2}f_{y,CH_{4},SMR,k,2}+f_{y,CH_{4},SMR,k,3}$

This method is beneficial for representing production assets that are linked
to multiple assets with varying time resolutions, as illustrated in the
example. However, all production assets don’t need to utilize this method.
There are other specific methods to restrict the flows from the production
asset to the grid if there is only one connection and the production assets
can offer reserves, as discussed in Section 3.2.

####  3.1.3 Balance method for consumption assets

Consumption assets are elements that demand energy from one or more energy
assets in the network. These assets can include energy demands and exports.
When the method is applied to consumption assets, equation (1) can be
simplified into the following expression.

$\displaystyle\sum_{\alpha\in\mathcal{A}^{in}(a)}\quad\sum_{\tau(\alpha)\in\Omega_{\alpha
p\tau}}M_{\alpha p\tau}f_{y\alpha ak\tau}\left\\{\begin{array}[]{l}=\\\
\geqslant\\\
\leqslant\end{array}\right\\}\sum_{\tau(a)\in\Omega_{ap\tau}}M_{ap\tau}D_{yak\tau}$
$\displaystyle\qquad\forall{y}\in\mathcal{Y}_{m},\forall{a}\in\mathcal{A}_{c}\cap\mathcal{A}_{b},\forall{k}\in\mathcal{K},\forall{p}\in\mathcal{P}$

On the LHS, we have the representation of the flows that are directed towards
the consumption asset, while on the RHS, we have the consumption demand
profile. This approach can prove helpful for various consumption assets, but
it’s particularly relevant for energy exports. As an example, take a look at
Figure 5, which showcases electricity exports that can take their demand from
two nodes/countries:

  * •

Input: Node 1 ($N1$) with a resolution of 1h, i.e.,
$\tau(\alpha=N1)=\\{1,2...12\\}$

  * •

Input: Node 2 ($N2$) with a resolution of 4h, i.e.,
$\tau(\alpha=N2)=\\{1,2,3\\}$

  * •

Consumption: Electricity exports ($E$) with a resolution of 6h, i.e.,
$\tau(\alpha=E)=\\{1,2\\}$

  * •

$\mathcal{A}^{in}(E)=\\{N1,N2\\}$

Figure 5: Consumption asset example - Electricity exports

Table 2 maps the variables to the periods in the representative period ($k$).
The $M_{\alpha p\tau}$ matrices can be derived from the table values, such as
before.

$p$ | Periods in | $N1$ | $N2$ | $D$
---|---|---|---|---
| the $k$ | variables | variables | profile values
1 | 1 | $f_{y,N1,E,k,1}$ | $1/4f_{y,N2,E,k,1}$ | $1/6D_{y,E,k,1}$
| 2 | $f_{y,N1,E,k,2}$ | $1/4f_{y,N2,E,k,1}$ | $1/6D_{y,E,k,1}$
| 3 | $f_{y,N1,E,k,3}$ | $1/4f_{y,N2,E,k,1}$ | $1/6D_{y,E,k,1}$
| 4 | $f_{y,N1,E,k,4}$ | $1/4f_{y,N2,E,k,1}$ | $1/6D_{y,E,k,1}$
| 5 | $f_{y,N1,E,k,5}$ | $1/4f_{y,N2,E,k,2}$ | $1/6D_{y,E,k,1}$
| 6 | $f_{y,N1,E,k,6}$ | $1/4f_{y,N2,E,k,2}$ | $1/6D_{y,E,k,1}$
2 | 7 | $f_{y,N1,E,k,7}$ | $1/4f_{y,N2,E,k,2}$ | $1/6D_{y,E,k,2}$
| 8 | $f_{y,N1,E,k,8}$ | $1/4f_{y,N2,E,k,2}$ | $1/6D_{y,E,k,2}$
| 9 | $f_{y,N1,E,k,9}$ | $1/4f_{y,N2,E,k,3}$ | $1/6D_{y,E,k,2}$
| 10 | $f_{y,N1,E,k,10}$ | $1/4f_{y,N2,E,k,3}$ | $1/6D_{y,E,k,2}$
| 11 | $f_{y,N1,E,k,11}$ | $1/4f_{y,N2,E,k,3}$ | $1/6D_{y,E,k,2}$
| 12 | $f_{y,N1,E,k,12}$ | $1/4f_{y,N2,E,k,3}$ | $1/6D_{y,E,k,2}$
Table 3: Mapping between variables and periods for the electricity exports $E$
example

$M_{N1,k,\tau(N1)}=\begin{array}[]{l}{\color[rgb]{.75,.75,.75}p=1}\\\
{\color[rgb]{.75,.75,.75}p=2}\end{array}\left(\begin{array}[]{cccccccccccc}{\color[rgb]{0,0,1}1}&{\color[rgb]{0,0,1}1}&{\color[rgb]{0,0,1}1}&{\color[rgb]{0,0,1}1}&{\color[rgb]{0,0,1}1}&{\color[rgb]{0,0,1}1}&0&0&0&0&0&0\\\
0&0&0&0&0&0&{\color[rgb]{0,1,0}1}&{\color[rgb]{0,1,0}1}&{\color[rgb]{0,1,0}1}&{\color[rgb]{0,1,0}1}&{\color[rgb]{0,1,0}1}&{\color[rgb]{0,1,0}1}\\\
\end{array}\right)$

$M_{N2,k,\tau(N2)}=\begin{array}[]{l}{\color[rgb]{.75,.75,.75}p=1}\\\
{\color[rgb]{.75,.75,.75}p=2}\end{array}\left(\begin{array}[]{ccc}{\color[rgb]{0,0,1}1}&{\color[rgb]{0,0,1}1/2}&0\\\
0&{\color[rgb]{0,1,0}1/2}&{\color[rgb]{0,1,0}1}\\\ \end{array}\right)$

$M_{D,k,\tau(D)}=\begin{array}[]{l}{\color[rgb]{.75,.75,.75}p=1}\\\
{\color[rgb]{.75,.75,.75}p=2}\end{array}\left(\begin{array}[]{cc}{\color[rgb]{0,0,1}1}&0\\\
0&{\color[rgb]{0,1,0}1}\\\ \end{array}\right)$

In this example, the equation (1) can be used with different options to
indicate minimum consumption, maximum consumption, or fixed consumption. We
will use the minimum consumption limit as an example in this case.

$p=6$:

$\displaystyle f_{y,N1,E,k,1}+f_{y,N1,E,k,2}+f_{y,N1,E,k,3}+f_{y,N1,E,k,4}$
$\displaystyle+f_{y,N1,E,k,5}+f_{y,N1,E,k,6}+f_{y,N2,E,k,1}+\frac{1}{2}f_{y,N2,E,k,2}\leq
D_{y,E,k,1}$

$p=12$:

$\displaystyle f_{y,N1,E,k,7}+f_{y,N1,E,k,8}+f_{y,N1,E,k,9}+f_{y,N1,E,k,10}$
$\displaystyle+f_{y,N1,E,k,11}+f_{y,N1,E,k,12}+\frac{1}{2}f_{y,N2,E,k,2}+f_{y,N2,E,k,3}\leq
D_{y,E,k,2}$

This approach helps depict consumption assets that are associated with several
assets having different time resolutions, as shown in the example. Rather than
using a node to combine the flows, we can simplify the connections by
directing all the flows towards the consumption asset that has a demand, which
ultimately reduces the number of variables.

####  3.1.4 Balance method for transport assets

The RHS of the equation represents all the production that goes into the
transport asset. The first term on the LHS is the flow going out of the
transport asset. Since transport assets do not convert energy, the efficiency
(represented by the symbol $\eta$) equals 1. Additionally, the constraint is
an equality constraint because there is no production or consumption of energy
at the transport assets. We can represent, for example, a nodal balance in an
intermediate point, such as in a transhipment problem.

$\displaystyle\sum_{\alpha\in\mathcal{A}^{in}(a)}\quad\sum_{\tau(\alpha)\in\Omega_{\alpha
p\tau}}M_{\alpha p\tau}f_{y\alpha
ak\tau}=\sum_{\alpha\in\mathcal{A}^{out}(a)}\quad\sum_{\tau(\alpha)\in\Omega_{\alpha
p\tau}}M_{\alpha p\tau}\frac{f_{ya\alpha k\tau}}{\eta_{a\alpha}}$
$\displaystyle\qquad\forall{y}\in\mathcal{Y}_{m},\forall{a}\in\mathcal{A}_{t}\cap\mathcal{A}_{b},\forall{k}\in\mathcal{K},\forall{p}\in\mathcal{P}$

####  3.1.5 Balance method for storage assets

For storage assets using a balance method, the LHS represents the charging
flows into the storage asset, typically only one. The LHS also has the change
in the storage level of the asset and inflows (if applicable). On the RHS,
there are the discharging flows (also typically only one). Each equation
component may have a different time resolution, allowing for situations where
the discharging and charging decision variables have a lower resolution than
the desired storage level variables. For example, hourly charging/discharging
decisions with a 24-hour resolution for the storage level help represent
pumped-hydro storage units.

$\displaystyle\sum_{\alpha\in\mathcal{A}^{in}(a)}\quad\sum_{\tau(\alpha)\in\Omega_{\alpha
p\tau}}M_{\alpha p\tau}f_{y\alpha ak\tau}$
$\displaystyle+\sum_{\tau(a)\in\Omega_{ap\tau}}M_{ap\tau}\cdot\left(s_{yak\tau}-s_{yak,\tau-1}+S_{yak\tau}^{I}\right)=$
$\displaystyle\sum_{\alpha\in\mathcal{A}^{out}(a)}\quad\sum_{\tau(\alpha)\in\Omega_{\alpha
p\tau}}M_{\alpha p\tau}\frac{f_{ya\alpha k\tau}}{\eta_{a\alpha}}$
$\displaystyle\qquad\forall{y}\in\mathcal{Y}_{m},\forall{a}\in\mathcal{A}_{s}\cap\mathcal{A}_{b},\forall{k}\in\mathcal{K},\forall{p}\in\mathcal{P}$

###  3.2 Operation limit constraints

Constraints (2) and (3) represent the maximum and minimum output flow for each
asset. Note that some assets may have only one output flow (e.g. renewable
assets), while others may have several (e.g., fuel cells). Unlike the balance
constraint in Section 3.1, these constraints allow for the consideration of
reserve ($r$) and unit commitment ($u$) variables. Furthermore, the parameters
of the profile, which are represented by $\tilde{F}^{max}$ and
$\tilde{F}^{max}$, can be set with a temporal resolution that is separate from
that of the variables. For example, the profiles could have hourly values
while the variable definition is set every 4 hours. To account for this, the
formulation includes a mapping matrix that multiplies the profiles, as seen in
the balance constraint.

$\displaystyle\sum_{\alpha\in\mathcal{A}^{out}(a)}f_{ya\alpha
kt}+\sum_{r\in\mathcal{R}_{a}^{+}}r_{yarkt}\leq$
$\displaystyle\overline{F}_{a}u_{yakt}\sum_{p\in\Omega_{apt}}M_{apt}\tilde{F}_{yakp}^{max}$
(2)
$\displaystyle\forall{y}\in\mathcal{Y}_{m},\forall{a}\in\mathcal{A},\forall{k}\in\mathcal{K},\forall{t}\in\tau(a)$

$\displaystyle\sum_{\alpha\in\mathcal{A}^{out}(a)}f_{ya\alpha
kt}-\sum_{r\in\mathcal{R}_{a}^{-}}r_{yarkt}\geq$
$\displaystyle\overline{F}_{a}u_{yakt}\sum_{p\in\Omega_{apt}}M_{apt}\tilde{F}_{yakp}^{min}$
(3)
$\displaystyle\forall{y}\in\mathcal{Y}_{m},\forall{a}\in\mathcal{A},\forall{k}\in\mathcal{K},\forall{t}\in\tau(a)$

As a general feature, the minimum output flow constraint helps represent
minimum production levels without a unit commitment method, e.g., a nuclear
plant with a minimum stable load of 80%; however, it can go up to 100% of its
capacity if needed. Notice that for assets with a unit commitment method, then
$\tilde{F}^{min}=0$, and for assets with a transport method, then
$\tilde{F}^{min}=-1$. The upwards/downwards spinning reserve variables will
appear or not in the constraint depending on whether the asset has the reserve
method assigned or not (e.g., transport assets do not have reserve variables,
whereas power plant units might have the possibility to provide different
types of reserves).

###  3.3 Reserve requirement constraints

The reserve limit ensures that all assets capable of providing reserve type
$r$, denoted by $\Omega_{ar}$, must sum up to at least the minimum reserve
requirement ($R$) at each time period ($p$) in the representative period
($k$). Note that the constraint is defined on the temporal definition of the
parameter $R$ (i.e., $p$). Therefore, we need to multiply the variable for
reserve with the mapping matrix between the temporal definitions of the
variables and the parameter.

$\sum_{a\in\Omega_{ar}}\quad\sum_{\tau(\alpha)\in\Omega_{\alpha
p\tau}}M_{\alpha p\tau}r_{yark\tau}\geq
R_{yrkp}\quad\forall{y}\in\mathcal{Y}_{m},\forall{r}\in\mathcal{R},\forall{k}\in\mathcal{K},\forall{p}\in\mathcal{P}$
(4)

###  3.4 Investment methods

The subset $\mathcal{A}_{i}$ includes all assets with an investment method.
Therefore, assets can have two types of constraints: with and without
investment.

####  3.4.1 Asset with investment

The following constraints relate the available energy asset units with the
investment decision variables and limit the maximum investment. For instance,
constraint (5) imposes that the available energy asset units at year $y$ must
be lower or equal to those available at the beginning plus those built in the
previous years. In addition, constraint (6) ensure that the cumulative
investment is lower or equal to the investment potential
$\overline{F}_{ya}^{max}$ at a specific year $y$ for each energy asset.

$u_{yakt}\leq\overline{U}_{ya}+\sum_{\psi=\max(y-LT_{a}+1,0)}^{y}\overline{u}_{\psi,a}\quad\forall{y}\in\mathcal{Y}_{m},\forall{a}\in\mathcal{A}_{i},\forall{k}\in\mathcal{K},\forall{t}\in\tau(a)$
(5)

$\overline{F}_{a}\sum_{\psi=\max(y-LT_{a}+1,0)}^{y}\overline{u}_{\psi
a}\leq\overline{F}_{ya}^{max}\quad\forall{y}\in\mathcal{Y}_{m},\forall{a}\in\mathcal{A}_{i}$
(6)

####  3.4.2 Asset without investment

For assets excluded from the investment method ($a\notin\mathcal{A}_{i}$), the
available units are only limited by their initial amount and $\overline{u}$ is
set to 0.

$u_{yakt}\leq\overline{U}_{ya}\quad\forall{y}\in\mathcal{Y}_{m},\forall{a}\notin\mathcal{A}_{i},\forall{k}\in\mathcal{K},\forall{t}\in\tau(a)$
(7)

##  4 Objective function

The objective function is to minimise the energy system’s total investment and
operational cost.

$\displaystyle\min\quad C^{I}+C^{O}$ (8) $\displaystyle
C^{I}=\sum_{y\in\mathcal{Y}_{m}}\sum_{a\in\mathcal{A}_{i}}\frac{1}{(1+IR)^{y}}\cdot(C_{ya}^{T}-SV_{ya})\cdot\overline{F}_{a}\overline{u}_{ya}$
$\displaystyle C^{O}=\sum_{y\in\mathcal{Y}_{m}}\frac{1}{(1+IR)^{y}}\cdot
W_{y}^{m}\cdot\sum_{a\in\mathcal{A}}C_{ya}^{op}\cdot\sum_{\alpha\in\mathcal{A}^{out}(a)}\sum_{k\in\mathcal{K}}W_{yk}^{op}\cdot\sum_{t\in\mathcal{T}}\sum_{\tau(\alpha)\in\Omega_{\alpha
t\tau}}M_{\alpha t\tau}f_{ya\alpha k\tau}$

The term $C^{I}$ represents the investment costs of energy assets. It is
multiplied by the interest rate $IR$ for each milestone year
($y\in\mathcal{Y}_{m}$). It considers the total investment cost $C_{y}^{T}$
and salvage value $SV_{y}$ for the end-of-horizon effect.

The term $C^{O}$ represents the total variable production cost $C_{ya}^{op}$,
also known as operational cost. The model will include other costs like start-
up costs, no-load costs, and shut-down costs once unit commitment methods are
defined. Operational cost includes the interest rate $IR$, the weight of the
milestone year $W_{y}^{m}$, and the weight of the representative period
$W_{yk}^{op}$. Each flow in the model is multiplied by the mapping matrix to
account for the variable decision’s suitable duration due to the temporal
structure’s flexibility.

The article referenced by [7] explores various methods for modeling multiyear
investments in energy systems, outlining the pros and cons of each one.

##  5 GitHub repository

The Tulipa Energy Model is implemented in Julia [8] using the JuMP [9] and
Graphs [10] packages. The code is available on GitHub and is based on the
formulation presented in this document.

GitHub repository: https://github.com/TulipaEnergy/TulipaEnergyModel.jl

## References

  * [1] Antonio J. Conejo, Luis Baringo Morales, S. Kazempour and Afzal S. Siddiqui  “Investment in Electricity Generation and Transmission: Decision Making under Uncertainty”  Springer Cham, 2016, pp. 384  DOI: https://doi.org/10.1007/978-3-319-29501-5
  * [2] Douglas B. West  “Introduction to Graph Theory”  Prentice Hall, 2000
  * [3] TNO  “ESDL, de Energy System Description Language” Accessed on August 7, 2023, 2021  URL: https://www.esdl.nl/
  * [4] Niina Helistö et al.  “Backbone—An Adaptable Energy Systems Modelling Framework”  In _Energies_ 12.17, 2019  DOI: 10.3390/en12173388
  * [5] Maren Ihlemann et al.  “SpineOpt: A flexible open-source energy system modelling framework”  In _Energy Strategy Reviews_ 43, 2022, pp. 100902  DOI: https://doi.org/10.1016/j.esr.2022.100902
  * [6] Zhi Gao et al.  “Fully Flexible Temporal Resolution for Energy System Optimization Using Index Aggregation” pre-print, 2023
  * [7] Diego A. Tejada-Arango  “Multi-year Investment Modelling in Energy Systems”, 2023  arXiv:2307.16842 [math.OC]
  * [8] Jeff Bezanson, Alan Edelman, Stefan Karpinski and Viral B Shah  “Julia: A fresh approach to numerical computing”  In _SIAM Review_ 59.1 SIAM, 2017, pp. 65–98  DOI: 10.1137/141000671
  * [9] Miles Lubin et al.  “JuMP 1.0: Recent improvements to a modeling language for mathematical optimization”  In _Mathematical Programming Computation_ , 2023  DOI: 10.1007/s12532-023-00239-3
  * [10] James Fairbanks et al.  “JuliaGraphs/Graphs.jl: an optimized graphs package for the Julia programming language”, 2021  URL: https://github.com/JuliaGraphs/Graphs.jl/


# Non-standard Neutrino and $Z^{\prime}$ Interactions at the FASER$\nu$ and
the LHC

Kingman Cheunga,b,d, C.J. Ousepha,b, TseChun Wangc a Department of Physics,
National Tsing Hua University, Hsinchu 300, Taiwan
b Center for Theory and Computation, National Tsing Hua University, Hsinchu
300, Taiwan
c Physics Division, National Center for Theoretical Sciences, Taipei 10617,
Taiwan
d Division of Quantum Phases and Devices, School of Physics, Konkuk
University, Seoul 143-701, Republic of Korea

(August 28, 2024)

###### Abstract

We study the impact of non-standard neutrino interactions in the context of a
new gauge boson $Z^{\prime}$ in neutral-current deep-inelastic scattering
performed in ForwArd Search ExpeRiment-$\nu$ (FASER$\nu$) and in monojet
production at the Large Hadron Collider (LHC). We simulate the neutral-current
deep-inelastic neutrino-nucleon scattering $\nu N\rightarrow\nu N$ at
FASER$\nu$ in the presence of an additional $Z^{\prime}$ boson, and estimate
the anticipated sensitivities to the gauge coupling in a wide range of
$Z^{\prime}$ mass. At the LHC, we study the effect of $Z^{\prime}$ on monojet
production, which can be enhanced in regions with large missing transverse
momenta. We then use the recent results from ATLAS with an integrated
luminosity of 139 fb-1 to improve the limits on the gauge coupling of
$Z^{\prime}$. We interpret such limits on $Z^{\prime}$ gauge couplings as
bounds on effective non-standard neutrino interactions. We show that the
FASER$\rm\nu$ and the LHC results cover the medium and high energy scales,
respectively, and complement one another.

††preprint: NCTS-HEPAP/2102

##  I Introduction

Determining the properties of neutrinos is important in the context of
searching for physics beyond the standard model and understanding the
universe. There are several motivations of these measurements. First, the
neutrino oscillation indicates that the neutrinos are light but not massless,
in contrast to the prediction by the standard model (SM). This rock-solid fact
motivates physicists to search for any other neutrino properties, which are
beyond the standard model (BSM). Second, the precision of parameters in the
neutrino physics needs improvements. This is because neutrinos are hard to be
detected, compared to the other charged fermions. In these measurements,
statistical uncertainties are usually dominating over other sources of
uncertainties. Therefore, even if there are truly some BSM features, our
precision might be not good enough to catch these features. And finally,
neutrino is one of the most abundant particles in the universe since the Big
Bang nucleosynthesis (BBN). As a result, the property of neutrinos is an
essential factor in the evolution of the universe, e.g. the neutrino-neutrino
self interaction might affect the measurement of Hubble constant $H_{0}$
Blinov:2019gcj ; Cyr-Racine:2013jua ; Kreisch:2019yzn ; Ghosh:2019tab .

To catch the BSM features in neutrino physics, physicists have considered a
variety of experimental configurations. One approach, which has the least
uncertainty, is to observe the neutrinos in collider experiments. Several
proposals have been discussed, such as: Search for Hidden Particle (SHiP)
DeLellis:2015usa ; SHiP:2015vad  and ForwArd Search ExpeRiment (FASER)
FASER:2019dxq ; FASER:2019aik , etc. Both SHiP and FASER are proposed to
detect neutrinos and long-live particles from CERN. In addition to the FASER
main detector, the sub-detector in the front is a $1.2$-ton tungsten detector
– FASER$\nu$ FASER:2019dxq . With high neutrino luminosity ($2\times
10^{11}\nu_{e}$, $6\times 10^{12}\nu_{\mu}$, $4\times 10^{9}\nu_{\tau}$ during
LHC-Run3), FASER$\nu$ provides an optimal window for precision measurements of
neutrino properties of all flavors at the medium-high energy scale ($600$ GeV
to $1$ TeV). This can also be used to search for BSM physics, e.g. new
interactions. The main detector in FASER is now taking data. When the LHC
restarts in 2022. FASER$\nu$ is expected to measure the neutrinos during the
next period of LHC operation from 2022 to 2024.

The most distinct feature of the FASER$\nu$ experiment is the unique energy
range of neutrinos that it can cover. The ICECUBE focuses on very high-energy
neutrinos with energy 10 TeV to 1 PeV, and the LHC covers from hundreds of GeV
to a few TeV. On the other hand, the short- and long-baseline experiments
cover mostly around MeV up to a few GeV. There are no precise measurements of
neutrino scattering in a few tens of GeV to a few hundreds of GeV region.
FASER$\nu$ based on the neutrino flux coming off the LHC opens such a unique
window in this energy range.

One of the simplest extensions to the SM is to add an extra $U(1)$ gauge
symmetry, which results in a new neutral gauge boson $Z^{\prime}$. Such a
$Z^{\prime}$ boson can couple to the SM fermions or simply hidden, depending
on the construction. It is also motivated by some theoretical models (e.g.
dark matter models He:1991qd ; Shepherd:2009sa ). However, without any hint
about the energy scale of this new physics, we can only treat the mass of
$Z^{\prime}$ ($M_{Z^{\prime}}$) as a free parameters. A systematical search in
a wide energy range is therefore very important. As this $Z^{\prime}$ might be
a generator of a new symmetry of flavor such as the $\mu$-$\tau$ symmetry, the
neutrino detection with all flavors is an advantage to test the $Z^{\prime}$
models of this kind. We see that FASER$\nu$, with the capability of
distinguishing the flavors of neutrinos, will play an important role in the
$Z^{\prime}$ search and test for the flavor structure in the $Z^{\prime}$
interactions. We will further investigate these features in this work.

On the other hand, the LHC monojet production can cover effectively the mass
range from a few hundred GeV up to a few TeV. The LHC monojet data can put
stringent constraints on the $Z^{\prime}$ gauge coupling. We use the most
updated monojet data with 139 fb-1 luminosity ATLAS:2021kxv , and obtain the
best limit on the gauge couplings $g_{q}g_{\nu}$, which can be translated to
the effective
$\epsilon(\bar{q}\gamma^{\mu}q)(\bar{\nu}_{L}\gamma_{\mu}\nu_{L})$.
Considerable improvement over previous works is demonstrated here.
Nevertheless, the monojet data is not sensitive to the flavor of the
neutrinos, and therefore the $\epsilon$ is the sum of contributions from all
three flavors, in contrast to low-energy oscillation experiments.

The organization of the work is as follows. In the next section, we briefly
introduce the theoretical aspects of $Z^{\prime}$, the relevant phenomenology,
and the current status. In Sec. III, we show the effects of $Z^{\prime}$
interactions on LHC monojet production and obtain the limits on the effective
NSI. In Sec. IV, we study the sensitivities of $Z^{\prime}$ interactions
achieved at FASER$\nu$. In Sec. V, we show the complementarity of LHC monojet
production and FASER$\nu$ in the coverage of mass range of $Z^{\prime}$.
Finally, we give our conclusions in Sec. VI.

##  II The $Z^{\prime}$ model and Non-standard Neutrino Interactions

Renormalizable interactions of the $Z^{\prime}$ with flavor-conserving quark
and neutrino interactions can be written as

${\cal
L}_{Z^{\prime}}=-\left(g_{\nu}\bar{\nu}\gamma^{\mu}P_{L}\nu\;+\;g_{q}\bar{q}\gamma^{\mu}q\right)\;Z^{\prime}_{\mu}\;.$
(1)

Here we assume that $q=u,d$ have equal coupling strength $g_{q}$ and
$\nu=\nu_{e}.\nu_{\mu},\nu_{\tau}$ have equal strength $g_{\nu}$. In this
simplified $Z^{\prime}$ model, we assume that the coupling strengths to the
left- and right-handed $u,d$ quarks are the same, and so are the coupling
strengths to the three flavors of neutrinos, as production of high-energy
neutrinos is not sensitive to the flavors of neutrinos. Nevertheless, the
results can be easily extended to non-universal coupling strengths.

Although we use a simplified $Z^{\prime}$ model in our working procedures,
there are still a number of existing constraints on general $Z^{\prime}$
models. We briefly discuss in the following.

  1. 1.

The Big Bang Nucleosynthesis (BBN) places constraints on the mass of the boson
$M_{Z^{\prime}}\lesssim 5$ MeV Huang:2017egl .

  2. 2.

Supernova cooling also leads to substantial effects in the observed supernova
neutrino spectrum, which implied the $Z^{\prime}$ coupling to be as small as
$g_{\nu}\sim 10^{-10}$ for $M_{Z^{\prime}}\lesssim 30\;{\rm MeV}$ Dent:2012mx
; Harnik:2012ni .

  3. 3.

The branching ratio for $K^{0}_{L}\rightarrow\pi^{0}Z^{\prime}$ leads to a
bound $g_{q}\lesssim 10^{-8}$ for $M_{Z^{\prime}}=100-200\,{\rm MeV}$
Nelson:1989fx .

  4. 4.

The measurement of $\eta\rightarrow\pi^{0}\gamma\gamma$ gives a bound
$g_{q}\lesssim 10^{-5}-0.01$ for $M_{Z^{\prime}}$ ranging from 200 to 600 MeV
Tulin:2014tya . Other measurements on the branching ratios of
$\eta^{\prime}\rightarrow\pi^{0}\pi^{+}\pi^{-}\gamma$, $\psi\rightarrow
K^{+},K^{-}$, and $\Upsilon\rightarrow$ hardons give a bound $g_{q}\lesssim
0.01-0.1$ for $M_{Z^{\prime}}=0.5,5.5$ and 9.8 GeV.

  5. 5.

BaBar put a constraint on the coupling strength of electron to $Z^{\prime}$
$g_{e}\lesssim 3.3\times 10^{-2}$ from the process
$e^{+}e^{-}\rightarrow\gamma Z^{\prime}$ with $M_{Z^{\prime}}\lesssim 10\,{\rm
GeV}$ BaBar:2008aby ; Essig:2013vha .

  6. 6.

Borexino placed a bound $g_{e,\mu}\lesssim\mathcal{O}(10^{-2})$ for
$M_{Z^{\prime}}\sim 1\,{\rm GeV}$. Furthermore, for a very light $Z^{\prime}$
of mass $M_{Z^{\prime}}\sim 1\,{\rm MeV}$ the constraint becomes more
stringent $g_{e,\mu}\lesssim\mathcal{O}(10^{-5})$ Harnik:2012ni .

###  II.1 Non-standard Neutrino Interactions

It is clear from the above discussion that the constraints on tree-level
couplings of $Z^{\prime}$ to neutrinos, charged leptons, and quarks are quite
stringent for $M_{Z^{\prime}}\lesssim 1{\rm GeV}$. In the following sections,
we investigate the effects of the non-standard neutrino interactions (NSI) or
$Z^{\prime}$ interactions on monojet production at the LHC and NC scattering
at FASER$\nu$, which signify a very high energy and a medium energy scale,
respectively.

The pursuit of NSI’s is one of the main goals of current and future neutrino
experiments. The NSI’s can be categorized into charged current and neutral
current ones. Specifically, we are looking at neutral current
NSI’sProceedings:2019qno

${\cal
L}_{NC}=-2\sqrt{2}G_{F}\,\sum_{f,P,\alpha,\beta}\,\epsilon_{\alpha\beta}^{f,P}\,\left(\bar{\nu}_{\alpha}\gamma^{\mu}P_{L}\nu_{\beta}\right)\,\left(\bar{f}\gamma_{\mu}Pf\right)\;,$
(2)

where $G_{F}$ is the Fermi constant, $\alpha,\beta$ are flavor indices,
$(f,f^{\prime})=(d,u)$, $P=P_{L}$ or $P_{R}$ is the chirality projection
operator. The parameters $\epsilon_{\alpha\beta}^{f,P}$ quantify the strength
of the NC NSI’s, $\alpha,\beta=e,\mu,\tau$, $f=u,d$. For simplicity we only
consider the flavor-conserving interactions on the quark leg, while the
neutrino leg allows for changes in neutrino flavors. Note that the neutrino
field $\nu_{L}$ originates from the lepton doublet $L$, such that the above
interactions can be generated from SM gauge invariant higher dimensional
operators, such as

$-\frac{1}{\Lambda^{2}}\,\left(\bar{L}_{\alpha}\gamma_{\mu}L_{\beta}\right)\,\left[\bar{Q}\gamma^{\mu}P_{L}Q+\bar{u}_{R}\gamma^{\mu}P_{R}u_{R}+\bar{d}_{R}\gamma^{\mu}P_{R}d_{R}\right]$
(3)

where $L$ is the lepton doublet, $Q$ is the quark doublet, $u_{R},d_{R}$ are
the quark singlets. We can then equate to obtain

$\epsilon^{f,P}_{\alpha\beta}=\frac{1}{2\sqrt{2}G_{F}\Lambda^{2}}\;.$ (4)

Since in this work we deal with the effects of the NC NSI’s at the FASER$\nu$
and the LHC, one may concern about the validity of the effective operators in
Eq. (2). The simplified $Z^{\prime}$ model in Eq. (1) converges back to Eq.
(2) when $m_{Z^{\prime}}\to\infty$.

Similarly, when the square of momentum transfer $\hat{s},|\hat{t}|$ are much
smaller than $M_{Z^{\prime}}$, the ratio $\epsilon^{f,P}_{\alpha\beta}$ can be
approximated by $\epsilon_{eff}$,

$\epsilon_{eff}=\frac{g_{q}g_{\nu}}{2\sqrt{2}G_{F}M^{2}_{Z^{\prime}}}\;.$ (5)

Straightly speaking here the $\epsilon_{eff}$ is not the same as
$\epsilon^{f,P}_{\alpha\beta}$ of Eq. (2), but for ease of comparison to those
limits obtained at low energies. In this work, we first work out the
sensitivity constraints in terms of coupling strengths of $Z^{\prime}$, and
then later translate back to the effective coupling $\epsilon_{eff}$’s.

##  III Effects of $Z^{\prime}$ on Monojet Production

A number of new physics models, such as large extra dimensions, invisibly
decaying scalar bosons, sterile neutrinos and dark matter models, can give
rise to missing-energy signals at the LHC, other than the active neutrinos.
The visible object in such events would be the single jet radiating off the
initial quark legs, giving rise to monojet events plus large missing energy.
In the current $Z^{\prime}$ model, the $Z^{\prime}$ boson can be produced
associated with a jet, followed by the $Z^{\prime}$ decay into neutrinos.
Thus, the signature is a single jet plus large missing energy. In the
following, we calculate the production rates of monojet production due to the
$Z^{\prime}$ interactions. Without loss of generality we assume the
$Z^{\prime}$ boson couplings to $u$ and $d$ are the same, and do not couple to
other generations. We can easily extend to different $Z^{\prime}$ couplings in
expense of more independent parameters. After computing the production rates
for monojet events, we can then use a recent experimental result on monojet
production ATLAS:2021kxv  to put bounds on the product of couplings
$(g_{q}g_{\nu})$. Note that production of monojet events has been studied to
test effective neutrino-quark interactions Friedland:2011za ; Pandey:2019apj ;
Babu:2020nna ; Liu:2020emq . An improvement on the constraints can be achieved
from previous works because we have used the most recent result on monojet
production ATLAS:2021kxv .

Figure 1:  The ratio of the cross sections $\sigma_{\rm Tot}/\sigma_{\rm SM}$
versus the mass of the $Z^{\prime}$ boson, in which we have used
$g_{q}=g_{\nu}=1$.

The process that we calculate is

$pp\to\bar{\nu}_{\beta}\nu_{\beta}+j,$

where we sum over all three neutrino flavors assuming their couplings
$g_{\nu}$ to be the same, and $j$ refers to either $q,\bar{q},g$. The
contributing Feynman diagrams include the SM $Z$ boson and the $Z^{\prime}$
boson exchanged in $s$-channel. When the mass $M_{Z^{\prime}}\to\infty$ the SM
result is restored. In principle, the $Z$ and $Z^{\prime}$ diagrams interfere
with each other, such that the interference term is proportional to the
couplings $(g_{q}g_{\nu})$ while the sole $Z^{\prime}$ contribution is
proportional to $(g_{q}g_{\nu})^{2}$.

In the calculation, we generate the aforementioned process using
MadGraph5aMC@NLO Alwall:2011uj ; Alwall:2014hca  with the model file generated
by the effective Lagrangian in Eq. (1), followed by parton showering and
hadronization with PYTHIA8 Sjostrand:2014zea , detector simulations carried
out by Delphes3 package deFavereau:2013fsa . The total cross-section
$\sigma_{\rm Tot}$ for $pp\to\nu\bar{\nu}+1j$ can be expressed as follows

$\sigma_{\rm Tot}=\sigma_{Z^{\prime}}+\sigma_{\rm Int}+\sigma_{\rm SM}\;,$ (6)

where $\sigma_{Z^{\prime}}$ is the cross-section of the aforementioned process
only with the $Z^{\prime}$ propagator, $\sigma_{\rm Int}$ is the interference
term and $\sigma_{\rm SM}$ is the standard model cross-section. We show the
ratio of $\sigma_{\rm Tot}/\sigma_{\rm SM}$ in Fig. 1. It is clear that the
$\sigma_{\rm Tot}$ approaches $\sigma_{\rm SM}$ as $M_{Z^{\prime}}$ becomes
very large. Note that the total decay width of the $Z^{\prime}$ boson is
assumed to be $\frac{\Gamma_{Z^{\prime}}}{M_{Z^{\prime}}}=0.1$.

###  III.1 Sensitivity reach on parameter space of the $Z’$ model and NSI’s

Here we derive the bounds on the product of the $Z^{\prime}$ couplings
$(g_{q}g_{\nu})$ as a function of $M_{Z^{\prime}}$ based on a recent result on
monojet production by the ATLAS experiment ATLAS:2021kxv . Later, our goal is
to translate such constraints into the conventional NSI parameters
$\epsilon_{eff}=\epsilon_{u}=\epsilon_{d}$ defined in Eq. (2).

We follow closely the experimental cuts outlined in the ATLAS paper
ATLAS:2021kxv  in order to directly use their upper limits on the monojet
production cross sections. Their results were based on the monojet search at
13 TeV with an integrated luminosity of 139 $\rm fb^{-1}$ ATLAS:2021kxv .
Events are selected with $E_{T}^{\rm miss}>200\,{\rm GeV}$, a leading jet with
$p_{T}>150\,{\rm GeV}$ and $|\eta|<2.4$ and upto three jet with
$p_{T}>30\,{\rm GeV}$ and $|\eta|<2.8$, as well as additional cuts specified
in ATLAS:2021kxv . Jets are defined with the anti-$k_{t}$ jet algorithm with a
cone size $R=0.4$.

With all the acceptance cuts the same as Ref. ATLAS:2021kxv , we still need
the overall efficiency in order to obtain the event rates to compare with the
experimental results. We rely on an information given in Ref. ATLAS:2021kxv .
A signal model with an axial-vector gauge boson $Z_{A}$, via which a pair of
dark matter particles $\chi$ can be produced in $s$-channel, was investigated.
The reported ”acceptance $\times$ efficiency” in the kinematic region EM0
111Here $p_{T}^{\rm recoil}$ is the same as $E_{T}^{\rm miss}$ for signal
models. ($p_{T}^{\rm recoil}=200-250$ GeV) was 13%. Since the event topology
of such a signal ($pp\to Z_{A}+j\to\bar{\chi}\chi+j$) is similar to our signal
($pp\to Z^{\prime}+j\to\bar{\nu}\nu+j$), we can then compare our acceptance to
their value of “acceptance $\times$ efficiency”. Therefore, we obtain an
efficiency of $0.582$, which is then applied to all our event rates. We have
calculated bounds of $\sqrt{g_{q}g_{\nu}}$ using the 95% C.L. upper limits on
the signal event rates in a number of kinematic regions defined in Ref.
ATLAS:2021kxv  (see Table 9 of Ref. ATLAS:2021kxv ). The resulting limits are
within a factor of two among one another. We show in the left panel of Fig. 2
the bounds on $\sqrt{g_{q}g_{\nu}}$ based on the 95% C.L. upper limit on the
observed event rate $S^{95}_{\rm obs}=11937$ in the kinematic region IM3
($p_{T}^{\rm recoil}>350$ GeV) ATLAS:2021kxv .

Figure 2:  Left: Monojet bounds on the product of couplings
$\sqrt{g_{q}g_{\nu}}$ versus the $Z^{\prime}$ mass. Right: Constraints on
$\epsilon_{eff}$ of the NSI’s translated from $(g_{q}g_{\nu})$ using Eq. (5)
versus $M_{Z^{\prime}}$. Here we have assumed
$\epsilon_{ee}=\epsilon_{\mu\mu}=\epsilon_{\tau\tau}=\epsilon_{eff}$.

Next we can translate the bounds on $\sqrt{g_{q}g_{\nu}}$ to $\epsilon_{eff}$
using Eq. (5). The bounds on $\epsilon_{eff}$ are shown on the right panel of
Fig. 2. We can see that the best limit on $\epsilon_{eff}$ appears around
$M_{Z^{\prime}}\sim 2$ TeV. The bound becomes less stringent as the
$Z^{\prime}$ mass increases, because the $Z^{\prime}$ becomes more difficult
to be produced directly.

##  IV Effects of $Z^{\prime}$ and Neutral-current NSI’s interactions at
FASER$\nu$

FASER FASER:2019dxq ; FASER:2019aik  is an approved experiment located about
480 m away from the interaction point (IP) of the ATLAS detector down along
the direction of the proton beam. It is well known that huge number of
hadrons, such as pions, kaons and other hadrons, are produced along the beam
direction. These hadrons will decay during the flight, thus producing a lot of
neutrinos of all three flavors at very high energy up to a few TeV.

There is a proposed new component, called FASER$\nu$ FASER:2019dxq , to be put
in front of the FASER detector. It is an $\rm 25cm\times 25cm\times 1.5m$
emulsion detector, consisting of 1000 layers of emulsion films interleaved
with 1-mm-thick tungsten plates with mass 1.2 tons. The main goal of
FASER$\nu$ is to distinguish various flavors of neutrinos. Indeed, it can
measure the flux of electron, muon, and tau neutrinos coming off from the IP
of the ATLAS detector, which can be done by detecting the charged lepton
coming off the charged-current (CC) scattering. Notably, muon nutrino is the
most abundant due to production of charged pions and kaons while tau neutrino
is the least as it requires at least the heavier mesons. like $D_{s}$ meson,
in order to decay into $\tau\nu_{\tau}$. On the other hand, it is also
feasible to measure the neutral-current (NC) scattering of the neutrinos
FASER:2019dxq  making use of the emulsion detector, although the detection of
NC interactions is somewhat more difficult than the CC one. The total cross-
section ($\sigma_{\nu N}$) of the NC scattering at the FASER$\nu$ detector can
be expressed as $\sigma_{\nu N}=n\sigma_{\nu n}+p\sigma_{\nu p}$, where
$\sigma_{\nu n}$ and $\sigma_{\nu p}$ are the neutrino-neutron and neutrino-
proton scattering cross sections, and $n$ and $p$ are the number of neutrons
and protons in the tungsten atom, respectively.

In this section we compute the sensitivity of FASER$\nu$ to the NC NSI’s due
to physics beyond the SM. Similar to the last section, we use the same
simplifed $Z^{\prime}$ model (see Eq. (1) ) to calculate the sensitivity reach
at FASER$\nu$. The effect of $Z^{\prime}$ is similar to that at the LHC, other
than the fact that the $Z$ and $Z^{\prime}$ bosons are exchanged in
$t$-channel in the NC deep-inelastic scattering, such that the most
significant effect of $Z^{\prime}$ appears in the small $M_{Z^{\prime}}$
region. The SM result is restored as $M_{Z^{\prime}}\to\infty$. We estimate
the 95% C.L. sensitivity reach on the parameter space of the $Z^{\prime}$
model. We show that the best sensitivity can be achieved in the small
$M_{Z^{\prime}}$ region that it is highly complementary to that obtained by
monojet production at the LHC.

###  IV.1 $Z^{\prime}$ Interactions at FASER$\nu$

The square of the Feynman amplitude for the subprocess
$\nu(p_{1})q(p_{2})\to\nu(k_{1})q(k_{2})$, where $q=u,d$ and the 4-momenta of
each particle is shown in parenthesis, is given by

$\sum\left|{\cal M}\right|^{2}=4\hat{u}^{2}|M^{\nu
q}_{LL}|^{2}+4\hat{s}^{2}|M^{\nu q}_{LR}|^{2}$ (7)

where the reduced amplitudes $M^{\nu q}_{L\beta}$ are given by

$M^{\nu
q}_{L\beta}\left(\nu(p_{1})q(p_{2})\to\nu(k_{1})q(k_{2})\right)=\frac{e^{2}g_{Z}^{\nu}g_{Z}^{q_{\beta}}}{\sin^{2}\theta_{w}\cos^{2}\theta_{w}}\,\frac{1}{\hat{t}-M_{Z}^{2}}+\frac{g_{\nu}g_{q_{\beta}}}{\hat{t}-M_{Z^{\prime}}^{2}}\;,$
(8)

where $\beta=L,R$, $g_{Z}^{f_{L}}=T_{3f}-Q_{f}\sin^{2}\theta_{w}$, and
$g_{Z}^{f_{R}}=-Q_{f}\sin^{2}\theta_{w}$ are the SM $Z$ couplings to the
fermion $f_{L}$ and $f_{R}$, and $\theta_{w}$ is the Weinberg angle. Here
$\hat{s},\hat{t},\hat{u}$ are the usual Mandelstam variables. In our
$Z^{\prime}$ model, the couplings $g_{q_{L}}$ and $g_{q_{R}}$ are the same. It
is easy to see that when $M_{Z^{\prime}}\to\infty$ the SM result is restored.
We used MadGraph5aMC@NLO Alwall:2011uj ; Alwall:2014hca  for fixed target
deep-inelastic neutrino-nucleon scattering computation. We build the model
file for Eq (1) using Feynrules Alloul:2013bka . We consider the $Z^{\prime}$
mass ranging from 0.01 GeV to 10 TeV and show the scattering cross section
normalized by $E_{\nu}$ in Fig. 3. The process cross-section $\sigma_{\nu{\rm
N}}$ decreases with $M_{Z^{\prime}}$. At heavy $M_{Z^{\prime}}$ mass regime
the total cross section approaches to the standard model values, which was
already reported in FASER:2019dxq ; IceCube:2017roe ; Ismail:2020yqc ; Cooper-
Sarkar:2011jtt . We used the values of neutrino flux FASER:2019dxq  for the
evaluation of neutrino-nucleus interaction with FASER$\nu$

Figure 3:  Left: Deep-inelastic neutral-current scattering cross section
normalized by the energy $E_{\nu}$ of the incoming neutrino beam. Right: the
same as the left but with anti-neutrino beam. Here N is the tungsten nucleus.
We have used the CTEQ6L1Pumplin:2002vw  for parton distribution functions. We
have set $g_{q}=g_{\nu}=0.1$.  Figure 4:  The energy spectrum of neutrinos
with NC interactions mediated with $Z$ and $Z^{\prime}$ in a 1-ton tungsten
detector with dimensions $\rm 25\,cm\times 25\,cm\times 1\,m$ centered on the
beam collision axis at the FASER location at the 14 TeV LHC with 150 $\rm
fb^{-1}$

Next the energy spectra for the neutral-current interactions of three flavors
of neutrinos are shown in Fig. 4 for a number of values for $M_{Z^{\prime}}$.
The expected number of NC events of three flavors of neutrinos versus
$E_{\nu}>10\rm\leavevmode\nobreak\ GeV$ for various $Z^{\prime}$ mass in
FASER$\nu$ can be obtained from the corresponding energy spectrum. The highest
number of NC events was reported in the $\nu_{\mu}$ channel, while the lowest
number of NC events in the $\nu_{\tau}$ channel. In Fig. 4, we sum up the
contributions from both neutrino and anti-neutrino events. Here we have
assumed a benchmark detector made of tungsten with dimensions $\rm
25\,cm\times 25\,cm\times 1\,m$ at the 14 TeV LHC with an integrated
luminosity of $L=\rm 150\,fb^{-1}$. We use the neutrino fluxes and energy
spectra obtained in FASER:2019dxq  to study the neutrinos that pass through
FASER$\nu$. We find that muon neutrinos are mostly produced from charged-pion
decays, electron neutrinos from hyperon, kaon, and $D$-meson decays, and tau
neutrinos from $D_{s}$ meson decays. With average energies ranging from 600
GeV to 1 TeV, the spectra of the three neutrino flavors cover a broad energy
range.

To estimate the sensitivity reach in the parameter space $(g_{\nu}q_{q})$ of
the $Z^{\prime}$ model, we first calculate the predicted number of events
$N_{\rm BSM}$ for the $Z^{\prime}$ model and the SM number of events $N_{\rm
SM}$, and treat the statistical error as $\sqrt{N_{\rm BSM}}$ and systematic
uncertainty $\sigma_{\rm norm}$ as a fraction ($\sigma_{\rm norm}=20\%,5\%$)
of the normalization of the SM predictions. We then define the measure of
$\chi^{2}$ as a function of $(g_{\nu}q_{q})$ and a nuisance parameter $\alpha$
as follows Davidson:2003ha :

$\displaystyle\chi^{2}(g_{q}g_{\nu},\alpha)$ $\displaystyle=$
$\displaystyle\min_{\alpha}\Biggr{[}\frac{(N^{\nu_{e}}_{BSM}-(1+\alpha)N^{\nu_{e}}_{SM})^{2}}{N^{\nu_{e}}_{BSM}}+\frac{(N^{\nu_{\mu}}_{BSM}-(1+\alpha)N^{\nu_{\mu}}_{SM})^{2}}{N^{\nu_{\mu}}_{BSM}}$
(9)
$\displaystyle+\frac{(N^{\nu_{\tau}}_{BSM}-(1+\alpha)N^{\nu_{\tau}}_{SM})^{2}}{N^{\nu_{\tau}}_{BSM}}+\left(\frac{\alpha}{\sigma_{norm}}\right)^{2}\Biggr{]}\;,$

Figure 5:  Left: Sensitivity reach on the product of couplings
$\sqrt{g_{q}g_{\nu}}$ versus the $Z^{\prime}$ mass achieved at FASER$\nu$.
Right: Sensitivity reach in terms of $\epsilon_{eff}$ of the NSI’s translated
from $(g_{q}g_{\nu})$ using Eq. (5). Systematic uncertainty $\sigma_{\rm
norm}=5,20\%$ and without systematic uncertainties are shown.

where $N_{\rm BSM}=N_{Z^{\prime}}+N_{\rm int}+N_{\rm SM}$ and the minimization
is over the nuisance parameter $\alpha$. Here $N_{Z^{\prime}}$ is the number
of events from the $Z^{\prime}$ diagram only, $N_{\rm int}$ is the
interference term. Here we have treated the systematic uncertainties in each
neutrino flavor to be the same and use only one nuisance parameter $\alpha$.
Physics-wise the systematic uncertainties come from theoretical calculations,
the flux of neutrinos from the ATLAS IP, detector response, etc. We show in
Fig. 5 (Left) the 95% C.L. sensitivity reach (corresponding to
$\chi^{2}=3.84$) of the product $\sqrt{g_{q}g_{\nu}}$ versus $M_{Z^{\prime}}$
at FASER$\nu$. The higher the systematic uncertainty the weaker the limit on
$\sqrt{g_{q}g_{\nu}}$ will be. Nevertheless, the differences among
$\sigma_{\rm norm}=5\%,20\%$ and without systematic uncertainties are
relatively small. The sensitivity reach on $\sqrt{g_{q}g_{\nu}}$ is the best
at very small $M_{Z^{\prime}}$ around $10^{-4}$ at $M_{Z^{\prime}}=0.01$ GeV
and reduces to about $1$ at $M_{Z^{\prime}}=1000$ GeV. Now we can translate
the bounds on $\sqrt{g_{q}g_{\nu}}$ to $\epsilon_{eff}$ using Eq. (5). The
bounds on $\epsilon_{eff}$ are shown on the right panel of Fig. 5. We could
see the best limit of $\epsilon_{eff}$ occurs at $M_{Z^{\prime}}\sim 100$ GeV
irrespective of the choice of $\sigma_{\rm norm}$. The curve without
systematic uncertainties is giving the best limit of $\epsilon_{eff}$ in the
whole $M_{Z^{\prime}}$ space. The limit on $\epsilon_{eff}$ is clearly getting
stronger as $M_{Z^{\prime}}$ increases from 0.01 to 100 GeV, but staying flat
after $M_{Z^{\prime}}=$100 GeV onward. The monojet study also shows similar
behavior of $\epsilon_{eff}$ at the higher $M_{Z^{\prime}}$ region.

Figure 6:  Sensitivity reach on the product of couplings
$\sqrt{g_{q}g_{\nu_{\beta}}}$ for each neutrino flavor $\beta=e,\mu,\tau$
versus the $Z^{\prime}$ mass at FASER$\nu$ (without systematic uncertainties).

FASER$\nu$ is primarily designed for the purpose of identifying the flavors of
neutrinos Bahraminasr:2020ssz ; Bakhti:2020szu . The expected sensitivity for
each flavor at FASER$\nu$ is shown in Fig. 6. We only show the curves with no
systematic uncertainty included. The curves can be compared to the
corresponding one “Without Systematic” of Fig. 5. For the evaluation of
$\chi^{2}$, we consider the special case of Eq.(9) with $\alpha=0$ and it
reduces to

$\chi^{2}(g_{q}g_{\nu_{\beta}})=[\frac{(N^{\nu_{\beta}}_{BSM}-N^{\nu_{\beta}}_{SM})^{2}}{N^{\nu_{\beta}}_{BSM}}]$
(10)

where $\beta=e,\mu,\tau$. The green and orange curves of Fig. 6 depict the
sensitivity reach of $\sqrt{g_{q}g_{\nu_{\mu}}}$ and $\sqrt{g_{q}g_{\nu_{e}}}$
($\chi^{2}=3.84$) versus $M_{Z^{\prime}}$, respectively, while the blue curve
represent the sensitivity reach of $\sqrt{g_{q}g_{\nu_{\tau}}}$ versus
$M_{Z^{\prime}}$.

It is clear from the Fig. 5 and Fig. 6 that overall sensitivity reach of
$\sqrt{g_{q}g_{\nu}}$ is dominated by $\sqrt{g_{q}g_{\nu_{\mu}}}$. For each
neutrino flavor the sensitivity reach on $\sqrt{g_{q}g_{\nu_{\beta}}}$ hitting
unity at different $Z^{\prime}$ mass: (i) $\sqrt{g_{q}g_{\nu_{\mu}}}$ (orange
curve) approaches to 1 at $M_{Z^{\prime}}\sim$ 1000 GeV, (ii)
$\sqrt{g_{q}g_{\nu_{e}}}$ (blue curve) approaches to 1 at $M_{Z^{\prime}}\sim$
800 GeV, and (iii) $\sqrt{g_{q}g_{\nu_{\tau}}}$ (green curve) reaches the
unity faster than the other two flavors at $M_{Z^{\prime}}\sim 500$ GeV.

##  V Complementarity of Monojet and FASER$\nu$ Results

Monojet production at the LHC and the NC deep-inelastic scattering at
FASER$\nu$ cover different energy scales. It would be useful to put both
results together. We show in Fig. 7 the future sensitivity reach at FASER$\nu$
and the most updated constraints due to monojet production at the LHC. It is
interesting to see that FASER$\nu$ is mostly sensitive to small
$M_{Z^{\prime}}$ region from $10^{-2}-O(100)$ GeV while monojet production is
more sensitive for $M_{Z^{\prime}}\gtrsim 100$ GeV to a few TeV. In Fig. 7, we
also include other existing constraints at 95% C.L., including (i) the CCFR
measurement of the neutrino trident cross-section Ballett:2019xoj , (ii) the
search of SM $Z$ boson decay to 4 charged leptons in CMS Altmannshofer:2014pba
and ATLAS Altmannshofer:2016jzy ; ATLAS:2014jlg  reinterpreted under the
hypothesis of $Z\rightarrow Z^{\prime}\mu\mu$, (iii) the search of
$e^{+}e^{-}\rightarrow\mu^{+}\mu^{-}Z^{\prime},$ followed by
$Z^{\prime}\rightarrow\mu^{+}\mu^{-}$ from BaBarCMS:2018yxg , (iv) bounds from
Borexino BaBar:2016sci ; Kamada:2015era , (v) $(g-2)_{\mu}\leavevmode\nobreak\
2\sigma$ band related to the anomalous magnetic moment of muon
Gninenko:2020xys , (vi) the constraint from the present COHERENT data
Cadeddu:2020nbr ; Denton:2018xmq ; Liao:2017uzy , (vii) the LMA-DARK solution
Cadeddu:2020nbr  with $x=0$ (with $x=2$), and (viii) the LEP II bounds on
couplings to electrons derived from Buckley:2011vc ; Carena:2004xs , where we
have assumed a single fermion helicity in the $Z^{\prime}$ coupling. The
constraints on the couplings of the $Z^{\prime}$ to leptons are significantly
more stringent than those to quarks. In particular, the process
$e^{+}e^{-}\rightarrow Z^{\prime}\rightarrow e^{+}e^{-}$ leads to a constraint
of $g_{ee}^{Z^{\prime}}\leq 0.044\times(M_{Z^{\prime}}/\rm 200GeV)$ for
$Z^{\prime}$ masses above roughly 200 GeV.

In the intermediate mass range ($1\,{\rm GeV}\lesssim M_{Z^{\prime}}<50\,{\rm
GeV}$), the FASER$\nu$’s sensitivities are comparable with the existing
constraints, except for the range $M_{Z^{\prime}}=5-50$ GeV, where the CMS and
ATLAS searches on SM $Z$ boson decay into 4 charged leptons are somewhat
better. In the low mass regime ($M_{Z^{\prime}}=0.01-1\rm GeV)$ the COHERENT
results are better than the FASER$\nu$ sensitivites. The LMA-DARK solution is
also better than FASER$\nu$ sensitivites in $0.01\,\rm GeV<M_{Z^{\prime}}\leq
0.1\,\rm GeV$ region, however in the higher $M_{Z^{\prime}}$ region
($M_{Z^{\prime}}>0.1\,\rm GeV$) FASER$\nu$ can constrain better than the LMA-
DARK. In the high mass regime (100 GeV $\lesssim M_{Z^{\prime}}$), the LHC
Monojet results constrain better than the sensitivities offered by the
FASER$\nu$, wheree we can see the crossover between FASER$\nu$ and LHC-monojet
results at $M_{Z^{\prime}}\sim$250 GeV.

Here we make a brief comparison with the sensitivity achieved at the DUNE
near-detector. The $\nu-e$ scattering sensitivity to the
$L_{e}-L_{\mu}\leavevmode\nobreak\ Z^{\prime}$ model at 90% C.L. was performed
in Ref. He:1991qd , and the dimuon neutrino trident sensitivity to the
$L_{\mu}-L_{\tau}$ model He:1991qd  with no kinetic mixing at 90% C.L. were
reported in Ref. Ballett:2019xoj . Sensitivity on $g^{\prime}$ with
$L_{e}-L_{\mu}\leavevmode\nobreak\ Z^{\prime}$ model reaches the best at very
small $M_{Z^{\prime}}$ around $\sim 5\times 10^{-5}$ at $M_{Z^{\prime}}$=0.01
GeV and rises to about 0.01 when $M_{Z^{\prime}}$=10 GeV. For the case of
$L_{\mu}-L_{\tau}$ model the $g^{\prime}$ value with $M_{Z^{\prime}}$=0.01 is
in the order of $\sim 2\times 10^{-4}$ and rises to $\sim$ 0.01 at
$M_{Z^{\prime}}$=10 GeV. On the other hand, the best FASER$\nu$ sensitivity
that we can achieve is $\sqrt{g_{q}g_{\nu}}\sim 10^{-4}$ at
$M_{Z^{\prime}}=0.01$ GeV and rises to about $0.01$ at $M_{Z^{\prime}}=10$
GeV. Therefore, we can see that the FASER$\nu$ sensitivity is comparable to
that of DUNE.

Figure 7:  Future sensitivity reach at FASER$\nu$ and the most updated
constraint due to monojet production at the LHC at 95% C.L. Other existing
constraints at 95% C.L. shown include (i) the CCFR measurement (red area) of
the neutrino trident cross-section Ballett:2019xoj , (ii) the search of SM $Z$
boson decay to 4 charged leptons in CMS Altmannshofer:2014pba  (gray area) and
ATLAS Altmannshofer:2016jzy ; ATLAS:2014jlg  (yellow area) reinterpreted under
the hypothesis of $Z\rightarrow Z^{\prime}\mu\mu$, (iii) the search of
$e^{+}e^{-}\rightarrow\mu^{+}\mu^{-}Z^{\prime},\;Z^{\prime}\rightarrow\mu^{+}\mu^{-}$
from BaBar CMS:2018yxg  (purple area), (iv) bounds from Borexino BaBar:2016sci
; Kamada:2015era  (green area), (v) $(g-2)_{\mu}\leavevmode\nobreak\ 2\sigma$
band related to the anomalous magnetic moment of muon Gninenko:2020xys  (cyan
area), (vi) the constraint from the present COHERENT data Cadeddu:2020nbr ;
Denton:2018xmq ; Liao:2017uzy  (blue area), (vii) the LMA-DARK solution
Cadeddu:2020nbr  with $x=0$ (with $x=2$) (black curve), and (viii) the
constrain from LEP II Buckley:2011vc ; Carena:2004xs  (orange area).

##  VI Conclusions

In this paper, we have studied the neutral-current scattering between
neutrinos and nuclei in the FASER$\nu$ detector and calculated the expected
sensitivity reach on possible NSI’s using a simplified $Z^{\prime}$ model. We
investigated the advantage of FASER$\nu$ in wide mass range search for
$Z^{\prime}$ and to determine the flavor dependence of the coupling between
neutrino and this new boson, for which we found that FASER$\nu$ is sensitive
to $g_{\nu_{\mu}}$ because of the larger statistics. We also found that the
impact of systematical uncertainty due to normalization is relatively small in
the smaller $M_{Z^{\prime}}$ region.

We have also investigated the effects of the simplified $Z^{\prime}$ model on
monojet production at the LHC, followed by an update on the existing bound
using the most recent results on monojet production at the LHC with 139 fb-1
luminosity. We have found substantial improvement over previous works.

While the FASER$\nu$ can achieve the best sensitivity at small
$M_{Z^{\prime}}$ regime, the sensitivity using monojet production, on the
other hand, is more profound at high mass region. Thus, complementarity in
mass range coverage is established. Overall, the FASER$\nu$ offers a
sensitivity reach better than the existing constraints at low mass region
($M_{Z^{\prime}}<0.1$) GeV, except for the COHERENT constraint and for the
DUNE near-detector $\nu-e$ scattering sensitivity. The FASER$\nu$ sensitivity
is comparable to existing constraints in the intermediate mass region
($0.1\lesssim M_{Z^{\prime}}\lesssim 10$ GeV). We explored the capability of
thee FASER$\nu$ detector to discern individual neutrino flavors. In both the
FASER$\nu$ experiment and LHC we obtained the best limit for $\epsilon_{eff}$
based on the translation of $\sqrt{g_{q}g_{\nu}}$. Further full detector
simulation at FASER$\nu$ is called for establishing the feasibilty.

More and more particle-physics experiments or cosmological observatories
provide the bounds at the lower mass region of $Z^{\prime}$. However, we have
not seen any signal so far. One may be more interested in the heavier
$Z^{\prime}$ models, for which FASER$\nu$ and monojet play an important role
in that search. We are looking forward to the upgrade of FASER/FASER$\nu$,
which is being discussed in the collaboration group. For a complete picture of
$Z^{\prime}$ search, our suggestion is to cover the mass range around $100$
GeV and those above $\sim 300$ GeV, which are still lack of constraints.

## Acknowledgement

Special thanks to Felix Kling, Shih-Chieh Hsu and Zeren Simon Wang for
enlightening discussion. Also thanks to Olivier Mattelaer for a wonderful
usage of Madgraph. TC acknowledges the support from National Center for
Theoretical Sciences. The work was supported in part by Taiwan MoST with grant
no. MOST-110-2112-M-007-017-MY3.

## References

  * (1) N. Blinov, K. J. Kelly, G. Z. Krnjaic and S. D. McDermott, “Constraining the Self-Interacting Neutrino Interpretation of the Hubble Tension,” Phys. Rev. Lett. 123, no.19, 191102 (2019) doi:10.1103/PhysRevLett.123.191102 [arXiv:1905.02727 [astro-ph.CO]].
  * (2) F. Y. Cyr-Racine and K. Sigurdson, “Limits on Neutrino-Neutrino Scattering in the Early Universe,” Phys. Rev. D 90, no.12, 123533 (2014) doi:10.1103/PhysRevD.90.123533 [arXiv:1306.1536 [astro-ph.CO]].
  * (3) C. D. Kreisch, F. Y. Cyr-Racine and O. Doré, “Neutrino puzzle: Anomalies, interactions, and cosmological tensions,” Phys. Rev. D 101, no.12, 123505 (2020) doi:10.1103/PhysRevD.101.123505 [arXiv:1902.00534 [astro-ph.CO]].
  * (4) S. Ghosh, R. Khatri and T. S. Roy, “Can dark neutrino interactions phase out the Hubble tension?,” Phys. Rev. D 102, no.12, 123544 (2020) doi:10.1103/PhysRevD.102.123544 [arXiv:1908.09843 [hep-ph]].
  * (5) G. De Lellis, “Search for Hidden Particles (SHiP): a new experiment proposal,” Nucl. Part. Phys. Proc. 263-264, 71-76 (2015) doi:10.1016/j.nuclphysbps.2015.04.013
  * (6) M. Anelli et al. [SHiP], “A facility to Search for Hidden Particles (SHiP) at the CERN SPS,” [arXiv:1504.04956 [physics.ins-det]].
  * (7) H. Abreu et al. [FASER], “Detecting and Studying High-Energy Collider Neutrinos with FASER at the LHC,” Eur. Phys. J. C 80, no.1, 61 (2020) doi:10.1140/epjc/s10052-020-7631-5 [arXiv:1908.02310 [hep-ex]].
  * (8) A. Ariga et al. [FASER], “FASER: ForwArd Search ExpeRiment at the LHC,” [arXiv:1901.04468 [hep-ex]].
  * (9) X. G. He, G. C. Joshi, H. Lew and R. R. Volkas, “Simplest Z-prime model,” Phys. Rev. D 44, 2118-2132 (1991) doi:10.1103/PhysRevD.44.2118
  * (10) W. Shepherd, T. M. P. Tait and G. Zaharijas, “Bound states of weakly interacting dark matter,” Phys. Rev. D 79, 055022 (2009) doi:10.1103/PhysRevD.79.055022 [arXiv:0901.2125 [hep-ph]].
  * (11) G. Aad et al. [ATLAS], “Search for new phenomena in events with an energetic jet and missing transverse momentum in $pp$ collisions at $\sqrt{s}$ =13 TeV with the ATLAS detector,” Phys. Rev. D 103, no.11, 112006 (2021) doi:10.1103/PhysRevD.103.112006 [arXiv:2102.10874 [hep-ex]].
  * (12) G. y. Huang, T. Ohlsson and S. Zhou, “Observational Constraints on Secret Neutrino Interactions from Big Bang Nucleosynthesis,” Phys. Rev. D 97, no.7, 075009 (2018) doi:10.1103/PhysRevD.97.075009 [arXiv:1712.04792 [hep-ph]].
  * (13) J. B. Dent, F. Ferrer and L. M. Krauss, “Constraints on Light Hidden Sector Gauge Bosons from Supernova Cooling,” [arXiv:1201.2683 [astro-ph.CO]].
  * (14) R. Harnik, J. Kopp and P. A. N. Machado, “Exploring nu Signals in Dark Matter Detectors,” JCAP 07, 026 (2012) doi:10.1088/1475-7516/2012/07/026 [arXiv:1202.6073 [hep-ph]].
  * (15) A. E. Nelson and N. Tetradis, “CONSTRAINTS ON A NEW VECTOR BOSON COUPLED TO BARYONS,” Phys. Lett. B 221, 80-84 (1989) doi:10.1016/0370-2693(89)90196-2
  * (16) S. Tulin, “New weakly-coupled forces hidden in low-energy QCD,” Phys. Rev. D 89, no.11, 114008 (2014) doi:10.1103/PhysRevD.89.114008 [arXiv:1404.4370 [hep-ph]].
  * (17) B. Aubert et al. [BaBar], “Search for Invisible Decays of a Light Scalar in Radiative Transitions $\upsilon_{3S}\to\gamma$ A0,” [arXiv:0808.0017 [hep-ex]].
  * (18) R. Essig, J. Mardon, M. Papucci, T. Volansky and Y. M. Zhong, “Constraining Light Dark Matter with Low-Energy $e^{+}e^{-}$ Colliders,” JHEP 11, 167 (2013) doi:10.1007/JHEP11(2013)167 [arXiv:1309.5084 [hep-ph]].
  * (19) P. S. Bhupal Dev, K. S. Babu, P. B. Denton, P. A. N. Machado, C. A. Argüelles, J. L. Barrow, S. S. Chatterjee, M. C. Chen, A. de Gouvêa and B. Dutta, et al. “Neutrino Non-Standard Interactions: A Status Report,” SciPost Phys. Proc. 2, 001 (2019) doi:10.21468/SciPostPhysProc.2.001 [arXiv:1907.00991 [hep-ph]].
  * (20) A. Friedland, M. L. Graesser, I. M. Shoemaker and L. Vecchi, “Probing Nonstandard Standard Model Backgrounds with LHC Monojets,” Phys. Lett. B 714, 267-275 (2012) doi:10.1016/j.physletb.2012.06.078 [arXiv:1111.5331 [hep-ph]].
  * (21) S. Pandey, S. Karmakar and S. Rakshit, “Strong constraints on non-standard neutrino interactions: LHC vs. IceCube,” JHEP 11, 046 (2019) doi:10.1007/JHEP11(2019)046 [arXiv:1907.07700 [hep-ph]].
  * (22) K. S. Babu, D. Gonçalves, S. Jana and P. A. N. Machado, “Neutrino Non-Standard Interactions: Complementarity Between LHC and Oscillation Experiments,” Phys. Lett. B 815, 136131 (2021) doi:10.1016/j.physletb.2021.136131 [arXiv:2003.03383 [hep-ph]].
  * (23) D. Liu, C. Sun and J. Gao, “Constraints on neutrino non-standard interactions from LHC data with large missing transverse momentum,” JHEP 02, 033 (2021) doi:10.1007/JHEP02(2021)033 [arXiv:2009.06668 [hep-ph]].
  * (24) J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer and T. Stelzer, “MadGraph 5 : Going Beyond,” JHEP 06, 128 (2011) doi:10.1007/JHEP06(2011)128 [arXiv:1106.0522 [hep-ph]].
  * (25) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro, “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations,” JHEP 07, 079 (2014) doi:10.1007/JHEP07(2014)079 [arXiv:1405.0301 [hep-ph]].
  * (26) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen and P. Z. Skands, “An introduction to PYTHIA 8.2,” Comput. Phys. Commun. 191, 159-177 (2015) doi:10.1016/j.cpc.2015.01.024 [arXiv:1410.3012 [hep-ph]].
  * (27) J. de Favereau et al. [DELPHES 3], “DELPHES 3, A modular framework for fast simulation of a generic collider experiment,” JHEP 02, 057 (2014) doi:10.1007/JHEP02(2014)057 [arXiv:1307.6346 [hep-ex]].
  * (28) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, “FeynRules 2.0 - A complete toolbox for tree-level phenomenology,” Comput. Phys. Commun. 185, 2250-2300 (2014) doi:10.1016/j.cpc.2014.04.012 [arXiv:1310.1921 [hep-ph]].
  * (29) M. G. Aartsen et al. [IceCube], “Measurement of the multi-TeV neutrino cross section with IceCube using Earth absorption,” Nature 551, 596-600 (2017) doi:10.1038/nature24459 [arXiv:1711.08119 [hep-ex]].
  * (30) A. Ismail, R. Mammen Abraham and F. Kling, “Neutral current neutrino interactions at FASER$\nu$,” Phys. Rev. D 103, no.5, 056014 (2021) doi:10.1103/PhysRevD.103.056014 [arXiv:2012.10500 [hep-ph]].
  * (31) A. Cooper-Sarkar, P. Mertsch and S. Sarkar, “The high energy neutrino cross-section in the Standard Model and its uncertainty,” JHEP 08, 042 (2011) doi:10.1007/JHEP08(2011)042 [arXiv:1106.3723 [hep-ph]].
  * (32) J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. M. Nadolsky and W. K. Tung, “New generation of parton distributions with uncertainties from global QCD analysis,” JHEP 07, 012 (2002) doi:10.1088/1126-6708/2002/07/012 [arXiv:hep-ph/0201195 [hep-ph]].
  * (33) S. Davidson, C. Pena-Garay, N. Rius and A. Santamaria, “Present and future bounds on nonstandard neutrino interactions,” JHEP 03, 011 (2003) doi:10.1088/1126-6708/2003/03/011 [arXiv:hep-ph/0302093 [hep-ph]].
  * (34) M. Bahraminasr, P. Bakhti and M. Rajaee, “Sensitivities to secret neutrino interaction at FASER$\nu$,” doi:10.1088/1361-6471/ac11c2 [arXiv:2003.09985 [hep-ph]].
  * (35) P. Bakhti, Y. Farzan and S. Pascoli, “Discovery potential of FASER$\nu$ with contained vertex and through-going events,” JHEP 04, 075 (2021) doi:10.1007/JHEP04(2021)075 [arXiv:2010.16312 [hep-ph]].
  * (36) P. Ballett, M. Hostert, S. Pascoli, Y. F. Perez-Gonzalez, Z. Tabrizi and R. Zukanovich Funchal, “$Z^{\prime}$s in neutrino scattering at DUNE,” Phys. Rev. D 100, no.5, 055012 (2019) doi:10.1103/PhysRevD.100.055012 [arXiv:1902.08579 [hep-ph]].
  * (37) W. Altmannshofer, S. Gori, M. Pospelov and I. Yavin, “Neutrino Trident Production: A Powerful Probe of New Physics with Neutrino Beams,” Phys. Rev. Lett. 113, 091801 (2014) doi:10.1103/PhysRevLett.113.091801 [arXiv:1406.2332 [hep-ph]].
  * (38) W. Altmannshofer, S. Gori, S. Profumo and F. S. Queiroz, “Explaining dark matter and B decay anomalies with an $L_{\mu}-L_{\tau}$ model,” JHEP 12, 106 (2016) doi:10.1007/JHEP12(2016)106 [arXiv:1609.04026 [hep-ph]].
  * (39) G. Aad et al. [ATLAS], “Measurements of Four-Lepton Production at the Z Resonance in pp Collisions at $\sqrt{s}=$7 and 8 TeV with ATLAS,” Phys. Rev. Lett. 112, no.23, 231806 (2014) doi:10.1103/PhysRevLett.112.231806 [arXiv:1403.5657 [hep-ex]].
  * (40) A. M. Sirunyan et al. [CMS], “Search for an $L_{\mu}-L_{\tau}$ gauge boson using Z$\to 4\mu$ events in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” Phys. Lett. B 792, 345-368 (2019) doi:10.1016/j.physletb.2019.01.072 [arXiv:1808.03684 [hep-ex]].
  * (41) J. P. Lees et al. [BaBar], “Search for a muonic dark force at BABAR,” Phys. Rev. D 94, no.1, 011102 (2016) doi:10.1103/PhysRevD.94.011102 [arXiv:1606.03501 [hep-ex]].
  * (42) A. Kamada and H. B. Yu, “Coherent Propagation of PeV Neutrinos and the Dip in the Neutrino Spectrum at IceCube,” Phys. Rev. D 92, no.11, 113004 (2015) doi:10.1103/PhysRevD.92.113004 [arXiv:1504.00711 [hep-ph]].
  * (43) S. Gninenko and D. Gorbunov, “Refining constraints from Borexino measurements on a light Z’-boson coupled to L$\mu$-L$\tau$ current,” Phys. Lett. B 823, 136739 (2021) doi:10.1016/j.physletb.2021.136739 [arXiv:2007.16098 [hep-ph]].
  * (44) M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. F. Li, E. Picciau and Y. Y. Zhang, “Constraints on light vector mediators through coherent elastic neutrino nucleus scattering data from COHERENT,” JHEP 01, 116 (2021) doi:10.1007/JHEP01(2021)116 [arXiv:2008.05022 [hep-ph]].
  * (45) P. B. Denton, Y. Farzan and I. M. Shoemaker, “Testing large non-standard neutrino interactions with arbitrary mediator mass after COHERENT data,” JHEP 07, 037 (2018) doi:10.1007/JHEP07(2018)037 [arXiv:1804.03660 [hep-ph]].
  * (46) J. Liao and D. Marfatia, “COHERENT constraints on nonstandard neutrino interactions,” Phys. Lett. B 775, 54-57 (2017) doi:10.1016/j.physletb.2017.10.046 [arXiv:1708.04255 [hep-ph]].
  * (47) M. R. Buckley, D. Hooper, J. Kopp and E. Neil, “Light Z’ Bosons at the Tevatron,” Phys. Rev. D 83, 115013 (2011) doi:10.1103/PhysRevD.83.115013 [arXiv:1103.6035 [hep-ph]].
  * (48) M. Carena, A. Daleo, B. A. Dobrescu and T. M. P. Tait, “$Z^{\prime}$ gauge bosons at the Tevatron,” Phys. Rev. D 70, 093009 (2004) doi:10.1103/PhysRevD.70.093009 [arXiv:hep-ph/0408098 [hep-ph]].