\section{Introduction}
\label{Sec:introduction}

A model of two types of two-level systems (TLSs), interacting weakly and strongly with the vacuum field, describes both magnetic insulators having electronic and nuclear spins, and orientational glasses, in which inversion symmetric and inversion asymmetric excitations have a weak and strong interaction with the phonon field, respectively\cite{SS09}.

Such orientational glasses, and similarly amorphous solids, show remarkable universality in their low temperature characteristics; their specific heat increasing as $T^a$ with $a \approx 1$, the thermal conductivity increasing as $T^b$ with $b  \approx 2$ (both quantities increase as $T^3$ in ordered lattices), and their internal friction being temperature independent\cite{ZP71,HR86,PLT02}. The universality is also quantitative, as all characteristics that are dictated by phonon attenuation suggest a rather similar ratio of $\sim 150$ between phonon mean free path to phonon wavelength, in systems ranging from amorphous solids, to disordered lattices, polymers and porous aerogels\cite{PLT02}. All the above phenomena are present below a temperature of $1-3$K, which is, again, rather  universal. However, these systems show also characteristic phenomena at higher temperatures, such as a distinct functional dependence of the thermal conductivity, including a plateau between $3-10$K, and a bump in the specific heat dubbed the "boson peak"\cite{Phi81}.

Much of the universal behavior of disordered systems at low temperatures can be explained by the "Standard Tunneling Model"\cite{AHV72,Phi72,Jac72}. This model introduces phenomenologically tunneling two-level systems (TLSs), which interact weakly with the phonon field and negligibly between each other. TLSs are of fundamental and practical interest. First, they are the dominant dynamic degrees of freedom determining low temperature thermodynamics and kinetics in amorphous solids\cite{DH81,BGH90,RNO96}. Second, they are important for quantum information since they restrict the quantum coherence in Josephson junction qubits\cite{SLH+04,SSMM05,MCM+05}. Understanding the TLSs nature can help to reduce their destructive effect.

Although the standard tunneling model has been successful in explaining many of the above mentioned characteristics of disordered solids at low temperatures, it left some central questions unanswered, including the explanation of the quantitative universality between different systems, and the origin of the energy scale of $\approx 1-3$K below which universality exists. These questions have been the subject of theoretical scrutiny\cite{KFAA78,KKI83,SC85,AJL88,BNOK98,BGG+92,Par94,LW01,SK94,Kuhn03,PSG07,VL11} for the last four decades, yet with limited success.

A novel attempt to resolve these questions was recently made within a "two-TLS" model\cite{SS09}, consisting of two types of TLSs, weakly interacting and strongly interacting with the phonon field. A central property of the Two-TLS model is the structure of the low energy density of states (DOS) of the weakly interacting and strongly interacting TLSs. The weakly interacting TLSs have energies smaller than $\approx 10$K, and  are only weakly affected by TLS-TLS correlations. In contrast, the strongly interacting TLSs have typical energies of $\approx 300$K, but are significantly gapped at low energies. As a result, it is the weakly interacting TLSs that dictate the universal phonon attenuation at low temperatures and it is the width of their energy distribution, and the consequent gapping of the strongly interacting TLS DOS at the same energy scale, that dictates the energy scale of $\approx 1-3$K below which universal properties are observed.

Whereas the two-TLS model explains the above mentioned phenomena as they appear in disordered lattices and amorphous solids alike, for disordered lattices the two-TLS model is microscopically derived\cite{SS09}, is validated numerically using bare atomic interactions Hamiltonian\cite{GS11,CBS14}, and found to agree with experiment\cite{GS11}. Furthermore, for $T < 1-3$ K, where other degrees of freedom are gapped, the two-TLS model can be used to calculate, for disordered lattices, not only the universal phonon attenuation properties, but also properties such as the dephasing of superconducting qubits by TLSs, and phonon attenuation in non-equilibrium\cite{SNB14}. For $T > 3$ K other degrees of freedom, e.g. molecular librations, become significant and dominate the physics of orientational glasses\cite{RS88,GRS88}. However, also in this temperature range it is to be expected that tunneling states affect physical properties, and thus the two-TLS model, and specifically the bimodal structure of the DOS of the bias energies as is calculated below, may be found relevant.

With regard to amorphous solids, the general applicability of the two-TLS model is an open question of current research, which is, however, beyond the scope of this paper. Still, the two-TLS model, and specifically the form of the DOS of the weakly and strongly interacting TLSs found here, were already found relevant in amorphous solids under certain conditions: the existence of two types of TLSs was recently demonstrated, and was suggested to explain dielectric loss data in amorphous ${\rm Al_2O_3}$ and ${\rm LaAlO_3}$ films\cite{KSG+13}; and power law dependence of the DOS of the TLSs as we find here for the strongly interacting TLSs was recently found experimentally in amorphous SiO films\cite{SKW+14} and was used to explain decoherence rates of superconducting qubits\cite{FI14}.

The gapping of the strongly interacting TLSs is a manifestation of the Efros-Shklovskii gap for a system of two types of TLSs interacting via a random dipolar-like interaction. The Efros-Shklovskii gap of low energy single particle excitations in glassy systems with long range interactions\cite{ES75} has been a subject of thorough theoretical and experimental studies in various systems such as spin glasses\cite{BKS08}, electron glasses\cite{DLR82,VOP02,BSK+06,AOI11}, proteins\cite{BBFK07}, graphene nanoribbons\cite{DKM+11}, and amorphous solids\cite{AJL88,Bur95,RNO96}. The particular behavior of the gap dictates thermodynamic and transport properties. Thus, and since the two-TLS model suggests that the low energy excitations in orientationsl glasses, and possibly also in amorphous solids, are generically given by the weakly and strongly interacting TLSs as described in the Hamiltonian (\ref{H-St}),(\ref{Jcutoff}) below, the rigorous calculation of their DOS is of fundamental interest.

In this paper we derive analytically, and verify numerically, the single particle DOS of the weakly and strongly interacting TLSs as are given by the Hamiltonian in Eqs. (\ref{H-St}),(\ref{Jcutoff}). The TLSs interact via the acoustic interaction, which is modeled as having a $1/r^3$ spatial dependence with a short distance cutoff, and a random angular dependence.
Whereas for a single species the Efros Shklovskii correlation gap is derived by a self-consistent calculation for the DOS\cite{ES75,BSE80}, the two-TLS structure as presented in the Hamiltonian (\ref{H-St}),(\ref{Jcutoff}) allows analytical derivation within controlled approximations, not invoking self-consistency. Our results are then confirmed using Monte-Carlo simulations. We find that the DOS of the weakly interacting TLSs is well described by that of a random field Ising model, with an effective field much larger than the interaction, as the former is a result of their interaction with the strongly interacting TLSs. The DOS of the strongly interacting TLSs has a {\it power law} dependence on energy at small energies, with a power that depends on the value of the short distance cutoff of the acoustic interaction. This dependence of the DOS on the form of the interaction at short distances is quite remarkable, as the physics of the correlation gap at low energies is that of long distances. We show that it is the product of the interaction energy constant at large distances and the DOS of the weakly interacting TLSs at low energies that dictates the gap behavior. This product depends on the short distance cutoff, because with enhanced cutoff the distribution of energies of the weakly interacting TLSs becomes narrower, and their DOS at low energies becomes larger.

The paper is organized as follows: In Sec.~\ref{Sec:Analytical} we present the model and the analytical derivation of the DOS of the weakly interacting and strongly interacting TLSs. Our numerical results are presented in Sec.~\ref{Sec:Numerical}. We then conclude with a discussion of our results and their consequences in Sec.~\ref{Sec:Discussion}. Some details of the analytical derivation of the DOS are deferred to App.~\ref{Sec:tauDOS} (weakly interacting TLSs) and App.~\ref{Sec:SDOS} (strongly interacting TLSs). The numerical calculation of the power law in the DOS of the $S$-TLSs is detailed in App.~\ref{Appnumerical}, and the effect of finite size is considered in App.~\ref{Appfinitesize}.

\section{Analytical theory}
\label{Sec:Analytical}

The two-TLS model assumes a bi-modality of the interaction of TLSs with the phonon field. This form of the interaction, which, for the disordered lattices, is rigorously derived\cite{SS09} and numerically validated\cite{GS11,CBS13}, leads to acoustic mediated interaction between tunneling TLSs in disordered systems that is described by the Hamiltonian\cite{SS08,SS09}

\begin{multline}
H_{S \tau} \;=\; -\sum_{i \neq j} \left[\frac{1}{2}J_{ij}^{SS} S_i^z S_j^z +
J_{ij}^{S\tau} S_i^z \tau_j^z + \frac{1}{2}J_{ij}^{\tau \tau} \tau_i^z
\tau_j^z \right] \, ,
 \label{H-St}
\end{multline}
where
\begin{equation}
J_{ij}^{ab}=c_{ij}^{ab} \frac{J_o^{ab}}{r_{ij}^3 + \tilde{a}^3} \, .
\label{Jcutoff}
\end{equation}
$a,b=S$ or $\tau$ stand for the strongly interacting and weakly interacting TLSs, and are denoted by pseudo Ising spins with $S^z=\pm 1, \tau^z=\pm 1$. $c_{ij}^{ab}$ is chosen randomly from a Gaussian distribution of width unity. The factor of $1/2$ in the first and third terms of the Hamiltonian (\ref{H-St}) accounts for the double summation.
$J_o^{ab}$ denote the interaction energy scales at nearest neighbor distance, and are related by $J_o^{\tau \tau} = g J_o^{S \tau} = g^2 J_o$\cite{SS09}, where we define $J_o \equiv J_o^{SS}$ and the dimensionless parameter $g \ll1$. This separation of energy scales, which distinguishes between the weakly interacting $\tau$-TLSs and the strongly interacting $S$-TLSs, is central to our analysis. In orientational glasses $J_o$ is typically of the order of the Debye energy, and $g \approx 0.01-0.03$\cite{SS09,GS11,CBS13}.
Both the elastic interaction in glasses and the magnetic dipolar hyperfine interaction decay at large distances as $1/r^3$, but have a different spatial dependence at distances of the order of the interatomic distance. We therefore introduce a short distance cutoff $\tilde{a}$ through Eq.(\ref{Jcutoff}). We expect other forms of the cutoff, that eliminate the divergence of the interaction at small distances, to give similar results. Here and throughout the paper we assume that the $S$-TLSs and $\tau$-TLSs are randomly placed in a lattice with concentration (density per site) $\rho$, and distances will be given in units of the lattice spacing $a_o$.

We emphasize that in this paper we are concerned only with the bias energies $E$ (usually denoted $\Delta$) between the two states of each of the TLSs, and not with the tunneling amplitudes $\Delta_o$ between these states, which can be approximately ignored for the dipole gap analysis. As such, the DOS calculated for both the weakly interacting TLSs and the strongly interacting TLSs is that of the bias energies. The full energy of a TLS is given by $\sqrt{E^2+\Delta_o^2}$. We then postulate the distribution functions in the standard form

\begin{equation}
P^{(\tau,S)}(E,\Delta_o)=P^{(\tau,S)}_o(E)/\Delta_o \, ,
\end{equation}
following the ansatz of the standard tunneling model for the independence of the bias energies and the tunneling amplitudes, but allowing an energy dependence of the DOS of the bias energies, which we calculate below. We note that while the typical values of the tunneling amplitudes is limited below a few Kelvin, their range
(i.e. $\Delta_o^{min}$ and $\Delta_o^{max}$) could differ between the $\tau$-TLSs and S-TLSs. However, the $1/\Delta_o$ dependence of the distribution function rests on the approximate homogeneity of the tunnel barriers, which is a generic consideration similarly applicable to both classes of tunneling states.

The Hamiltonian we consider [in Eq. (\ref{H-St})] is therefore classical, and the bias energies are given by the single particle excitation (spin flip) energies for the $S$-TLSs and for the $\tau$-TLSs as

\begin{equation}
E_S^j \equiv 2\sum_{i(\neq j)}[J_{ij}^{SS} S_i^z S_j^z +
J_{ij}^{S\tau} S_j^z \tau_i^z]
\label{Senergy}
\end{equation}
and
\begin{equation}
E_{\tau}^j \equiv 2\sum_{i(\neq j)}[J_{ij}^{S\tau} S_i^z \tau_j^z + J_{ij}^{\tau \tau} \tau_i^z
\tau_j^z] .
\label{tauenergy}
\end{equation}
The factor of $2$ is a result of the spin changing by $2$ upon flipping.
We then calculate for the $\tau$-TLSs and for the S-TLSs the functional dependence of the single particle DOS at zero temperature, $n_{\tau}(E)$, $n_S(E)$. The uncorrelated DOS for the $S$-TLSs and $\tau$-TLSs is given by a Gaussian distribution with typical energy and peak values dictated by the dominating interaction ($S-S$ for the $S$-TLSs and $S-\tau$ for the $\tau$ TLSs). However, correlations in between the $S$-TLSs, between the $S$-TLSs and the $\tau$ TLSs, and in between the $\tau$-TLSs result in a depression in $n_S(E)$ and in $n_{\tau}(E)$ at low energies. This depression is governed by the Efros-Shklovskii stability criterion\cite{ES75}, which for the $S-S$ correlations reads:

\begin{equation}
E_S^i + E_S^j - 4J_{ij}^{SS} S_i^z S_j^z > 0 \,
\label{ESSS}
\end{equation}
for any two $S$-TLSs in the system, while for the $S-\tau$ correlations one requires

\begin{equation}
E_S^i + E_{\tau}^j - 4J_{ij}^{S\tau} S_i^z \tau_j^z > 0 \,
\label{ESStau}
\end{equation}
for any $S$ and $\tau$ TLSs, and for the $\tau-\tau$ correlations one has

\begin{equation}
E_{\tau}^i + E_{\tau}^j - 4J_{ij}^{\tau\tau} \tau_i^z \tau_j^z> 0 \, .
\label{EStautau}
\end{equation}

These conditions are a manifestation of the requirement that the ground state must be stable to flips of any pair of spins. Due to the many-body nature of spin-spin interactions the energy to flip such a pair can be lower than the sum of the single TLS flip energies. The difference is represented by four times the intra pair interaction term (a factor of $2$ for each single spin flip as in Eqs. (\ref{Senergy}),(\ref{tauenergy}), and another factor of $2$ comes from the fact that this interaction energy is counted once for each single flip, but not at all for the double flip). For a single species of TLSs interacting via the dipolar interaction, e.g. for a system described by the Hamiltonian given by the first term in Eq.(\ref{H-St}) with $\tilde{a}=0$, correlations lead to logarithmic depression of the DOS, which for low energies $E$ is proportional to $1/\log{(E_o/E)}$\cite{BSE80,Bur95}. The same of course would be true for the $\tau$-TLSs in the absence of the $S$-TLSs. However, the presence of two types of TLSs, with significant difference in their coupling as is given in Eq.(\ref{H-St}), changes things remarkably.

\subsection{DOS of the weakly interacting TLSs}

Let us consider first $n_{\tau}(E)$. As is shown below, $S$-TLS excitations are scarce at low energies, and thus, for the consideration of the $\tau$-TLSs, to first approximation, and up to an overall constant (resulting from the $S-S$ interactions),
the Hamiltonian in Eq.(\ref{H-St}) reduces at low temperatures to the effective Hamiltonian

\begin{multline}
H_{\tau}^{RF} \;=\; -\sum_{i \neq j} \left[J_{ij}^{S\tau} \langle S_i^z \rangle \tau_j^z + \frac{1}{2} J_{ij}^{\tau \tau} \tau_i^z
\tau_j^z \right] \\ \equiv -\frac{1}{2} \sum_{i \neq j} J_{ij}^{\tau \tau} \tau_i^z \tau_j^z + \sum_j h_j \tau_j^z \, ,
\label{tauRFIM}
\end{multline}
which is equivalent to the random field Ising model with the random field $h_j = -\sum_{i(\neq j)} J_{ij}^{S\tau} \langle S_i^z \rangle$ being $\approx 1/g$ times larger than the $\tau-\tau$ interaction.
The large random field leads to a substantial reduction of the effect of the $\tau-\tau$ correlations.  The effect of the correlations with the $S$-TLSs [present in the full Hamiltonian (\ref{H-St})] is even much weaker because of the smallness of $n_S(E)$ at low energies, as is derived below.
In Ref.\cite{SS09} it was shown that $n_{\tau}(E)$ dips because of correlations only at energies smaller than $g^2 J_o$, and that the relative decrease in DOS is small, proportional to $g$. Apart from this decrease we show in App.\ref{Sec:tauDOS} that for intermediate TLS spatial concentrations $\rho$, i.e. $[z(1+\tilde{a}^3)]^{-1} < \rho \lesssim 1$, $n_{\tau}(E)$ is well approximated by a Gaussian distribution\cite{Lorenote}.
Here $z$ is the coordination number of the lattice (or the underlying lattice in the amorphous state), and all distanced are in units of interatomic spacing. The left inequality ensures that the few largest contributions (coming from $\tau-S$ interactions) to any given $\tau$-TLS are of similar magnitude. The second inequality assures strong enough disorder. The width of the distribution, $E_{\tau}^{typ}$, is dictated by the Hamiltonian (\ref{tauRFIM}). In the limit of $\tilde{a} \gg a_o$ we find

\begin{equation}
n_{\tau}(E) = \frac{2 \rho}{\sqrt{2\pi}E_{\tau}^{typ}}e^{\frac{E^2}{2(E_{\tau}^{typ})^2}}
\end{equation}
with
\begin{equation}
E_{\tau}^{typ}=\sqrt{\frac{16\pi}{3} \frac{\rho}{\tilde{a}^3}} g J_o \, .
\label{E-typ}
\end{equation}
The analysis of the distribution in detail is given in App.\ref{Sec:tauDOS}. Both the Gaussian functional dependence of $n_{\tau}(E)$ and the analytical results for the widths of the distribution as function of $\tilde{a}$ are verified numerically, see Sec. \ref{Sec:Numerical}.

\subsection{DOS of the strongly interacting TLSs}

We now turn to the $S$-TLSs. At energies much larger than $gJ_o$, i.e. energies larger than the maximum energy of a $\tau$-TLS excitation, we expect the $S$-TLS DOS to have the same functional behavior, i.e. logarithmic depression, as in the absence of the $\tau$-TLSs. This is confirmed numerically in Sec.\ref{Sec:Numerical}. However, at low energies, where $S$-TLSs are scarce and $\tau$-TLSs are abundant, $S-\tau$ correlations as manifested in Eq.(\ref{ESStau}) dominate the functional dependence of $n_S(E)$. We note also that it is the low energy regime which is of most interest to us, as in the related regime of low temperatures TLSs determine the physical properties of glasses.
We thus calculate the $S$-TLS DOS resulting from Eq.(\ref{ESStau}) analytically for energies $E \ll gJ_o$.

\begin{figure*}
\subfigure[]{
  \includegraphics[width=8.6cm]{Fig1a_r.eps}
   \label{fig:Fig1subfig1}
   }
 \subfigure[]{
  \includegraphics[width=8.6cm]{Fig1b_r.eps}
   \label{fig:Fig1subfig2}
   }
 \subfigure[]{
  \includegraphics[width=8.6cm]{Fig1c_fit.eps}
   \label{fig:Fig1subfig3}
   }
 \subfigure[]{
  \includegraphics[width=8.6cm]{Fig1d_fit.eps}
   \label{fig:Fig1subfig4}
   }
\caption{ Single particle DOS for cube of volume $12\times12\times12$ with $\rho=0.5$ and $\tilde{a}=2$. [Here and in all subsequent figures DOS is plotted per unit cell (i.e. per site)](a) $\tau$-TLSs DOS with a zero free parameter Gaussian fit, using the standard deviation calculated in App. \ref{Sec:tauDOS}. Inset shows the small shift of states from the lowest energies upwards. (b) $S$-TLSs DOS.  (c)  Zoom into low energies of the $S$-TLSs DOS in the absence and presence of $\tau$-TLSs. The shift from log to power law behavior as a result of the $S-\tau$ interactions is presented. Solid lines denote a fit to the function $A/Log(B/E)$ [$A=250, B=0.00038$, upper curve], and to the analytical result for the $S$-TLS DOS at low energies, [Eq.(\ref{PE6}) with $\alpha=2.65$, lower curve]. (d) The ratio between the $S$-TLSs DOS with $\tau$ TLSs and
the $S$-TLSs DOS  in the absence of the $\tau$ TLSs, zoomed into low energies.}
\label{fig:Fig1}
\end{figure*}

We define $P(E)$ by

\begin{equation}
n_S(E) = P(E) \bar{n}_S(E)
\label{PEdef}
\end{equation}
where $\bar{n}_S(E)$ is the $S$-TLS DOS neglecting $S-\tau$ correlations, and calculate explicitly $P(E)$.
We use the Efros-Shklovskii criterion of stability

\begin{equation}
E_{\tau}+E_S-\frac{4\eta gJ_o}{R^3+\tilde{a}^3} > 0
\label{E_S}
\end{equation}
which, for each $S$-TLS, should be fulfilled for all $\tau$-TLSs. Thus, for a given $S$-TLS, the reduction factor in its probability to be at an energy $E$ coming from fulfilling Eq.(\ref{E_S}) for all $\tau$-TLSs is given by

\begin{multline}
P(E)=\prod_{d{\bf R}} \biggl[ 1-d{\bf R}\int_{\frac{E\tilde{a}^3}{4gJ_o}}^{\infty} d\eta \frac{1}{\sqrt{2\pi}} e^{-\eta^2/2} \\ \cdot
\int_0^{4\eta gJ_o/\tilde{a}^3} dE_{\tau} n_{\tau}(E) \Theta\left(-E-E_{\tau}+\frac{4\eta gJ_o}{R^3+\tilde{a}^3}\right) \biggr]  \, .
\label{PE1}
\end{multline}

Since we are interested in low energies we approximate $n_{\tau}(E)$ by a constant given by the $\tau$-TLS DOS at zero energy $n_{\tau}(0)$. Substituting $z\equiv R^3$ we obtain

\begin{multline}
P(E)=exp\biggl[-\frac{4\pi}{3}\int_0^{\infty}dz\int_{\frac{E\tilde{a}^3}{4gJ_o}}^{\infty} d\eta \frac{1}{\sqrt{2\pi}} e^{-\eta^2/2} \\
\cdot \int_0^{4\eta gJ_o/\tilde{a}^3} dE_{\tau} n_{\tau}(0) \Theta\left(-E-E_{\tau}+\frac{4\eta gJ_o}{z+\tilde{a}^3}\right)\biggr] \, .
\label{PE2}
\end{multline}

The detailed evaluation of $P(E)$ is given in App.\ref{Sec:SDOS}. We find that for all values of $\tilde{a}$ the DOS of the $S$-TLSs has a power law dependence on energy at low energies

\begin{equation}
n_S(E) = \bar{n}_S \left( \frac{cE}{gJ_o} \right)^{\alpha} \, ,
\label{PE6}
\end{equation}
where $\alpha$ is given by

\begin{equation}
\alpha=\frac{8\sqrt{2\pi} n_{\tau}(0) gJ_o}{3} \, ,
\label{alphaanalytical}
\end{equation}
and $c$ is a constant of order unity, a result of the ability to calculate the exponent in Eq.~(\ref{PE2}) only within logarithmic accuracy.
The dependence of $\alpha$ on the interaction cutoff $\tilde{a}$ is given through $n_{\tau}(0)$. A general derivation of $n_{\tau}(0)$, along with the values of $\alpha$ for various cutoffs $\tilde{a}$ and concentrations $\rho$, is given in App. \ref{Sec:tauDOS}, see also Tab. \ref{tab:Table2}.

For large short distance cutoffs, $\tilde{a} \gg a_o$ we find

\begin{equation}
n_{\tau}^{\tilde{a}}(0) = \frac{\sqrt{3/2}\sqrt{\rho \tilde{a}^3}}{4\pi gJ_o} \, .
\label{ntauzero}
\end{equation}
We assume here, for simplicity, $\rho_S=\rho_{\tau}=\rho$, but generalization to $\rho_S \neq \rho_{\tau}$ is straightforward, see App. \ref{Sec:tauDOS}. We can then write the $S$-TLS exponent as

\begin{equation}
\alpha=\frac{2 \sqrt{\rho \tilde{a}^3}}{\sqrt{3\pi}} \, .
\end{equation}
With decreasing $\tilde{a}$, $\alpha$ decreases monotonically, reaching the value of $\alpha=0.92 \sqrt{\rho}$ for $\tilde{a}=0$.

\section{Numerical calculations}
\label{Sec:Numerical}

\begin{figure*}[!tbh]
\centering
\subfigure[]{
  \includegraphics[width=8.6cm]{Fig2a.eps}
   \label{fig:Fig2subfig1}
   }
 \subfigure[]{
  \includegraphics[width=8.6cm]{Fig2b_fits.eps}
   \label{fig:Fig2subfig2}
   }
\caption{$S$-TLSs DOS for a cube of volume $14\times14\times14$ with $\rho=0.25$ and different values of interaction cutoff, $\tilde{a}=$0, 1, 2. $g=1/30$ for all cutoffs. $J_o$ is chosen for each value of $\tilde{a}$ such that the peak of the distribution is cutoff independent (see text). (a) $S$-TLSs DOS - entire plot. (b)  Zoom into low energies: 0K - 30K. The increase of the power $\alpha$ with increasing values of interaction cutoff is clearly observed (fits done with $\alpha$ values calculated analytically, as are given in Table~\ref{tab:Table2}).}
\label{fig:Fig2}
\end{figure*}

Our analytical derivation of the DOS relies on the Efros Shklovskii criterion for the low energy DOS of the $S$-TLSs, and on the assumption of a Gaussian distribution for all but the low energy excitations of the $\tau$-TLSs. In order to check our analytical results and verify the validity of the above assumptions we perform a numerical calculation of the DOS of both the $\tau$-TLSs and the $S$-TLSs starting from the Hamiltonian in Eqs.(\ref{H-St}),(\ref{Jcutoff}).

\subsection{DOS of $S$ and $\tau$ TLSs}

Calculations are performed on cubic lattices of size $L$, with $L=6-14$ and periodic boundary conditions, where the interaction between each pair of TLSs is calculated according to their shortest separation in the extended lattice. $S$ and $\tau$ TLSs are placed randomly in the lattice with concentration $\rho$ where the cases of $\rho=0.25$ and $\rho=0.5$ are analyzed. In accordance with the model in Ref. \cite{SS09} we take each occupied site to contain both an $S$ and a $\tau$ TLS, and avoid on site interactions, but taking the positions of the $S$ and $\tau$ TLSs to be uncorrelated produces similar results. Using Monte Carlo simulation with the Hamiltonian (\ref{H-St}),(\ref{Jcutoff}) we lower the free energy of the system at temperatures decreasing from $300$K to $0.01$K. As in realistic systems, our simulations reach very low energy states at the lowest temperature, yet the system does not equilibrate. Such low energy states produce the correct Efros Shklovskii gap once they are stable to single and double spin flips\cite{BSE80}. Both conditions are explicitly checked, and are very well satisfied in our simulations at $T=0.01$K. Once we reach the final state of the simulation at a given size, dilution, and interaction cutoff, we measure the excitation energy of each $S$ and $\tau$ TLS.

In Fig.~\ref{fig:Fig1}(a),(b) we show the single particle DOS of the $\tau$-TLSs and of the $S$-TLSs, calculated for a cube of volume $12\times12\times12$ with impurity concentration $\rho=0.5$ and short distance cutoff $\tilde{a}=2$, with the parameters $J_o=300$K and $g=1/30$. The typical energy scale for $S$ excitations is indeed $J_o$, a result of the interactions between different $S$-TLSs. At lower energies $n_S(E)$ is logarithmically reduced. The typical energy scale of the $\tau$-TLSs is of order $gJ_o$, a result of their interactions with $S$-TLSs.
In Fig.\ref{fig:Fig1}(a) we compare our data to the analytical results presented in Sec.~\ref{Sec:Analytical}, i.e. to a Gaussian with standard deviation of $1.008gJ_o$, see App.\ref{Sec:tauDOS}, Table~\ref{tab:Table1}. We obtain an excellent fit, with no free parameters,
except at very low energies, smaller than $g^2 J_o$, where it is slightly diminished (by relative magnitude $\approx g$) because of $\tau-\tau$ correlations\cite{SS09}. Furthermore, the peak value of $n_{\tau}(E)$ is $\approx 1/g$ times larger than the peak value of $n_S(E)$, as expected\cite{SS09}. This validates the fact that $\tau$-TLSs can be described by the random field Ising model, see Eq.~(\ref{tauRFIM}), with the exception of a small correction at low energies.

The same energy scale, $gJ_o$, which marks the onset of $\tau$-TLSs, marks also the sharp decay in $n_S(E)$.
In Fig.\ref{fig:Fig1}(c) we zoom into low energies. We see that at the energy where the $\tau$ TLSs appear ($\approx 10$K) the reduction of the $S$ TLS DOS changes its functional form to a power law,
as is demonstrated by a fit of the low energy data to Eq.(\ref{PE6}), with $\alpha=2.65$ (see Table~\ref{tab:Table2}).
The fact that this power law gap is a result of the $S-\tau$ correlations is further demonstrated by the calculation of the $S$-TLS DOS in the absence of the $\tau$ TLSs [i.e. taking $J_o^{S \tau} = J_o^{\tau \tau} =0$ in the Hamiltonian (\ref{H-St})]. This graph is shown for comparison in Fig.\ref{fig:Fig1}(c). Indeed, in the absence of $\tau$-TLSs the logarithmic gap continues to low energies. In Fig.\ref{fig:Fig1}(d) we plot the ratio between the DOS of the $S$-TLSs in the absence and in the presence of the $\tau$-TLSs, singling out the effect of the $S-\tau$ correlations on the $S$-TLS DOS.

\subsection{Dependence of the DOS on the short distance cutoff}

In Fig.~\ref{fig:Fig2} we plot $n_S(E)$ for a cube of volume $12\times12\times12$ with different values of interaction cutoffs. For each value of the cutoff $J_o$ is chosen in the way that the position of the peak of $n_S(E)$ is cutoff independent, and we keep $g=1/30$ independent of the cutoff.
At low energies we indeed find a deepening of the gap with the power $\alpha$ increasing with increasing cutoff $\tilde{a}$, with the excellent fit with the values of $\alpha$ obtained analytically in Sec. \ref{Sec:Analytical}. We emphasize that the value of the power $\alpha$ dictating $n_S(E)$ at low energies does not depend on our choice of $J_o$, as can be inferred from Eq.~(\ref{E-typ}) and Eq.~(\ref{alphaanalytical}) above (noting that $n_{\tau}(0) \propto 1/J_o$).

\begin{figure}[h!]
\centering
\includegraphics[width=8.6cm]{notau.eps}
\caption{$S$-TLSs DOS plots in the absence of the $\tau$ TLSs for cube of volume $12\times12\times12$ with $\rho=0.5$ and $\tilde{a}=0,2,4$. Enhanced cutoffs diminish the low energy DOS, but do not change its functional dependence.}
\label{fig:notau}
\end{figure}

The power $\alpha$ can also be extracted from our numerical results by analyzing the integral plots of the $S$-TLS  DOS. This is done in App.\ref{Appnumerical} for various cutoff parameters, with $\rho=0.25$ and $\rho=0.5$.
Despite finite size effects (see App.\ref{Appfinitesize}) our numerical results for the value of $\alpha$ are in reasonable agreement with the analytical results both in absolute value (see Tab.~\ref{tab:Table2}), and in the functional dependence of $\alpha \propto \sqrt{\rho \tilde{a}^3}$ predicted analytically for $\tilde{a} \gg a_o$.

The specific parameters chosen for $J_o$ and $g$ are in accordance with reasonable values for amorphous solids and disordered lattices that express the low energy universal characteristics\cite{ZP71,HR86,PLT02}, and specifically with calculated values for KBr:CN\cite{GS11,CBS13}. Our results show that with these parameters $N_S(E) = g^2 N_{\tau}(E)$ at $E \approx 3$K for cutoffs $\tilde{a}$ of order unity, in agreement with experiment.

The short distance cutoff affects also a system consisting of a single species of TLSs, i.e. having only the first term in the Hamiltonian~(\ref{H-St}). However, in this case the short distance cutoff only changes the magnitude of the DOS at low energies, but not their functional form.  This can be seen in Fig.\ref{fig:notau}, where $N_S(E)$ in the absence of $\tau$ TLSs is plotted for $\tilde{a}=0,2,4$.

\section{Discussion}
\label{Sec:Discussion}

We have derived analytically and numerically the DOS of the bias energies of the weakly ($\tau$) and strongly (S) interacting TLSs within the two-TLS model. We find that the $\tau$-TLSs are confined to energies smaller than $J_o^{S \tau}$ $(\approx 10$ K), whereas the $S$-TLSs are spread to much larger energies, and have a power law gap at low energies.
At energies smaller than $\approx 3$ K, where the $S$-TLSs are scarce, the DOS of the $\tau$ TLSs is, to a very good approximation, constant, and in a wide regime of energies $J_o^{\tau \tau} < E < J_o^{S \tau}$ the $\tau$ TLSs are practically non interacting.
As such, the characteristics of the $\tau$-TLSs at low energies are equivalent to those introduced phenomenologically in the Standard Tunneling Model\cite{AHV72,Phi72}. However, within the Two-TLS model [Eq.~(\ref{H-St})] these characteristics are derived, and the relation between the density of states and the coupling to the phonon field is obtained in terms of the small parameter $g$ (Ref. \cite{SS09}). This allows for the derivation of the magnitude $C_o \approx 10^{-3}$ of the tunneling strength\cite{PLT02} and its universality, and the energy scale of $gJ_o\approx 3$ K dictating the temperature below which universality is observed.

Our results for the DOS of the $S$-TLSs can be used to facilitate predictions for the behavior of various properties of orientational glasses beyond the phenomena of the low temperature universality, both below and above the temperature of $3$ K. Below $3$ K, whereas the number of thermal $S$ TLSs is small, their strong interaction with the phonon field can lead to their dominance of properties in which this interaction contributes with a power larger than two\cite{SNB14}. At temperature larger than $3$ K it is expected that TLS contribution will be dominated by the $S$-TLSs, resulting in added contribution to that of other prevalent excitations (e.g. librations) at these energy scales, both directly by $S$-TLSs, and through their interactions with e.g. librational modes. Our results here could therefore be useful to the discussion of phenomena such as the plateau in the thermal conductivity at $3-10$ K and the boson peak, as the introduced $S$-TLSs enhance phonon attenuation and the specific heat in the relevant temperature regime.

The two-TLS model was rigorously derived\cite{SS09} and thoroughly validated\cite{GS11,CBS14,CBS13} for the disordered lattices, which constitute a significant subclass of systems showing universality. At the same time, it provides an explanation for crucial aspects of the low temperature universality as are exhibited in both disordered and amorphous solids. Given the equivalence of the phenomenon as observed experimentally in all systems showing the low temperature universality\cite{YKMP86,LVP+98}, we believe it is plausible that the Hamiltonian in Eqs. (\ref{H-St}),(\ref{Jcutoff}) below describes the DOS of the low energy excitations in all amorphous solids. A direct way to check this possibility would be by detecting the existence of the strongly interacting TLSs, and comparing their DOS to our results here. This could be done at high energies using e.g. TeraHertz absorption and spectroscopy experiments\cite{PHC+13}, and at low energies using the recently acquired ability to study single TLSs via their interaction with phase qubits and with strain\cite{LMC+10,GPL+12}, which allows distinction between weakly and strongly interacting TLSs.
Such verification of the applicability of our results to amorphous solids
could give strong support to the general applicability of the Two-TLS model in describing the low temperature universality in glasses\cite{SS09}. It could also lead to an enhancement of our understanding of the microscopic nature of amorphous solids, and its relation to the detailed characteristization of the TLSs, and to the properties of amorphous solids at low temperatures.

It is exciting that the presence of weakly coupled TLSs leads to the dramatic reduction in the low energy DOS of the strongly coupled TLSs. In principle, this opens the opportunity to control and reduce the number of relatively strongly coupled TLSs at low energies by somehow adding TLSs which are more weakly coupled, and in this way reduce the destructive absorption and decoherence effects caused by the TLSs. Although this suggestion does not seem to be easy practically because the nature of TLSs is unclear, one way of introducing weakly coupled TLSs can be associated with the hydrogenation of the material. Hydrogen atoms are expected to easily participate in tunneling and indeed they introduce the weakly coupled TLSs identified in aluminum and beryllium oxide glasses in Ref.\cite{KSG+13}. Moreover, hydrogenation of silicon oxide\cite{LWP+97} and silicon nitride\cite{PO10} results in the remarkable reduction of TLS induced absorption of sound or electromagnetic waves which can be due to the TLS gapping by hydrogen induced tunneling defects as described in this work. The investigation of this problem is a matter of current research.

Our results can also be carried through to magnetic insulators, as our model in Eqs.(\ref{H-St}),(\ref{Jcutoff}) describes the interactions of electronic and nuclear spins (the electron spin-spin interaction, the electron nuclear hyperfine interaction, and the nuclear spin-spin interaction are described by the first, second and third terms of Eq.~(\ref{H-St}), respectively). In particular, our results suggest a remarkable reduction in electronic spin flip rate in random magnetic systems at temperatures corresponding to the thermal energy being smaller than the typical hyperfine interaction. This implies a corresponding reduction in the decoherence of spin qubits at very low temperatures, and is a subject of a separate study.

{\it Acknowledgements ---} We would like to thank Ariel Amir, Juan Carlos Andresen, Danny Barash, Doron Cohen, Yoseph Imry, and Helmut Katzgraber for useful discussions. M. S. acknowledges financial support from the ISF and from the Marie Curie Grant PIRG-GA-2009-256313. A. B. acknowledges financial support from the Louisiana Board of Regents PFund project no. PFUND-330, Army Research Office grant number vv911NF-13-1-0186, and LINK Program of NSF and Louisiana Board of Regents for supporting his travel to Israel.


\section{Introduction}

In Ref.~\cite{Carrazza:2013wua} we have shown, using the latest NNPDF
technology~\cite{Ball:2012cx}, that the photon parton distribution
function (PDF) obtained from deep inelastic scattering (DIS) data has
large uncertainties. Furthermore, by computing with {\tt HORACE}
~\cite{CarloniCalame:2007cd} the $Z\rightarrow\mu^{+}\mu^{-}$
production in proton-proton collision, we have confirmed that the
photon PDF from DIS is not sufficiently precise for the level of
accuracy required by the LHC physics.

The large uncertainties of the photon PDF from DIS are due to the lack
of photon-induced contributions to DIS processes, when including QED
corrections up to leading order (LO) in $\alpha$ (i.e. $O(\alpha)$) to
PDFs. A possible cure for this issue is the inclusion of Drell-Yan
data where there are photon-induced contributions at $O(\alpha)$: for
example in $Z/\gamma^{*}$ production the $\gamma\gamma \rightarrow
l^{+} l^{-}$ process.

Here we present the determination of the photon PDF, for a PDF set
which is determined at NLO in QCD and LO in QED, by including a subset
of Drell-Yan data from LHC. Currently, fast interface codes for QCD
computations do not include contributions from photon PDFs, so we
decided to include these data through a reweighting procedure. In
short, the reweighting assigns weights to PDF replicas which satisfies
the new data by computing the $\chi^{2}$ between the new experimental
data and theoretical predictions obtained from each replica. Technical
details about the reweighting and unweighting procedure are presented
in Ref.~\cite{Ball:2010gb,Ball:2011gg} and will not be discussed here.

\section{Reweighting with ATLAS data}

We start the reweighting procedure by creating a prior PDF set with
500 replicas, combining the photon PDF from DIS with five copies of
the NNPDF2.3 NLO set~\cite{Ball:2012cx} at the initial scale
$Q^{2}_{0}=2$ GeV$^{2}$ and then evolving with the combined QCD+QED
DGLAP evolution equations. The prior PDF set thus defined preserves
the prediction accuracy of a global set for pure QCD computations:
indeed we have explicitly verified that PDF sets with and without QED
corrections are statistically equivalent which means that the
corrections to the quark and gluon due to their mixing with the photon
are indistinguishable from statistical fluctuations. Also, the
violation of the momentum sum rule due to the inclusion of
the photon PDF is less than 1\%.

We perform reweighting with two LHC data sets: ATLAS $W$,
$Z/\gamma^{*}$ rapidity distributions~\cite{Aad:2011dm} and ATLAS
Drell-Yan high mass differential
cross-section~\cite{Aad:2013iua}. These data sets provide constraints
to the photon PDF at medium and large-$x$ thanks to the presence of
photon-induced contributions in the $Z/\gamma^{*}$ and $W$
production. The QCD computation is performed with {\tt
  APPLGRID}~\cite{Carli:2010rw} for ATLAS $W$, $Z/\gamma^{*}$ rapidity
distributions data and {\tt DYNNLO}~\cite{Catani:2009sm} for ATLAS
Drell-Yan high mass data, while, in both cases we have computed the
Born and NLO photon-induced contributions with {\tt
  HORACE}~\cite{CarloniCalame:2007cd}, knowing that NLO contributions
have a very small effect.

Figure~\ref{fig:weights} shows the weights and $P(\alpha)$
distributions, defined in~\cite{Ball:2010gb}, for the present
reweighting. Recall that $\alpha$ is the rescaling of the covariance
matrix: $\alpha=1$ means that uncertainties are correctly
estimated. From these distributions we conclude that the procedure is
consistent and that the data is compatible with the reweighted PDF
set. After performing the reweighting we obtain a total number of
effective replicas of $N_{\rm eff}=266$ and the $\chi^{2}$ to the
ATLAS data is reduced from 1.98 to 1.14. In Table~\ref{table:chi2} we
collect the $\chi^{2}$ and $P(\alpha)$ values for the individual
reweighting procedure, looking at the $\chi^{2}$ values before and
after reweighting we confirm that data far from the $W$,
$Z/\gamma^{*}$ peaks, where the photon contribution is substantial,
contributes strongly to constrain the photon PDF.

\begin{table}
  \centering
  \begin{tabular}{c||c||c|c}
    \hline
    NLO &  ATLAS all  & ATLAS $W$, $Z/\gamma^{*}$  & ATLAS DY high mass  \\
    \hline
    \hline
    $\chi^2_{\rm in}$  & 1.98 & 1.20 & 3.78 \\
    \hline
    $\chi^2_{\rm rw}$  &  1.14 &  1.15 & 1.01 \\
    \hline
    $N_{\rm eff}$  & 266 & 364 & 326 \\
    \hline
    $\langle \alpha \rangle$  & 1.48 & 1.24 & 1.53 \\
    \hline
  \end{tabular}
  \caption{\label{table:chi2} Reweighting properties: $\chi^{2}_{\rm in}$ and $\chi^{2}_{\rm rw}$ are before and after reweighting, respectively; $N_{\rm eff}$ is the number of effective replicas; $\langle \alpha \rangle$ is the average rescaling of the covariance matrix.}
\end{table}

\begin{figure}
  \begin{centering}
    \includegraphics[scale=0.4]{whist-rw}\includegraphics[scale=0.4]{palpha-rw}
    \par\end{centering}
  \caption{\label{fig:weights} On the left, the distribution of
    weights after performing the reweighting with ATLAS $W$,
    $Z/\gamma^{*}$ rapidity data and ATLAS Drell-Yan high mass
    data. The prior set of $N_{\rm rep}=500$ replicas after the
    reweighting procedure is reduced to $N_{\rm eff}=266$ effective
    replicas. On the right, the associated $P(\alpha)$ distribution.}
\end{figure}

Figure~\ref{fig:pheno} shows a comparison between ATLAS data and
theoretical predictions computed with NNPDF2.3 NLO, with and without
photon-induced contributions using the photon PDF obtained after the
reweighting procedure. Again, the comparison shows that the
$Z/\gamma^{*}$ high mass production provides the most relevant
information to constrain the photon PDF.

\begin{figure}
  \begin{centering}
    \includegraphics[scale=0.39]{atlaswz-wpcomparison}\includegraphics[scale=0.39]{atlaswz-wmcomparison}
    \includegraphics[scale=0.39]{atlaswz-zcomparison}\includegraphics[scale=0.39]{atlashm-zcomparison}
    \par\end{centering}
  \caption{\label{fig:pheno} Comparison between experimental data and
    theoretical predictions computed with NNPDF2.3 NLO for each data
    set, with (blue lines) and without (green lines) photon-induced
    contributions.}
\end{figure}

\section{Unweighting the photon PDF}

Following the procedure presented in Ref.~\cite{Ball:2011gg} we
construct an unweighted photon PDF set. Figure~\ref{fig:rw} shows the
comparison between the photon PDF before and after performing the
reweighting with ATLAS data, in logarithmic scale (left) and linear
scale (right). We observe a significant reduction of the uncertainties
which essentially overlaps with the coverage in $x$ of the ATLAS data
used in the reweighting. It is interesting to observe that the central
value is also substantially reduced. This means that the central value
from the DIS fit is mostly determined by statistical fluctuations.

Figure~\ref{fig:photon} shows the details of the unweighted photon PDF
set with $N_{\rm rep}=100$ replicas. Also in this example, the current
reweighted photon PDF in comparison to the photon PDF from DIS has
small central values and smaller uncertainties at medium/small-$x$,
improving the quality of theoretical predictions. On the other hand,
we observe that at very small-$x$ the photon PDF has large
uncertainties, because the ATLAS data does not cover this region.

\begin{figure}
  \begin{centering}
    \includegraphics[scale=0.4]{photon_rw_log}\includegraphics[scale=0.4]{photon_rw}
   \par\end{centering}
 \caption{\label{fig:rw} Direct comparison between the photon PDF
   before (black lines) and after (green lines) the reweighting
   procedure. The reweighted photon PDF uncertainties are reduced at
   medium/small-$x$. The 68\% c.l. band is defined from the PDF mean value.}
\end{figure}

\begin{figure}
  \begin{centering}
    \includegraphics[scale=0.4]{qedcomparisonlog}\includegraphics[scale=0.4]{qedcomparison}
   \par\end{centering}
 \caption{\label{fig:photon} The unweighted photon PDF at
   $Q^{2}_{0}=2$ GeV$^{2}$ obtained after the reweighting procedure
   with ATLAS $W$, $Z/\gamma^{*}$ rapidity data and
   ATLAS Drell-Yan high mass data. The PDF set includes 100 replicas.}
\end{figure}

\section{Outlook}

In conclusion, the Drell-Yan LHC data from ATLAS reduces the photon
PDF uncertainties at medium/small-$x$, providing a much more
constrained and reliable photon PDF in comparison to the results
obtained by fitting DIS-only data.

Here we have shown only the impact of ATLAS data to the photon PDF. In
a the final photon PDF set~\cite{Carrazza:2013} we will also include
the LHCb low mass data, which reduces uncertainties at very
small-$x$. This will provide the first extraction of a modern photon
PDF with uncertainties and constrained by LHC data. Further
improvements of the photon PDF uncertainties will be possible when
more LHC data will be available, specifically, measurements of $WW$
production which at high energy scales constrains the photon PDF at
large-$x$.