text:
\section{introduction}

The investigation of quantum lattice models is one of the
central themes in modern condensed matter physics. A crucial
  step is to develop novel numerical methods to
efficiently simulate the interesting and complex phenomena of quantum
many-body systems. In particular, the tensor network states and the
related renormalization group methods, including the tree tensor
network state (TTN),\cite{Shi, Tagliacozzo, Murg} the multi-scale
entanglement renormalization ansatz,\cite{MERA} the projected
entangled pair state,\cite{Verstraete} the tensor renormalization
group (TRG),\cite{Levin, Xiang, Gu} and the second renormalization
group (SRG),\cite{Xie, Zhao} are now under rapid development. These
methods provide promising numerical tools for studying strongly
correlated systems, especially for the frustrated magnetic systems and
fermion models, and can be regarded as an extension of the fruitful
density matrix renormalization group (DMRG)\cite{White} in two or
higher dimensions.

In the study of the tensor network methods, one needs to first
determine the tensor network wavefunction for the ground state.  In
Refs. [\onlinecite{Xiang, Zhao}], a simple update scheme is proposed
to determine the ground state tensor network wavefunction in two
dimensions. This scheme is efficient and robust.  It proceeds in three steps:
(1) apply the imaginary time projection operators
simultaneously on bonds of the same type, for
example the $x$ directional bonds in Fig. \ref{fig-bethe-latt}a,
and enlarge the bond dimension; (2) construct a local evolving
block matrix and simulate the environment contribution by the diagonal
matrices on the external bonds [$\lambda_y$ and $\gamma_z$ in
Fig. \ref{fig-bethe-latt}b]; (3) decompose the evolving block
matrix by singular value decomposition (SVD) and decimate the vector
space of the enlarged geometric bond according to the singular values
in the updated diagonal matrix $\theta'_x$.  This technique has been
combined with the TRG/SRG to evaluate the ground state properties of
two-dimensional (2D) Heisenberg models.\cite{Zhao, Gu, Chen, Li,
  Zhao2} It is an accurate numerical method for evaluating local
physical quantities, but it is less accurate in evaluating the
long-range correlation functions.\cite{Zhao} This is the major
drawback of this simple update scheme.  It results from a mean-field
approximation for the environment tensor.  A way to go beyond this
approximation is to enlarge the size of the cluster that is used for
evaluating the environment tensor. This, as shown by Wang and
Verstraete\cite{Wang}, can indeed improve the accuracy for the long
range correlation function.

In this work, we apply the simple update scheme to infinite
TTN (iTTN) states on the Bethe lattice. We will show that
this is a quasi-canonical approach for treating an iTTN. Here by the
word ``quasi-canonical'' we mean that with increasing the number of iteration
steps and decreasing the Trotter error, the tree tensor network state
obtained by the simple update scheme would become asymptotically
canonical [i.e.\ the tensors satisfy certain canonical
  orthonormality conditions, see \Eqs{eq-can3} below]. Thus the
simple update scheme provides an accurate and efficient approach for
evaluating the ground state wavefunction on the Bethe lattice.

The Bethe lattice, as shown in Fig. \ref{fig-bethe-latt}a, has a
self-similar structure with an infinite Hausdorff dimension. The size
of the lattice is infinite, hence the boundary effects do not need to be explicitly
considered. The Bethe lattice was first used in the study of classical
statistical mechanics.\cite{Bethe, Plischke, Baxter} It has attracted
broader interest since a number of chemical compounds with the Bethe
lattice structures, such as the dendrimers,\cite{Delgado} have been
synthesized in the laboratory.\cite{Astruc}

A finite Bethe lattice is called a Cayley tree.  Soon after White's
invention of DMRG,\cite{White} the DMRG algorithm
for the quantum lattice models defined on the Cayley tree was
proposed. \cite{Otsuka, Friedman} Based on the DMRG calculation of
local physical quantities in the central part of the Cayley tree,
Otsuka\cite{Otsuka} claimed that the anisotropic S=1/2 Heisenberg
model (i.e. the XXZ model) on the Bethe lattice should exhibit a
first-order phase transition at the isotropic point. Later
Friedman\cite{Friedman} proposed an improved DMRG scheme and evaluated
the spin-spin correlations in the ground state. Based on the DMRG
result, he suggested that long-range magnetic order might exist at
the isotropic Heisenberg point. Recently, Kumar et al. calculated the
magnetization with a further improved DMRG algorithm, \cite{Kumar} and
showed that such long-range magnetic order does exist at that point.

The above DMRG calculations were done on the Cayley tree lattice, not
on the true infinite Bethe lattice.  Furthermore, it should be pointed out
that the boundary effect is very strong on a finite Cayley tree since
more than one-half of the total sites reside on the lattice edge. This
may strongly affect the properties of the system. In some cases, the
results obtained on a Cayley tree lattice can be completely different from
those for the corresponding Bethe lattice. For example,
the classical Ising model shows a phase transition on the Bethe
lattice, but not on the Cayley tree lattice.\cite{Ostilli}

To unambiguously resolve the above problems, it
is necessary to calculate the spin models directly on the Bethe
lattice.  The recent development of the TTN algorithms \cite{Shi,
  Tagliacozzo, Murg} has indeed made this feasible.\cite{Nagaj, Nagy}
In particular, Nagaj \textit{et al} in Ref. \onlinecite{Nagaj}
extended the infinite time evolving-block decimation\cite{Vidal}
technique to the Bethe lattice and determined the ground state
wavefunction by imaginary time evolution. For the transverse Ising
model on the Bethe lattice, it was found that a second-order quantum
phase transition exists at a critical transverse field. An interesting
result revealed in this calculation is that even at the second-order
critical point, the correlation length remains finite. For the Bethe
lattice with coordination $q=3$, the correlation length is shown to be
less than $1/\ln2$.  However, in the calculation Nagaj \textit{et al}
used a three-site projection operator to simultaneously evolve the two
equivalent incoming legs of the tensors, the computational cost is
thus very high. The computational time scales as O($D^8$) with $D$ the
tensor dimension, which limits the value of $D$ that can be handled to
$D\le 8$.

Recently, Nagy\cite{Nagy} proposed a different algorithm to reduce the computational cost by making use of the $C_3$ rotational symmetry of $q=3$ Bethe lattice. This algorithm reduces the computational cost to O($D^4$) hence greatly improves the efficiency. It can be used for studying the spin-1/2 quantum lattice models. However, the application of this method is restricted to the translation invariant spin-1/2 system.

\begin{figure}[tbp]
\includegraphics[angle=0,width=1.0\linewidth]{fig1.eps}
\caption{(Color online) (a) The $q=3$ Bethe lattice. Every site has 3 nearest neighbors, and the three bonds are labeled according to their directions as $x$, $y$, and $z$, respectively.  (b) The two-site cluster used in the single-bond projection of the simple update scheme. The diagonal matrices $\lambda_y$ and $\gamma_z$ on the dangling bonds should be included in the projection to mimic the entanglement renormalization of the environment to this two-site system. (c) A minimum cluster that is used in the Bethe approximation. It consists of one A tensor and three nearest-neighbor B tensors (or vice versa).
} \label{fig-bethe-latt}
\end{figure}

As will be shown below, the simple update scheme is very efficient. Its computational costs scale as O($D^4$), similar as for the algorithm proposed by Nagy\cite{Nagy}. But it is much more flexible. It can be applied to treat arbitrary TTN states, with or without translation invariance. Here we studied two spin models defined on the Bethe lattice. One is the transverse Ising model and the other is the antiferromagnetic XXZ Heisenberg model. The quantum phase transitions and the ground state phase diagrams of these models are studied.

The rest of the paper is arranged as follows. An introduction to the
simple update scheme and its relationship with the Bethe approximation
is presented in Sec. \ref{sec:2}. The study of the quantum phase
transitions of transverse Ising and XXZ Heisenberg models are presented
in Secs. III and IV, respectively. In Sec. V, the present scheme is
generalized to larger tree tensor clusters, in order to provide more
accurate approximations for 2D lattices. Finally, Sec VI is devoted to a summary.

\section{The canonical form and the simple update scheme}
\label{sec:2}

The iTTN state on the Bethe lattice comprises 4-indexed tensors $A_{x,y,z}^m$ and $B_{x,y,z}^m$ defined on the vertices, and the diagonal matrices $\theta$, $\lambda$, $\gamma$ defined on the bonds along the $x$, $y$, and $z$ directions as shown in Fig. \ref{fig-bethe-latt}a, respectively.
The bond indices represent the quantum numbers of the virtual basis states. The physical index $m$ runs over the $d$ basis states of the local Hilbert space at each lattice site. The diagonal matrices store the entanglement information, and play an important role in the simple update scheme.

In order to determine the ground state wavefunction, the imaginary time evolving operators $U(\tau) = \exp{(-\tau h_{i,j})}$ are applied to the iTTN iteratively.
At each step, the dimension of the evolved bond is increased by a factor of $d^2$.
Thus the tensor dimensions will proliferate exponentially with the increasing number of projection steps.
In order to sustain the projections until the iTTN converges to the true ground state wavefunction, one needs to truncate the bond dimension after each projection step. This needs a proper consideration of the renormalization effect of the environment tensor.

An accurate and full determination of the environment tensor is
computationally costly.  This limits generally the tensor dimension
$D$ that can be handled to a rather small value, say $D \le 6$.  The
simple update scheme\cite{Xiang, Xie, Zhao}, on the other hand, takes
the product of the dangling bond matrices as a mean field
approximation to the environment tensor.  It converts the complicated
global optimization problem into a local one, and
hence greatly simplifies the calculation.  On the
regular 2D lattice, the bond matrix is an approximate measure of the
entanglement between the system and environment tensors.  However, as
will be shown later, the square of the diagonal bond matrix on the
Bethe lattice is the eigenvalue of the reduced density matrix if the
iTTN is canonicalized, i.e., if the tensors in the network are
always kept in canonical form by some transformations.
Thus the simple update scheme is an accurate treatment for
the renormalization of the iTTN on the Bethe
lattice.

The simple update scheme is also closely related to the famous Bethe
approximation.\cite{Bethe, Plischke, Baxter} To understand this, we
show in Fig. \ref{fig-bethe-latt}c a 4-site cluster, which contains
one A tensor and three B tensors.  In the simple update calculation,
the two local tensors, A and B, and the three inner bond matrices
($\theta_x$, $\lambda_y$, $\gamma_z$) should be evaluated and updated
iteratively. After each single projection step on the inner bonds, to
keep the scheme self-consistent, one should also update all the
dangling bonds of the cluster, by replacing the bond matrices with the
corresponding ones on the inner bonds.

This cluster structure and the self-consistent scheme is in fact the
Bethe approximation that was first proposed by Bethe in 1930's, in the
context of statistical mechanics.\cite{Bethe} The key idea is to treat
the correlations between the central spin and its nearest neighbors in
the cluster exactly, and to use an effective mean field to approximate
the interactions between the cluster and the rest lattice spins.  By
solving this simple cluster problem, and assuming that all the spins
in the cluster have exactly the same local magnetization, one can
determine the spontaneous magnetization self-consistently.  For the
quantum cases, the six diagonal matrices on the dangling bonds of the
cluster are taken as the mean fields acting on the inner block. The
self-consistent condition requires that the matrices $\theta$,
$\lambda$, and $\gamma$ on dangling bonds are equal to the corresponding
matrices on the inner bonds between A and B tensors (see
Fig. \ref{fig-bethe-latt}c).

\begin{figure}[tbp]
\includegraphics[angle=0,width=1.0\linewidth]{fig2.eps}
\caption{(Color online) One iteration step in the simple update scheme:
(a) A and B tensors are connected by the bond $x$, which will be involved in the following projection steps.
There are diagonal matrices on the dangling bonds ($y$ and $z$ bonds) of A and B tensors.
(b) Absorb the four dangling matrices into A and B, and define the block matrix $M_{a(b)}$.
Then take the QR decomposition for $M_{a(b)}$, obtaining $Q_{a(b)}$ and $R_{a(b)}$ matrices.
(c) Project $U(\tau)$ onto the bond by contractions, and obtain the block matrix $G$ [see Eq. (\ref{eq-G})].
(d) Take singular value decomposition of $G$ to find the unitary matrices $U$ and $V^{T}$,
and the new diagonal matrix $\theta'$.
(e) Truncate the $x$-bond dimension to $D$ according to the diagonal values of $\theta'$.
Merge $U$ ($V^T$), $\gamma_y^{-1}$, and $\lambda_z^{-1}$ together into $Q_{a(b)}$,
finally we arrive at the updated $A'(B')$ tensors.}
\label{fig-CTN}
\end{figure}

A tensor network state contains redundant gauge degrees of freedom on
each bond.  It is invariant if one inserts a product of two reciprocal
matrices on a bond and absorbs separately each of them to a local
tensor at the two ends of the bond.  This gauge invariance of a tensor
network state can be used to simplify the calculation of local
tensors, especially for the iTTN states on the Bethe lattice, where a
special gauge, called canonical form, can be introduced.

To be specific, the local tensors of canonical iTTN states satisfy the
following orthonormality conditions
\begin{eqnarray}
\sum_m
\sum_{x,y} \theta_x^2 \lambda_y^2 (T^{m}_{x,y,z'})^* T^{m}_{x,y,z}
&=& \delta_{z',z},
\label{eq-can1}
\\
\sum_m
\sum_{y,z} \lambda_y^2 \gamma_z^2 (T^{m}_{x',y,z})^* T^{m}_{x,y,z}
&=& \delta_{x',x}, \\
\sum_m
\sum_{z,x} \gamma_z^2 \theta_x^2 (T^{m}_{x,y',z})^* T^{m}_{x,y,z}
&=& \delta_{y',y},
\label{eq-can3}
\end{eqnarray}
where T represents the A or B tensor. If we cut an arbitrary bond to
divide the Bethe lattice into two parts, denoted as a system and an
environment subblock, one can then define the reduced density matrix
of the system block by integrating out all the degrees of freedom
of the environment block. For the tensors that satisfy
Eqs. (\ref{eq-can1}-\ref{eq-can3}), the square of the diagonal bond matrices are the eigenvalues, and
the renormalized bond bases are the eigenvectors of the corresponding
reduced density matrix, which are orthonormal to each other. Thus, in
terms of the Schmidt decomposition, the square of the diagonal matrix
elements represent the probability amplitudes of the corresponding
eigenvectors appearing in the wavefunction.

The existence of this simple canonical form of the iTTN,
i.e. Eqs. (\ref{eq-can1}-\ref{eq-can3}), is very useful in the
calculations.  First, the diagonal bond matrix describes the
entanglement spectrum between the system and environment subblocks.
Thus to select the virtual bond basis states according to the values
of these diagonal matrix elements provides an optimal scheme to
truncate the bond dimension.  Second, the contribution of the
environment tensors can be faithfully represented by the 4 diagonal
matrices on the dangling bonds surrounding the central bond under
projection (see Fig. \ref{fig-CTN}a).  It means that the imaginary
time evolution on each bond can be done rigorously and locally.
Furthermore, we can also evaluate the expectation value of a local
operator simply by contracting a small cluster comprising those
tensors and bond matrices on which the operator acts.
This significantly reduces the computational cost.

Bearing in mind the benefits of the canonical iTTN states, one
can perform explicitly the canonical transformations during the projection
processes. However, to further save computational costs,
in practical calculations we choose to carry out the
imaginary time evolution just using the simple update scheme,
and gradually reduce the Trotter step $\tau$, which would bring
the iTTN states into its canonical form step by step.
This scheme works because the diagonal bond matrix provides an
approximate measure for the entanglement between the two sides of the
bond and can be used to substitute approximately the environment
tensor.  Therefore, it can stabilize the algorithm of the imaginary time
evolution, provided that the Trotter step $\tau$ is small enough so
that the bond projection operator $U(r)$ is nearly
unitary.\cite{Vidal} This near unitary evolution can
modify the wavefunction and reshape it in order to satisfy the canonical conditions.
The simple update scheme hence provides a quasi-canonical
evolution of the iTTN state, which will finally converge to the ground state
and become canonical in the limit $\tau\rightarrow 0$. In practical calculations,
the Trotter step $\tau$ is gradually decreased from $10^{-1}$ to $10^{-4}$, and the total number
of projections steps varies from 20000 to 200000.

Now let us consider how to implement the simple update scheme
efficiently.  A simple approach is to do directly the singular value
decomposition the evolving block tensor, which is a matrix of $D^2
\times D^2$.  The computational cost for doing this singular value
decomposition is high, since it scales as
$O(D^6)$. This cost can in fact be reduced to $O(D^4)$ if we carry
out this singular value decomposition in the following steps (again,
projection on $x$ bond is taken as an example):

(1) Define the following two $D^2 \times D d$ block
matrices (Fig. \ref{fig-CTN}b)
\begin{eqnarray}
(M_a)_{y,z; x,m} & = & \lambda_y \gamma_z A_{xyz}^m, \\
(M_b)_{y,z; x,m} & = & \lambda_y \gamma_z B_{xyz}^m,
\end{eqnarray}
by absorbing the diagonal matrices $\lambda_y$ and $\gamma_z$ into the tensors A and B, and calculate their QR decomposition
\begin{equation}
(M_\alpha)_{y,z; x, m} =  \sum_{k} Q^\alpha_{y,z; k} R^\alpha_{k;m,x},
\end{equation}
where $\alpha = a$ or $b$. $Q^\alpha$ is a $D^2 \times Dd$ column orthonormal matrix. $R^\alpha$ is a $Dd \times D d$ upper diagonal matrix.

(2) Apply the bond projection operator $U(\tau)$ to
the system and define the gate matrix (Fig. \ref{fig-CTN}c)
\begin{equation}
G_{m_1k_1;m_2k_2} = \sum_{x, m_1^\prime, m^\prime_2}
\langle m_1 m_2 | U(\tau ) |m_1^\prime m^\prime_2\rangle
R^a_{k_1;m_1^\prime x} \theta_x R^b_{k_2;m_2^\prime x}
\label{eq-G}
\end{equation}
(3) Take the singular value decomposition for this matrix (Fig. \ref{fig-CTN}d),
\begin{equation}
G_{m_1k_1;m_2k_2} = U_{m_1k_1;l} \theta^\prime_l V_{m_2k_2; l},
\end{equation}
where $U$ and $V$ are two $Dd\times Dd$ unitary matrices, and
$\theta^\prime$ is a semi-positive defined matrix.

(4) Truncate the inner bond
dimension by keeping the largest $D$ matrix elements of
$\theta^\prime$, and update the local tensors by
the formula (Fig. \ref{fig-CTN}e)
\begin{eqnarray}
{A^\prime}^m_{xyz} &=& \sum_k  \lambda_y^{-1} \gamma_z^{-1} Q^a_{y,z;k} U_{m,k;x}, \\
{B^\prime}^m_{xyz} &=& \sum_k  \lambda_y^{-1} \gamma_z^{-1} Q^b_{y,z;k} V_{m,k;x}.
\end{eqnarray}

Combining this efficient simple update scheme and the local
determination of physical observables using the canonical form, we can
keep the computational cost in a low level. In practice, this
allows us to keep a relative large bond dimension.

\begin{figure}[tbp]
\includegraphics[angle=0, width=0.85\linewidth]{fig3.eps}
\caption{(Color online)  (a) The longitudinal and transverse magnetizations $m_z$ and $m_x$ versus the transverse fields $h_x$. The transverse magnetization $m_x$ increases monotonously with $h_x$, while the longitudinal magnetization $m_z$ decreases and vanishes at the transition point. (b) The second-order derivative of the ground state energy per site $e$ with respect to $h_x$, $d^2 e/ dh_x^2$, obtained by taking the first-order derivative of $m_x$, $dm_x/dh_x$.}
\label{fig-magn-q-Ising}
\end{figure}

\section{The transverse Ising model}

The transverse Ising model is defined by the Hamiltonian
\begin{equation}
H_{\rm{TI}} = - \sum_{<i,j>} J S_i^z S_j^z - \sum_{i} h_x S_i^x - \sum_{i} h_z S_i^z,
\end{equation}
where the spin-spin exchange constant $J$ is set as the energy scale ($J=1$, ferromagnetic coupling).
The second term represents the transverse-field along $S^x$ direction, and the last term is the longitudinal-field along $S^z$ direction.

\begin{figure}[tbp]
\includegraphics[angle=0, width=0.85\linewidth]{fig4.eps}
\caption{(Color online) (a) The entanglement entropy $S_{\rm{E}}$ of the ground state for the transverse Ising model. The cusp at $h_c \simeq 1.115$ corresponds to the second-order phase transition point. (b) The correlation length $\xi$ of the ground state, which also shows a cusp at the transition point.}
\label{fig-ent-q-Ising}
\end{figure}

Figure \ref{fig-magn-q-Ising}a shows the longitudinal and
transverse magnetizations, $m_z = \langle S_z \rangle$ and $m_x =
\langle S_x \rangle$, as a function of the transverse field $h_x$. A
continuous order-disorder phase transition is found at $h_c$. For $h_x
< h_c$, the ground states undergo a spontaneous $Z_2$ symmetry
breaking with a finite longitudinal magnetization $m_z$, which
decreases with increasing $h_x$ and vanishes at the critical field. By
utilizing the Hellmann-Feynman theorem, the second-order derivative of
the ground state energy can be calculated by $d^2e/dh_x^2 =
-dm_x/dh_x$. As shown in Fig. \ref{fig-magn-q-Ising}b, $d^2e/dh_x^2$
exhibits a discontinuity at $h_c$, indicating that $h_c$ is a
second-order phase transition point.  The critical field is found to
be $h_c \simeq 1.115$, in agreement with previous
calculations.\cite{Nagaj, Nagy} It is also close to the critical field
$h_c = 1.06625(2)$ for the transverse Ising model on the honeycomb
lattice.\cite{Deng}

\begin{figure}[tbp]
\includegraphics[angle=0, width=0.9\linewidth]{fig5.eps}
\caption{(Color online) The field dependence of the longitudinal susceptibility $\chi^z$. It shows a divergent peak at the transition point $h_c \simeq 1.116$. The susceptibility is calculated by $\chi^z = [m_z(h_z)-m(0)]/\delta h$ with $h_z = 10^{-4}$.} \label{fig-sus-q-Ising}
\end{figure}

Figure \ref{fig-ent-q-Ising} shows the bipartite entanglement entropy $S_{\rm{E}}$ \begin{equation}
S_{\rm{E}} = - \rm{Tr}[\Lambda^2 \log_2(\Lambda^2)],
\label{eq:entropy}
\end{equation}and the correlation length $\xi$ for the ground state. In Eq. (\ref{eq:entropy}), $\Lambda $ = $\theta$, $\lambda$ or $\gamma$ is a diagonal matrix that satisfies the canonical condition. $S_{\rm{E}}$ shown in Fig. \ref{fig-ent-q-Ising}a is obtained by taking the average over the three bonds.

The correlation length $\xi$ is evaluated from the ratio of the the largest ($a_0$) and the second largest eigenvalue ($a_1$) of the transfer-matrix for the iTTN state,
\begin{equation}
\xi = 1/\ln{\frac{a_0}{a_1}} \; .
\end{equation}
For the $q=3$ Bethe lattice, there are six kinds of transfer-matrices, depending on
the site and the bond direction. For example the transfer matrix along
the $yz$-direction is defined by (for A sublattice site)
\begin{equation}
T^{a}_{y,z; y',z'}  =  \sum_{m,x} \sqrt{\lambda_y \gamma_z} (A^m_{x,y,z})^* \theta^2_x A^m_{x,y',z'} \sqrt{\lambda_{y'} \gamma_{z'}}.
\end{equation}
The other five transfer matrices including $T^b_{y,z;y',z'}$, $T^{a,b}_{z,x; z',x'}$, and $T^{a,b}_{x,y; x',y'}$ can be similarly defined. The results of the correlation length shown in Fig. \ref{fig-ent-q-Ising}b are evaluated from the product of the six transfer matrices along a specific path.

A distinctive feature revealed by Fig. \ref{fig-ent-q-Ising} is that
the correlation length $\xi$, as well as the entanglement
$S_{\rm{E}}$, does not diverge at the critical point. $\xi$ is found
to be upper bounded by $1/\ln{2}$, in agreement with the published
results\cite{Nagaj, Nagy}. This peculiar behavior is not observed in
the ordinary continuous phase transition systems, where the
correlation length is always divergent at the critical point.

We will now show that this non-critical behavior of the correlation
length at the critical point is due to the peculiar geometry of the Bethe lattice.
For the Bethe lattice, the number of sites on the boundary of
a finite connected region is roughly equal to the number of
internal sites within that region. It means that the lattice sites are highly
non-uniformly distributed as a function of lattice distance away from
a given center. This is a feature of the Bethe lattice that differs
from a regular lattice.

In order to understand why the correlation length is finite at the
critical point, let us take a scaling transformation to convert the
Bethe lattice to a ``regular'' 2D lattice whose lattice sites are
uniformly distributed in space. To do this, we first choose an arbitrary site,
to be viewed as ``center'' of the lattice, and define the distance $R$ for a given layer $r$ to the
center as
\begin{equation}
R \propto \sqrt{ \frac{N(r) }{ \pi } } \sim \sqrt{ \frac{(q-1)^r }{ \pi} }
\end{equation}
where $N(r) \propto (q-1)^r$ is the number of sites enclosed by the
$r$ layer.  In this rescaled lattice, $ r \sim 2\ln{R} + \mathrm{const}
$, and the exponentially decaying correlation
function $C(r) \sim \exp(-r/\xi)$ in the original Bethe lattice
corresponds to a power-law decaying function of $R$
\begin{equation}
C( R ) \sim R^{-2/\xi}.
\end{equation}
The algebraic decay of this correlation function
suggests the spins on the Bethe lattice are actually long range
correlated in the rescaled framework, even away from the critical
point. This is the reason why the system can undergo a phase
transition without exhibiting a divergent correlation length
at the critical point in the original Bethe lattice.

To understand why the correlation length is upper bounded, let us
consider the longitudinal magnetic susceptibility $\chi_z$ of the ground state (Fig. \ref{fig-sus-q-Ising})
\begin{equation}
\chi_z = \frac{d m_z}{d h_z} = \sum_{i} \langle S_0^z S_i^z \rangle_{\rm{GS}} - \langle S_0^z \rangle_{\rm{GS}} \langle S_i^z \rangle_{\rm{GS}},
\end{equation}
where $S_0^z$ is the spin at a reference center, $i$ runs over all the sites on the lattice.
$\langle \hat{O} \rangle_{\rm{GS}}$ is the expectation value of operator
$\hat{O}$ in the ground state. Exploiting the $C_3$ rotational symmetry
of the Bethe lattice, we define $C^z(r) = \langle S_0^z S_i^z \rangle_{\rm{GS}} - \langle S_0^z
\rangle_{\rm{GS}} \langle S_i^z \rangle_{\rm{GS}}$, with $r$ the layer number
where site $i$ resides. Thus $\chi_z$ can be rewritten as
\begin{equation}
\chi_z = \sum_r n(r) C^z(r),
\label{eq-chiz}
\end{equation}
where $n(r)$ is the number of spins on layer $r$.  On a regular
lattice, $n(r) \propto r^{\nu-1}$, where $\nu$ is the spatial
dimension of the lattice, the susceptibility is always finite if the
spin-spin correlation function $C^z(r)$ decays exponentially.
However, in the Bethe lattice, $n(r) \propto (q-1)^{r-1}$. Now if we
assume $C^{z}(r) \propto \exp(-r/\xi )$, then
\begin{equation}
\chi_z \propto \sum_r (q-1)^r e^{-r / \xi} = \sum_r e^{r \left[\ln (q-1) - 1/\xi\right]},
\end{equation}
which diverges if $\xi$ approaches the threshold value $1/\ln(q-1)$. This shows that the
susceptibility can diverge even if $C^{z}(r)$ decays exponentially
with $r$. The critical point occurs when $\xi = 1/\ln (q-1)$,
and the correlation length $\xi$ is therefore upper-bounded by
$1/\ln (q-1)$ on the Bethe lattice. \cite{Semerjian}

\begin{figure}[tbp]
\includegraphics[angle=0, width=0.85\linewidth]{fig6.eps}
\caption{(Color online) (a) The ground state energy and (b) the bipartite entanglement entropy as a function of $\delta$ for the XXZ model. The appearance of the cusp in the energy,\Y{as well as the discontinuity in the first-order energy derivative [inset in (a)]}, suggests that this is a first-order phase transition point.}
\label{fig-1st-QPT-XXZ}
\end{figure}

\begin{figure}[tbp]
\includegraphics[angle=0, width=0.9\linewidth]{fig7.eps}
\caption{(Color online) The staggered magnetization along the z-axis ($m_z^s$) and $x$-axis ($m_x^s$). The spin orientations flip suddenly at the transition point $\delta=1$. For $\delta < 1$, the spins are ordered within the $xy$-plain, while for $\delta > 1$, the ordering is along the $z$-axis.}
\label{fig-magn-XXZ}
\end{figure}

\section{Anisotropic Heisenberg model}

The anistropic Heisenberg model, i.e. the XXZ model, is defined by the Hamiltonian
\begin{equation}
H_{\rm{XXZ}} = \sum_{<i,j>} (S_i^x S_j^x + S_i^y S_j^y + \delta S_i^z S_j^z),
\end{equation}
where $\delta$ is the anisotropy parameter.

The above model has been intensively studied on the honeycomb and
square lattices by different numerical methods, which include
exact diagonalization,\cite{Lin} quantum Monte
Carlo,\cite{Lin, Barnes} coupled cluster
methods,\cite{Bishop1,Bishop2} and tensor network
algorithms.\cite{Chen, Bauer} It is found that the system possesses
magnetic long-range orders for all values of $\delta$.  The
antiferromagnetic ordering vector points within the easy $xy$-plane
for $\delta < 1$ or along the $z$-axis for $\delta > 1$. There is a
first-order phase transition at $\delta = 1$, the Heisenberg
point.\cite{Hajj}

This model was also studied on the Bethe lattice (more precisely, on the Cayley tree lattice)
by DMRG,\cite{Otsuka, Friedman, Kumar, Changlani}.
It was found that there exists a long-range magnetic order
at the isotropic point $\delta=1$. It was also suggested that a quantum phase
transition occurs at this point. However, the properties of this
transition and the phases on the two sides of the critical point have
not been clarified.

Figure \ref{fig-1st-QPT-XXZ} shows the $\delta$-dependence of the
ground state energy per bond and the entanglement entropy for
the XXZ model.  A clear first order quantum phase transition is
observed at $\delta=1$.  The energy per bond shows a change of slope at $\delta=1$ (\Y{the first-order
energy derivative is shown in the inset}), which suggests
that there is an energy level crossing. The entanglement entropy
varies continuously across the transition point, but exhibits a
cusp.

Figure \ref{fig-magn-XXZ} shows the staggered magnetization $m_x^s$
and $m_z^s$ around the critical point. The ground state is found to
possess in-plane antiferromagnetic order with a finite
$m^s_x$ for $\delta < 1$, and $z$-axis
antiferromagnetic order with a finite $m^s_z$ for $\delta > 1$. At the transition point,
the two order parameters change suddenly, displaying a spin flip
transition. This result verifies the conjecture made by
Otsuka. \cite{Otsuka}

At the transition point, we find that the ground state energy per bond
has the value $e_b = -0.359817(3)$ and the spontaneous
magnetization has the value $m_s = 0.34736(1)$ for $D=40$. The
errors in the parentheses are estimated by comparing the results for
different bond dimension $D$ and different Trotter slices $\tau$. Our
results agree well with the DMRG data published in
Ref. [\onlinecite{Kumar}], where the local magnetization is found to
be $m=0.347$ on the central lattice site and the bond energy between
the central spin and a spin on the first layer is $e=-0.359$. This satisfactory agreement
suggests that by calculating the Bethe lattice, we can reproduce the
results of local properties in the very center of a large Cayley tree.

\begin{figure}[tbp]
\includegraphics[angle=0, width=1.0\linewidth]{fig8.eps}
\caption{(Color online) The one-ring cluster with 12 sites. The diagonal matrices $\Lambda$ are defined on the dangling bonds. The dashed lines represent removed bonds (the physical couplings on the dashed line still exist).}
\label{fig-one-ring}
\end{figure}

\begin{figure}[tbp]
\includegraphics[angle=0, width=1.0\linewidth]{fig9.eps}
\caption{(Color online) The 4-ring cluster with 26 sites. The inequivalent lattice sites are numbered from 1 to 16. $\Lambda_{i,j}$ labels the diagonal matrix on the bond linking sites $i$ and $j$.}
\label{fig-four-ring}
\end{figure}

\begin{figure}[tbp]
\includegraphics[angle=0, width=1.0\linewidth]{fig10.eps}
\caption{(Color online) The ground state energy per site for the Heisenberg model on the honeycomb lattice. The results are evaluated in the central area of the clusters. The dashed line is the recent quantum Monte Carlo result.\cite{Low} The number of hexagonal rings is used to label the cluster size.}
\label{fig-eng-clust}
\end{figure}

\section{Cluster update scheme}

In the previous sections, the simple update has been applied to study the quantum spin models on the Bethe lattice, leading to very accurate results. What is more, in terms of the Bethe approximation, these results can also be regarded as approximations for the corresponding 2D lattice models. Actually, the simple update scheme has already been used to study regular 2D lattices, such as the honeycomb or square lattices. Combined with the TRG/SRG techniques, it can achieve rather accurate results.\cite{Xiang,Xie,Zhao} Nevertheless, in this section, we would provide a different way of using the simple update to calculate 2D lattices. Inspired by the generalization of the Bethe approximation to larger clusters in classical statistical mechanics,\cite{Plischke} we apply simple update to various tree tensor clusters. In this way, the advantage (efficiency) in treating a tree tensor network, namely the fact that it can be readily canonicalized, is utilized to improve the calculation accuracy on a regular 2D lattice.

To start, as a first-order approximation, let
us compare the results on the $q=3$ Bethe lattice (which has no loops)
with those on the 2D honeycomb lattice (whose coordination number is also $q=3$, and it does have loops).
Our result for the ground state energy of the Heisenberg
model on the $q=3$ Bethe lattice is $e_b = -0.359817(3)$, while the
corresponding energy on the honeycomb lattice obtained by the recent
quantum Monte Carlo calculation is
$e_{\rm{QMC}}=-0.36303(14)$. \cite{Low} The relative difference
between these two energies is less than $0.9\%$.  However, the
spontaneous magnetization for the ground state of the Heisenberg model
on the Bethe lattice, namely $m_s = 0.34736(1)$, is much larger
than the corresponding value on the honeycomb lattice, which is about
0.27 as obtained by the quantum Monte Carlo.\cite{Low} Notice, some
other results for the magnetization on the honeycomb
or square lattice obtained with the tensor network algorithms are also
found to be higher than the Monte Carlo ones.\cite{Xiang,Bauer}

As a next step, our approximate treatment of the 2D honeycomb
lattice can be improved by using tensor networks that include rings.
In Fig. \ref{fig-one-ring}, the cluster with one hexagonal ring is
shown, some geometric bonds are removed (dashed lines in
Fig. \ref{fig-one-ring}) to form a tree tensor cluster. Note although
the tensor network does not have geometric
bonds on the dashed lines, in the Hamiltonian the couplings along
these bonds nevertheless exist. Therefore, the projections by
imaginary time evolution should be executed also on the dashed lines.
This cannot be done directly as on usual bonds, but can be accomplished as follows with the help of the swap gates.
The swap gates are used to exchange the physical indices of two tensors,
which proceeds similarly as the projection scheme illustrated in Fig. \ref{fig-CTN},
with the minor revision that the imaginary tim evolving operator $U(\tau)$
is now replaced with a swap operator $U_s$,
that conducts $U_s | m_i, m_j \rangle = | m_j, m_i \rangle$.

In Fig. \ref{fig-one-ring}, take the dashed bond between site A and F as
an example, swap gates moves the physical index on site A in the order
A $\to$ B $\to$ C, and the physical index on site F as F $\to$ E $\to$
D. After that, the two spins are linked by the solid bond between C and D,
then we can take the projection and update processes as on an usual
bond. After that, we have to move the two spin indices back to their
original positions by reversed swap operations, which accomplishes the
special projection step on a dashed bond. Through iterative and
self-consistent projection processes on the solid and the dashed bonds
of the tree cluster, an approximation for 2D lattices can be obtained.
Compared with the simple Bethe lattice, this tree tensor cluster
approach can provide better approximation for 2D. On the other hand,
it can also be regarded as an ideal method for evaluating the ``super Bethe
lattice'', of which each ``single site'' is now placed with a
hexagonal ring, instead of a single site, and the coordinate number $q=6$.

Beyond the one-ring cluster, more rings can be included to further
enlarge the cluster. As an example, Fig. \ref{fig-four-ring} shows a
cluster with 4 hexagonal rings. The accuracy of energy versus
different cluster size (labelled by the number of rings included) are
shown in Fig. \ref{fig-eng-clust}, which verifies that the accuracy
could be improved consistently with enhancing the cluster size. To
obtain better approximation for true 2D lattices, the local
observables are detected in the center area of the cluster. In
practice, for the 4-ring tree tensor cluster in
Fig. \ref{fig-four-ring}, the results are obtained by averaging over sites 3 and 4.
We find an energy per site of $e \simeq
-0.54441$ (bond energy $e_b \simeq -0.36294$) and a local
magnetization of $m = [e(h_s)-e(h_s=0)]/h_s \simeq 0.3147$ (with a staggered magnetic field $h_s=0.01$).
Hence, the inclusion of rings clearly improves the agreement with
QMC data. For the transverse Ising model, through the 4-ring cluster calculations, the phase transition
point is estimated as $h_c \simeq 1.1$, which is also more accurate
than the simple Bethe approximation. More numerical results with
larger clusters and further details of this cluster Bethe
approximation will be published separately.

\section{Conclusion}

In summary, the simple update scheme is employed to study two spin
models on the Bethe lattice, i.e., the transverse Ising and the Heisenberg
XXZ model. For the Ising model, it is found that the correlation
length, as well as the entanglement entropy, does not diverge at the
second-order transition point. Through a scale transformation, we
have given an intuitive explanation of this peculiar ``critical''
phenomenon. Moreover, by studying the magnetic susceptibility,
we show that the correlation length is upper bounded.
For the Heisenberg XXZ model, the existence of a
first-order phase transition at the isotropic point is clearly
verified, and the two different magnetic
ordered phases are identified as the easy-plain and easy-axis phases,
respectively. Furthermore, in terms of the Bethe approximation, we
obtain accurate and scalable approximations for the 2D lattice
models by applying the simple update to tree tensor clusters.

\section{Acknowledgement}

The authors would like to thank Tomotoshi Nishino and Zhong-Chao Wei for stimulating discussions.
WL is also indebted to Gang Su, Cheng Guo, Guang-Hua Liu, Ming-Pu Qin, Li-Ping Yang,
and Hui-Hai Zhao for helpful discussions. TX was supported by the National Natural
Science Foundation of China (Grants No. 10934008 and No. 10874215) and
the MOST 973 Project (Grant No.  2011CB309703). WL was supported
by the DFG through SFB-TR12.


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

Stripped-envelope (SE) core-collapse (CC) supernovae (SNe) mark the terminal endpoint in the lives of massive stars that have shed their hydrogen (Type IIb/Ib)  and possibly helium (Type Ic) layers over their evolutionary lifetimes.
The majority of mass loss for single massive stars is
 due to line-driven winds, while  massive stars   in binary systems the majority of mass loss is  likely due to  mass transfer to  the companion via Roche lobe overflow  and/or through common envelope evolution.
Recent studies of samples of SE SNe \citep{drout11,cano13,taddia15a,lyman16,prentice16} indicate
 Type~IIb/Ib/Ic exhibit  a range of  $^{56}$Ni mass ($M_{Ni}$) and
 explosion energy ($E_{K}$), while ejecta mass ($M_{ej}$) estimates are found to span a  more narrow range.
 For example, \citet{lyman16} found from an extended sample of SE SNe that key explosion parameters range as:
 $<M_{\rm Ni}>~=~$0.14--0.29~M$_\odot$, $<M_{ej}>~=~$2.6--3.9~M$_\odot$ and  $<E_{K}>~=~$2.7--5.2$\times 10^{51}$~erg. Only the more extreme  SNe~Ic BL (broad-line) exhibit significantly higher values of $E_K$ and $M_{Ni}$.

The  low ejecta masses typically inferred for the current samples of SE~SNe, along with the direct progenitor detection for a handful of events (e.g., SN~1993J, SN 2008ax, SN~2011dh, SN~2013df, iPTF13bvn, see \citealp{smartt15} and references therein) and the rates of the various sub-types \citep{smith11rates}, suggest their progenitors are linked to  binary star systems \citep[see][for a review]{yoon15}.
However, single massive stars may account for the brightest and most energetic SE~SNe, including those associated with long-duration gamma-ray bursts and some superluminous supernovae. To date no progenitor detections exist of  any SN~Ic.

An important requirement to accurately estimate the luminosity and explosion parameters of any type of SN is an estimate of its reddening due to dust along the line-of-sight to the observer.
Assuming the extinction of any given SN can be defined by a dust law that has a wavelength dependence, the level of reddening experienced by two passbands will  differ.
 Reddening  is typically parameterized by the  $E(B-V)$ color excess and the  absorption coefficient, $R_V$.
 The absorption coefficient dictates the total-to-selective absorption and is  defined as $R_V = A_V / E(B-V)$.
 In this parameterization $E(B-V)$ sets the optical depth of the intervening dust and gas  and $R_V$ is dependent on the average dust grain size  \citep{weingartner01}.
In practice large values of $R_V$ lead to shallow extinction curves and are associated with larger than normal size dust grains, while low values of $R_V$ lead to steep extinction curves due to smaller than normal size dust grains.

Reddening is typically split between two components. The first component accounts for reddening due to dust within the Milky Way (defined by $R_V^{MW}$ and $E(B-V)_{MW}$) and the second component accounts  for reddening external to the Milky Way (defined by $R_V^{host}$ and  $E(B-V)_{host}$).
Reddening along any particular line-of-sight in the Milky Way can be approximately corrected for using the reddening maps of \citet{schlafly11} and \citet{schlafly16}.
However,  accounting for reddening  external to our own galaxy  is significantly more challenging.
In principle light traveling from an extragalactic source may be reddened  by dust located in its  immediate vicinity, dust located in the  host galaxy of the source, or potentially by dust encountered within the intergalactic medium.
Typically the combined effect of reddening due to dust associated with these locations is lumped together and labeled as $E(B-V)_{host}$.
In principle, this could lead to potential problems as dust at any of these locations could have  different properties, leading to significantly different values of $R_V$ and hence various levels of reddening.
Indeed even in the MW the dust in the bulge shows  significant variations of $R_V^{MW}$ \citep{nataf16}.

Efforts to devise robust methods to estimate $A_V^{host}~=~R_V^{host} \cdot E(B-V)_{host}$ for SE SNe
have largely  been hampered by a lack of well-observed and homogeneous samples, and as a result, the inability to identify well-defined intrinsic SN colors.
Typically, host extinction values of SE~SN are inferred from empirically-derived relations between the measured equivalent width of \ion{Na}{i}~D  (hereafter $EW_{\ion{Na}{i}~D}$) and $E(B-V)_{host}$, and assuming a typical Galactic reddening law characterised by $R_V = 3.1$.
These relations often rely on calibrations determined from observations within the Milky Way \citep[e.g.,][]{munari97,poznanski12} and their use to determine reddening in the hosts of extragalactic SNe is suspect as dust properties are known to vary among galaxies \citep{poznanski11}.
Indeed, these relations are associated with significant scatter, providing extinction values  with uncertainties of  $\sim$68\% of  the inferred value of  $A_V$ \citep{phillips13}.
 When combined with extinction estimates computed from Balmer decrement
measurements of the immediate environment of a SN \citep[cf.][]{xiao12} and/or  other less direct indicators, such as the position within the respective host galaxy and/or the SN peak colors, one can arrive at an estimate of $A_V^{host}$, but with significant uncertainty.

\citet{drout11} constructed an  intrinsic $V-R$ color curve template from a small sample of SE SNe.
With this template  they were then  able to define an intrinsic $V-R$ color at $+$10d relative to $V$-band maximum with a scatter of only a few tenths of a magnitude.
This intrinsic color allows for the inference of the $E(V-R)_{host}$ color excess for reddened objects, and when adopting a reddening law, the host extinction.
Unfortunately the photometry of the objects in the \citeauthor{drout11}  sample was computed from images obtained with various telescopes, which is known to lead to considerable calibration issues \citep[e.g.][]{stritzinger02}.
They also adopted  host extinction values for some of the literature-based objects that were  estimated from different \ion{Na}{i}~D vs. $E(B-V)_{host}$ relations.
A similar approach  was employed by \citet{taddia15a} to estimate the host-galaxy reddening for a sample of SE SNe, but in this case  a  $g-r$  intrinsic color template was constructed from photometry of a number of objects  from the literature.

Ideally, one should identify intrinsic SE SN colors from photometry obtained on a single, stable and well-understood photometric system, in addition to expanding the analysis to include various color combinations extending over a range of wavelength in order to facilitate a robust  estimate of both $A_V^{host}$ {\em and} $R_V^{host}$.

In this paper we build upon the work of \citeauthor{drout11} and in doing so present methods to  infer  $E(B-V)_{host}$, $A_V^{host}$, and in some instances, a preferred value of $R_V^{host}$ for SE SNe.
To realize these aims we make use of high-quality, multi-color light curves of over thirty SE SNe  observed by the first phase of the {\em Carnegie Supernova Project} (CSP-I; \citealt{hamuy06}).
The full CSP-I sample of SE SNe is presented in a companion paper \citep{stritzinger17}\footnote{Published CSP-I photometry is readily available in electronic format at our Pasadena-based webpage: \href{http://csp.obs.carnegiescience.edu/}{http://csp.obs.carnegiescience.edu/}}, and consists of optical ($ugriBV$) photometry  of 34\ objects, with a subset of 26\ of these  having at least some near-infrared (NIR, i.e. $YJH$) photometry.
The CSP-I sample is ideal for the present purpose as it has been obtained on a well-understood photometric system \citep{stritzinger11,krisciunas17}, and consists of photometry measured from high signal-to-noise images typically taken under excellent observing conditions.
Furthermore, the CSP-I has put forth a considerable effort to define accurate
local sequences of stars, calibrated relative to standard star fields observed over multiple photometric nights.
Combined with host-galaxy template-image subtraction of each science image, our SN magnitudes have been measured differentially to the local sequence, yielding  photometry with typical  (statistical and systematics) uncertainties ranging between 0.01 to 0.03~mag.
 Analysis of the full the CSP-I photometric light curves and spectroscopy are  presented in additional companion papers by \citet{taddia17} and \citet{holmbo17}, respectively.
 From these works we adopt  for each object (i) estimates  of the epoch of maximum  for each filtered light curve, (ii)  $K$-corrections, and (iii) the spectral classification  of each object.

In Sect.~\ref{sec:color} we examine the optical and NIR colors of the CSP-I SE SNe sample.
This is followed by  Sect.~\ref{sec:intrinsic} where we define intrinsic color-curve templates for the different spectroscopic SE SN sub-types, i.e., Type~IIb, Type~Ib, and Type~Ic.
In Sect.~\ref{sec:extinction} the intrinsic color templates are compared to the observed colors of reddened SE SNe to  compute optical and NIR color excesses. These inferred color excesses are then used to derive  host-extinction  $A_V^{host}$ for all reddened SE SNe in our sample and to constrain $R_V^{host}$ for the most reddened objects.
Finally in Sect.~\ref{sec:discussion}  the results obtained from various  techniques and color-excess combinations are compared and  discussed.

\section{Color curves of CSP stripped-envelope supernovae}
\label{sec:color}

Shown in each of the sub-panels of Fig.~\ref{fig:colors} are  various optical and optical/NIR colors plotted vs. time relative to maximum light for the CSP-I SE SN sample.
All colors have been corrected for Milky Way extinction
 using values listed in NED that originate from the \citet{schlafly11} re-calibration of the \citet{schlegel98} extinction maps of the Milky Way, assuming a \citet[][hereafter F99]{fitzpatrick99} reddening law characterized by $R_V^{MW} = 3.1$.
The plotted SN colors also include time-dilation corrections and $K$ corrections (see \citealp{taddia17}).
To facilitate the construction of optical/NIR colors the NIR light curves were first interpolated as described by \citet{taddia17} and \citet{stritzinger17},  and then evaluated on  the epochs that optical observations were obtained.

The overall shape of the color evolution for each color combination is quite similar.
In the epochs preceding maximum brightness the colors of SE SNe  are typically at their bluest values \citep[e.g.,][]{stritzinger02, galyam04}.
As the SN evolve through maximum light and the temperature of their photosphere drops, the colors evolve towards the red almost linearly.
The increase in colors observed in  Fig.~\ref{fig:colors}  during this time can exceed over two magnitudes for nearly all color combinations except for the $V-i$ color curves which exhibit  a color change of  no more than 1.5 mag.
Upon reaching their maximum red color between two to three weeks past maximum, the color curves make an abrupt turn and evolve  back towards the blue, following a linear evolution with a shallow slope out to later times.

Turning to NIR colors, shown in Fig.~\ref{fig:NIRcolors} (from top to bottom) are the $Y-J$, $Y-H$, and $J-H$ color curves for the  CSP-I SE SN sample.
Inspection of  the  $Y-J$ and $Y-H$ colors reveals an evolution similar to the optical and optical/NIR colors shown in Fig.~\ref{fig:colors}, albeit with more complex behaviour.
At early phases the colors are blue, and as the SN evolve over the next fortnight their colors move towards  the red. This is  followed by a turnover to the blue that is significantly more prevalent than in the case of the optical and optical/NIR color curves.
The $J-H$ color curves, on the other hand,  show blue SN colors at early times, which then evolve consistently towards the red over the entire duration of our photometric coverage.

To obtain a continuous representation of the color evolution the optical and optical/NIR color curves of each SN are fit with an analytic function.
The color curves shown in Fig.~\ref{fig:colors} largely mimic the color evolution of normal SNe~Ia.
We therefore opted to fit each color curve with the same analytical function employed  by \citet{burns14} to study the color properties of the CSP-I sample of SNe~Ia.
The \citet{burns14} analytical function is nearly linear with positive slope at early times, subsequently reaches a maximum value, and then makes a transitions back to a linear function with negative slope.
The functional form is as follows:
\begin{eqnarray}
y(t) = &\frac{\left(s_{0}+s_{1}\right)}{2}t+ \frac{\tau}{2}\left(s_{1}-s_{0}\right)\ln\left[\cosh\left(\frac{t-t_{1}}{\tau}\right)\right]
   \nonumber \\
   & +c+f(t,t_{0}),
   \label{eq:gifit}
\end{eqnarray}

\noindent where $s_0$ is the initial slope, $s_1$ is the final slope, $\tau$ is the length scale over which the transition occurs, and $c$ sets the overall normalization of the function.
The location of the color maximum is defined as $t_{max}=t_1+\tau~{\rm tanh}^{-1}((s_0+s_1)/(s_0-s_1))$.
In order to capture the color evolution at the earliest epochs,
 a second-order polynomial, $f(t,t_0)$,  for $t < t_0$ and equal to 0 for $t  > t_0$,  completes Eq.~\ref{eq:gifit}.
 An example fit of Eq.~\ref{eq:gifit} to the $B-V$ color curve of the Type~IIb SN~2006T is shown in Fig.~\ref{fig:colorfit}.
 The function provides an objective representation of the color evolution, and the  various fit parameters are well-determined for all optical and optical/NIR color combinations.

For some objects the photometric coverage is poor and does not cover the pre-peak phase or does not extend  to the epoch(s) when maximum color is reached.
In these cases the polynomial describing the early color evolution is removed from Eq.~\ref{eq:gifit} before the observed colors are fit.
 In some other instances the color coverage only extends over  several  days, and when so,  the colors are fit with a low-order polynomial.

\section{Intrinsic color-curve templates and temperature constraints from hydrodynamical explosion models}
\label{sec:intrinsic}
\subsection{Color-curve templates}

The construction of intrinsic color-curve templates requires the identification of a  sub-set of objects believed to suffer minimal to no reddening.
To identify such objects we implemented selection criteria based on several factors.
The main two criteria in identifying a minimally reddened SE SN is little to no \ion{Na}{i}~D absorption in  the optical spectrum  and relatively blue optical colors.
If an object fulfilled these two criteria we then also sought confirmation by considering the location of the SN in its host galaxy and the host's inclination.
Objects that fulfilled the first two criteria and were  located in the outskirts of their hosts and/or in minimally inclined galaxies are considered to suffer minimal to no reddening.

In principle, $EW_{\ion{Na}{i}~D}$  can  provide an indication of the host-galaxy color excess \citep[e.g.,][]{munari97,poznanski12}, although the relations used are associated  with large scatter \citep{phillips13}.
With this caveat in mind, plotted in Fig.~\ref{fig:NaID_bmv} (left panel) is $EW_{\ion{Na}{i}~D}$ measured from near-maximum light  spectra for each member of our sample \citep{holmbo17} vs. the $B-V$ color at $+$10d relative to $V_{\rm max}$ (hereafter $B-V_{+10d}$).
We report the  $EW_{\ion{Na}{i}~D}$ values in Table~\ref{tab:NaID} and note no time variations in the measured $EW_{\ion{Na}{i}~D}$ values are found for  objects with multiple epochs of spectra available.
For objects with minimal to no $EW_{\ion{Na}{i}~D}$ a 3$\sigma$ upper limit was computed based on the noise of the spectrum in the region of the \ion{Na}{i}~D line. These limits are plotted as arrows in Fig.~\ref{fig:NaID_bmv}.
The two quantities exhibit a linear trend accompanied with  scatter.
 The best fit linear function between the $EW_{\ion{Na}{i}~D}$ and  $B-V_{+10d}$ (including the upper limits) is
 given by
 \begin{equation}
 EW_{\ion{Na}{i}~D} = 1.568(\pm0.119)\cdot(B-V_{+10d}) - 0.709(\pm0.220).
 \end{equation}

The fact that those SE SNe with minimal to no $EW_{\ion{Na}{i}~D}$ are also the bluest objects within our sample is encouraging in the perspective of using SE SN colors to estimate host extinction.
The  small region in the lower left corner of Fig.~\ref{fig:NaID_bmv} (left panel) defined by  two dashed lines ($B-V_{+10d}~<$~1.05~mag,  $EW_{\ion{Na}{i}~D}~<$~0.4~\AA) contains the objects characterized by both blue color and  minimal to no $EW_{\ion{Na}{i}~D}$.
These objects can be used to construct intrinsic color-curve templates for each of the SE SN  sub-types.
Specifically, contained within the region are three objects for each spectroscopic sub-type and they correspond to
the Type~IIb SN~2005Q, SN~2008aq, and SN~2009Z; the
Type~Ib SN~2004ew (this SN is actually not observed at +10d but its color at epochs earlier than +20d are as blue as those of the other two unreddened SNe~Ib), SN~2007Y, and SN~2007kj; and the
Type~Ic SN~2004fe, SN~2005em, and SN~2008hh.
These SNe are also marked in the right-panel of
 Fig.~\ref{fig:NaID_bmv}, where plotted are the de-projected and normalized SN distances from their respective host galaxy centers vs. their respective host galaxy inclinations (cf. \citealp{taddia17}).
 For eight of the nine objects we have the necessary host-galaxy information  to place them in the right panel of  Fig.~\ref{fig:NaID_bmv} (black circles).
 Each of these SNe appear relatively far from the center of their host, and none of their respective hosts are significantly inclined.
 This provides an additional measure of confidence these objects  suffer minimal host reddening.

With a minimally-reddened subsample  identified,  intrinsic color-curve templates for each of three canonical SE SN sub-types was constructed.
An example is provided in Figure~\ref{fig:templates}  where the $B-V$ color curves of each of the  minimally-reddened Type~IIb, Type~Ib and Type~Ic SNe (hereafter SNe~IIb, SNe~Ib, SNe~Ic) are plotted, along with the intrinsic $B-V$ color templates derived from fitting simultaneously the colors of each SE SN sub-type with Eq.~\ref{eq:gifit}.
Also shown is a comparison  between each of the three template color curves.
 The templates are slightly different among the three sub-types, with the SNe~Ic showing redder colors during the 20 days post maximum and bluer colors after the peak of $B-V$. This is most likely due to differences in the spectral features.

Plotted in Fig.~\ref{fig:spectra_unreddened} are visual-wavelength spectra for 8 of the 9 objects selected as representative of the intrinsic color evolution of the  three sub-types.
This includes post maximum spectra of the Type~IIb SN~2005Q, SN~2008aq and SN~2009Z, the Type~Ib SN~2007Y and SN~2007kj, and the Type~Ic SN~2004fe, SN~2005em and SN~2008hh.
The spectra of the various spectral types are  broadly similar though clear differences between line strengths and line ratios are apparent between the different SN sub-types, particular around the location of H$\alpha$ and the \ion{He}{i} features, with the latter feature being absent in  SNe~Ic \citep[see][]{holmbo17}.

We next proceed to construct  intrinsic color-curve templates  extending over various optical and optical/NIR  color combinations for each SE SN sub-type.
Plotted in  Fig.~\ref{fig:all_templates} are the intrinsic color-curve templates for eight color combinations: $u-V$, $B-V$, $g-V$, $V-r$, $V-i$, $V-Y$, $V-J$, and $V-H$, along with their accompanying 1$\sigma$ uncertainties computed by adding in quadrature the photometry errors.
Color-curve templates were constructed by fitting  Eq.~\ref{eq:gifit} to the various color combinations of each of the SE SN subtypes.
In general  between 0d and $+$20d the different SE SN sub-types exhibit similar --though not identical-- colors and temporal morphology, while at earlier and later epochs color differences are more pronounced.
The average dispersion of our $V-X$ color templates between 0d and $+$20d is in the range 0.02--0.15~mag for SNe~Ib, 0.01--0.14~mag for SNe~IIb, 0.02--0.09~mag for SNe~Ic. The intrinsic $V-X$ colors at +10~d since $V_{max}$ are given in Table~\ref{tab:colors_10d}.  \citet{drout11} report a $V-R$ color at +10 d since $V_{max}$ of 0.26$\pm$0.06 mag, which is very close to our SN~IIb value for $V-r$ at that phase,  and as well is consistent within the errors with the $V-r$ colors of the SN~Ib and SN~Ic minimally reddened sub-samples at the same phase.
The intrinsic $V-X$ color templates for each subtype as well as the $B-X$, $g-X$ and $r-X$ templates (see Sect.~\ref{sec:othercolors}).
\footnote{The intrinsic color-curve templates are available in electronic format on our CSP Pasadena-based webpage: \href{http://csp.obs.carnegiescience.edu/}{http://csp.obs.carnegiescience.edu/}}

\subsection{Minimal post maximum color dispersion}

SE~SNe  potentially arise from multiple progenitors ranging from single to binary star systems and do exhibit a  (modest) range of inferred key explosion parameters.
Nevertheless, the low dispersion in colors around $+$10d
as demonstrated by \citet{drout11}, \citet{taddia15a}, and the present CSP-I sample is observed and an explanation is warranted if we are to
 define and employ intrinsic colors to estimate host-galaxy reddening.
For guidance  we turn to  hydrodynamical explosion models.
In \citet{taddia17} a grid of hydrodynamical models are presented \citep[see][for details on the code]{bersten11} that reproduces the bolometric light curves and velocity evolution  of our SE SNe sample.
The explosion models represent  He-rich progenitor  stars that cover a range of mass, yielding SNe  with $M_{Ni}$ values  ranging between 0.034 and 1.8~$M_{\odot}$,  $M_{ej}$ values ranging  between 1.0 and 6.2~$M_{\odot}$, and explosion energies ranging between $0.6-4\times10^{51}$~erg.
Plotted as diamonds in the left panel of Fig.~\ref{fig:explanation_drout}  are the  $+$10d photospheric temperatures for our series of models  vs.  the ratio between each model's $M_{Ni}$ and $M_{ej}$ values.
The parameters clearly correlate and are well fit  by a  power law as indicated with a dashed black line.
 Also shown  in Fig.~\ref{fig:explanation_drout} (left panel) is the  range  of the $M_{Ni}$ to $M_{ej}$ ratio inferred from simple analytical model fits to the bolometric light curves of several dozen SE SNe  \citep{lyman16}.
These regions are defined by the median and standard deviation values of the $M_{Ni}$ to $M_{ej}$ ratio for each sub-type  \citep[][see their Table~4]{lyman16}, and when compared to the hydrodynamical modeling, imply for each sub-type a rather narrow range in the $+$10d photospheric temperature (e.g., for SN~Ib we find $0.027<~M_{Ni}/M_{ej}~<0.076$).

Assuming  SE SNe can be treated as black-body sources, we compute the $+$10d $B-V$ color for each model and plot it versus the corresponding $M_{Ni}$ to $M_{ej}$ ratio in the right panel of Fig.~\ref{fig:explanation_drout}.
The parameters show a trend that is well fit by a power law function.
Over-plotted as vertical lines is the range of the  $M_{Ni}$ to $M_{ej}$ ratio as inferred by \citet{lyman16}, while the empty dotted points represent the corresponding  mean $B-V$ color value for each sub-type.
Taking the standard deviation of the $B-V$ colors for each sub-type in the allowed range of the $M_{Ni}$ to $M_{ej}$ ratio, the uncertainty in the intrinsic $B-V$ color is estimated to be $\sigma_{B-V} =  0.09, 0.09, 0.05$ and 0.09 mag for SNe~IIb, Ib, Ic and Ic-BL, respectively.
In the case of  SNe~Ic-BL, given the dearth of available models, $\sigma_{B-V}$ was taken to be half the difference between the $B-V$ values of the best fit function (black dashed curve) at the two extremes (magenta dashed lines) of the $M_{Ni}$ to $M_{ej}$ ratio.

If we consider the standard deviation of $B-V$ without distinguishing among SN sub-types, we obtain $\sigma_{B-V} =0.19$~mag. This suggests it is prudent to consider different intrinsic colors for each SE SN sub-type.
Note that the average uncertainty for our three $B-V$ color templates between 0 and $+$20d is 0.02--0.03~mag and that the uncertainty quoted for the $+$10d $V-R$ color by \citet{drout11} is 0.06~mag, which is consistent with our estimates of $\sigma_{B-V}$ for each SE SN sub-type.
This implies that the intrinsic color-curve templates can provide $E(B-V)_{host}$ color excess estimates with a systematic accuracy of $\leq$0.1 mag.
Finally, we note that recently \citet{dessart16} has  published an extended set of SE SN models which also exhibit $+$10d $V-R$ colors characterized by minimal dispersion.

Given the clear correlation between the $M_{Ni}$ to $M_{ej}$ ratio and the $B-V$ color at +10d, we also checked to see if  the minimally-reddened objects identified in this section are  also characterized by a high  $M_{Ni}$ to $M_{ej}$ ratio as compared to the reddened  objects.
If this were the case the minimally-reddened objects should appear blue due to a higher
$M_{Ni}$ to $M_{ej}$ ratio and not for being minimally reddened.
However,  in \citet{taddia17}
the minimally reddening objects  are found to cover the full phase space of the M$_{Ni}$to M$_{ej}$ ratio for each sub-type, suggesting the  color templates are representative of  the  intrinsic colors of  SE SNe.

\section{Host-galaxy color excess, host-galaxy extinction and constraints on $R_V$}
\label{sec:extinction}
With  intrinsic color-curve templates in hand  various color excesses can be inferred.
In the following section we describe how this is done and how color excesses are used to   estimate
the host-extinction  $A_V^{host}$, and constrain $R_V^{host}$ for the most heavily reddened objects.

\subsection{Color excess estimates via intrinsic color-curve templates}
\label{sec:colormini}

To compute the color excess for any particular optical or optical/NIR color combination the observed SN color curve (corrected for Galactic reddening)  is compared to its corresponding intrinsic color-curve template.
As already noted, \citet{drout11} used $+$10d as their fiducial epoch to infer the $E(V-R)_{host}$ color excess.
We on the other hand have elected to compute color excesses by taking the difference between the observed  and  intrinsic colors between 0d and $+20$d.
This temporal interval was chosen for several reasons.
Firstly, the majority of our light curves have their best time sampling in  the two to three weeks just after maximum.
Secondly, during this period the different SE SN sub-types follow a similar evolution in most colors, and thirdly, the overall scatter in the various intrinsic colors is minimal (see Fig.~\ref{fig:all_templates}).
To confirm this is an appropriate manner to proceed we experimented with making templates extending   from 0d to $+$10d and from $+$10d to $+$20d.  In short, these two sets of templates give  results consistent with one another and with the results obtained using the 0d to $+$20d templates.

The technique used to infer the various $E(V-X{_\lambda})_{host}$ (where $X{_\lambda}= u$, $B$, $g$, $r$, $i$, $Y$, $J$, $H$)  color excess values is demonstrated in
Fig.~\ref{fig:get_color_excess_06T}.
Here the observed $V-X{_\lambda}$ color curve of SN~2006T (yellow triangles) is fit with Eq.~\ref{eq:gifit} (black line) and this fit is compared  to the corresponding  SN~IIb  intrinsic-color template (green line).
Specifically, the difference between the two curves in the interval between 0 and $+$20d (where the line representing the template is plotted more thick in Fig.~\ref{fig:get_color_excess_06T}) represents the  $E(V-X{_\lambda})_{host}$ color excess.
 $E(V-X{_\lambda})_{host}$ is computed by taking a weighted average of the difference between the intrinsic color template and the observed colors of the reddened SNe. The adopted weights correspond to the uncertainties associated with the measured observed colors and the estimated epoch of maximum light.

Armed with three sets of intrinsic color-curve templates, one for each of the SE SN  sub-types, the various  $E(V-X{_\lambda})_{host}$ color-excess combinations are computed for the entire sample of reddened SE SNe.
The results are  provided in Table~\ref{tab:EVMX} with accompanying uncertainties computed by adding in quadrature (i) the weighted standard deviation of the difference between the observed and template colors, (ii) the average color uncertainty of the template, (iii) the average uncertainty of the observed color, and (iv) the   uncertainty related to the temporal phase.
Below in Sect.~\ref{sec:othercolors} this approach is extended to estimate color excess values for the entire suite of $B-X{_\lambda}$, $g-X{_\lambda}$ and $r-X{_\lambda}$ color combinations.

 \subsection{Optical-NIR color excesses}
 \label{sec:colorexcesses}

Assuming the dust contained within the host-galaxies of our SE SN sample follows a
F99 reddening law,  a universal value of $R_V^{host}$ can be derived for each of the SE SN spectral sub-types by comparing the color excesses between the optical and NIR bands.
 Plotted in Fig.~\ref{fig:evmx_evmx} are the measured $E(B-V)_{host}$  color excess values  compared with the measured $E(V-X_{\lambda})_{host}$ color excess values (for $X_{\lambda} = i, Y, J, H$) of 18 SE SN.
Over-plotted as lines are the best-fit linear functions for each of the SE SN sub-types where
the corresponding  slopes have been converted  to the $R_V^{host}$ values reported in each panel.
To convert the $E(V-i, V-Y, V-J, V-H)_{host}$ / $E(B-V)_{host}$  slopes to $R_V^{host}$, we adopt the method described by
\citet{krisciunas06} and \citet{folatelli10}.
This requires the calculation of the appropriate $a_{\lambda}$ and $b_{\lambda}$ coefficients \citep[see][their Appendix B]{folatelli10}, which is done using a $+$10d spectral template constructed for each SE SN  sub-type \citep[see][]{holmbo17}.
A summary of this work and a list of the resulting coefficients is provided Appendix~\ref{appendixB}.

The best-fit $R_V^{host}$  value for each SE SN sub-type is found to be strongly dependent on the color-excess combination considered and this is particularly the case of the SNe~Ic.
As the majority of the reddened objects in the sample  suffer relatively moderate amounts of reddening only weak constraints can be placed on $R_V^{host}$.
Nevertheless,  the comparison of optical with NIR color excess measurements in Fig.~\ref{fig:evmx_evmx}  suggests  that the  various SE SN sub-types are  not well represented by a universal $R_V^{host}$ value, although this simple analysis suggests  SNe~Ic favor
$R_V^{host}$ values  that are on average larger than those inferred form the  SNe~IIb and SNe~Ib.
These results are based on small number statistic and conclusions should be taken with some caution, however this type of analysis demonstrates the potential to understand the relationship between SE SN and dust reddening.

\subsection{Host-extinction $A_V^{host}$  from multi-color fits and preferred  values of $R_V^{host}$}
\label{sec:reddening}

The comparison of optical and NIR color excess values in Sect.~\ref{sec:colorexcesses} highlights the difficulty to fit  for $R_V^{host}$ separately for  objects suffering little-to-moderate amounts of reddening.
In principle accurate estimates of both $A_V^{host}$ and $R_V^{host}$ are best obtained from optical and NIR color excess measurements of SN  suffering appreciable levels of reddening, i.e.,
$E(B-V)_{host}$ $\gtrsim 0.20$ mag.
{Seven objects meet the criteria and have both optical and NIR photometry. These include: SN~2004gt (Ic),  SN~2005aw (Ic), SN~2006T (IIb), SN~2006ep (Ib),  SN~2007C (Ib), SN~2009bb (Ic-BL), and SN~2009dt (Ic).
The $E(V-X_{\lambda})_{host}$ color excess values of each of these SN are fit with  the F99 reddening law with $R_V^{host}$ set as a free parameter.
The best non-linear least-squares  fit for each object is plotted in Fig.~\ref{fig:EVMX_vs_lambda_onlyRV} and the corresponding $R_V^{host}$ and $A_V^{host}$ values  are listed in each sub-panel and in Table~\ref{tab:veryred}.
 In general, the shape of the color excess  vs. wavelength relations for the most highly reddened SE SN exhibit similar morphology, and this holds irrespective of spectroscopic subtype.
Figure~\ref{fig:EVMX_vs_lambda_onlyRV} is encouraging and demonstrates the importance of NIR photometry  to estimate the reddening parameter. As of today we are limited by small number statistics and an expanded sample is required to confidently determine the intrinsic properties of SE SNe at NIR wavelengths.

Among the  objects with significant reddening  a range of  $R_V^{host}$ values are obtained.
This includes  an $R_V^{host} = 2.6^{+1.2}_{-0.5}$ for the Type~Ib SN~2007C,  an $R_V^{host} = 3.3^{+0.4}_{-0.3}$ for the Type~Ic BL SN~2009bb and an $R_V^{host} = 4.3^{+0.4}_{-0.3}$
for the  normal Type~Ic SN~2009dt.
Interestingly, this  is  consistent with the findings in Sect.~\ref{sec:colorexcesses} suggesting SNe~Ib and SNe~IIb tend  to occur in  low  $R_V^{host}$ hosts compared to SNe~Ic.
The only SN~IIb with sufficient reddening to determine $R_V^{host}$
is SN~2006T, and as indicated in Fig.~\ref{fig:EVMX_vs_lambda_onlyRV}, it is found to exhibit the  low value of $R_V^{host}=1.1^{+0.2}_{-0.2}$.

The $E(B-V)_{host}$ color excess values presented in
Sect.~\ref{sec:colorexcesses} reveal that the majority of the objects in our sample  suffer relatively low amounts of reddening (i.e., $E(B-V)_{host} < 0.20$ mag).
For these events it is not possible to obtain reliable estimates of $R_V^{host}$ from our fitting method as the associated error  is exceedingly large.
This also holds true for  objects  lacking NIR photometry.
For these objects we therefore assume  their
 host-galaxy dust properties are similar to
those characterizing the dust of the hosts of the most reddened objects in each of the SE SN sub-type samples.
In doing so we adopt $R_V^{host} = 1.1$ for all SNe~IIb (as derived from SN~2006T), $R_V^{host} = 2.6$ for all SNe~Ib (as derived from SN~2007C), and $R_V^{host} = 4.3$ for all SNe~Ic (as derived from in SN~2009dt).
With $R_V^{host}$ set for each SE SN sub-type  each of the low reddened  object's  optical and NIR $E(V-X_{\lambda})_{host}$ sequences  are fit with  the F99 reddening law.

Plotted in Fig.~\ref{fig:EVMX_vs_lambda} are the $E(V-X_{\lambda})_{host}$ color excess values for the SE SNe suffering low extinction (where $X_{\lambda} = u, B, g, R, i, Y, J, H$) vs.  the effective wavelength of passband $X_{\lambda}$.\footnote{The method used to compute an  average effective wavelength for each of the CSP-I passbands is provided in the Appendix~\ref{appendixA}.}
 Over-plotted as solid lines in each panel are the best-fit  F99 models for the $R_V^{host}$ values associated to each sub-type with the associated  1$\sigma$ uncertain indicated by the dash lines.
The corresponding best-fit values of $A_V^{host}$ are reported in each sub-panel and listed in Table~\ref{tab:EBMV}.
From close inspection of the various panels one is able to see the  adopted $R_V^{host}$ values generally provide a good fit to the various color excess values plotted as a function of wavelength.

We next compute  the peak absolute magnitude cumulative distribution function (CDF)  for our sample. To do so CDFs are computed with and without the inclusion of host-extinction corrections. Plotted in  Fig.~\ref{fig:reduceddispersion} are the peak luminosity CDFs for each of the nine CSP-I bandpasses, which we fit assuming a kernel distribution to guide the eye (the best fits are shown by curves of the same color).
Inspection the CDFs  of the blue bands (i.e., $u$, $B$, $g$)  reveals   that host-extinction corrected  CDFs (in blue) are  steeper compared to the uncorrected CDFs (in red),  indicating  their distributions are narrower.
Following the discussion of  \citet{faran14}, this is an encouraging results as a reduction in the dispersion of peak luminosities of the sample population after application of reddening corrections suggests the validity of the extinction corrections.

We now examine the distribution of inferred extinction values for the different SE SN subtypes.
Plotted in  Fig.~\ref{fig:AV} (left-hand panel) are the best-fit $A_V^{host}$ cumulative distribution functions (CDFs) based on the fits shown  in
 Fig.~\ref{fig:EVMX_vs_lambda_onlyRV} and Fig.~\ref{fig:EVMX_vs_lambda}.
Comparing the different CDFs suggests  that on average SNe~Ic are  characterized by  larger $A_V^{host}$ values than SNe~IIb and SNe~Ib.
A Kolmogorv-Smirnov test indicates the significance of these
differences is larger than 2$\sigma$ for the comparison between SNe~Ic and SNe~IIb (p-value~$=$~0.0013). This is true also when we assume $R_V=3.1$ for all the objects suffering low extinction, rather than different $R_V$ based on their subclass. This is demonstrated in the right-hand panel of Fig.~\ref{fig:AV}, where SNe~Ic and SNe~IIb are different in their $A_V$ distributions, with p-value~$=$~0.0465.

\subsection{Host-extinction $A_V^{host}$ from Goobar model}
\label{G08}
Up to here we have assumed the color excess values of our SE SN sample are well reproduce by a F99 reddening law.
However, close inspection of the $E(V-u)_{host}$ color excess values and the best fits shown in Fig.~\ref{fig:EVMX_vs_lambda_onlyRV} and Fig.~\ref{fig:EVMX_vs_lambda} reveal that in a number of cases (e.g., SN~2004gt, SN~2004gv, SN~2005aw, SN~2006ir,  SN~2006lc, and SN~2009dp) the $E(V-u)_{host}$ color excess values are not perfectly reproduced by the F99 law.
This inspired us to fit an alternative reddening law to the observed optical and NIR color excess values. In this instance we adopt the \citet[][hereafter G08]{goobar08} reddening law, which corresponds to a power law taking the functional form:
\begin{equation}
\label{equation:GO08}
 E(V-X)_{host} = A_V^{host} \cdot a(1-(\lambda/\lambda_V)^{-p}).
\end{equation}
\noindent Here $A_V^{host}$ is the host-dust extinction, $\lambda_V$ is the central wavelength of the $V$ filter, and $a$ and $p$ correspond to the power-law fit parameters.
The power-law function described by Eq.~\ref{equation:GO08} was originally crafted to  describe the effects of reddening of SN light produced by multiple scatters off a shell of dust located in the circumstellar environment of  SNe~Ia.
To break the degeneracy between $A_V^{host}$ and $a$ we adopt $a = 1$  \cite[see][]{amanullah15}. As a caveat we note the inferred $A_V^{host}$ values from the power-law fits cannot be taken explicitly as the extinction in the $V$ band, and as well in such a model where part of the extinction is due to scattering of light off of CSM one can infer a large range of $R_V^{host}$ values. If such scattering is common among the SE SN population would then suggest SE SN subtypes do not favor a preferred $R_V^{host}$ value.

Plotted in  Fig.~\ref{fig:goobar}  are the best fit power-law functions (blue line) to 17 SE SNe where the G08 reddening law is found to  provide a  good match to the $E(V-X_{\lambda}$) color excess values.
In each of the instances where the $E(V-u)_{host}$ color excess points are poorly fit by a F99, the G08 model provides a superior fits.

We now examine the $p$ values obtained from the best G08 model fits  of the seven  (see Fig.~\ref{fig:EVMX_vs_lambda_onlyRV}) most highly reddening objects.  Of these objects the best-fit $p$ values of SN~2005aw, SN~2006T, SN~2006ep, and SN~2009dt are found to be inconsistent with the expected $p$ values between  $-1.5$ to $-2.5$ \citep{goobar08}, while SN~2004gt, SN~2007C, and SN~2009bb are found to have best-fit $p$ value confidence intervals that overlap with the range between $-1.5$ and $-2.5$ (see Fig.~\ref{fig:goobar}).
This suggests these later three objects might have CSM dust driving multiple scatterings of light. We note, however, that the best F99 model fits to the color excesses of these objects  also fit quite well (see Fig.~\ref{fig:EVMX_vs_lambda_onlyRV}).

\subsection{Reddening constraints from additional color excess combinations: $E(B-X_{\lambda})_{host}$, $E(g-X_{\lambda})_{host}$ and $E(r-X_{\lambda})_{host}$}
\label{sec:othercolors}

 Most of the current SN follow-up programs do not necessarily observe objects in both SDSS and Johnson passbands. Therefore, we want to consider alternatives to the $V-X_{\lambda}$ color combinations to infer $A_V^{host}$ and $R_V^{host}$. This will, for example, allow one to use our intrinsic color templates also for SNe only observed in the SDSS passbands and not in Johnson $B$ and $V$. Specifically,  in the following  the $B-X_{\lambda}$, $g-X_{\lambda}$, and $r-X_{\lambda}$  color combinations are considered.

Following the procedure used to construct  $V-X_{\lambda}$  intrinsic color-curve templates, intrinsic color-curve templates are constructed for each of the  $B-X_{\lambda}$,  $g-X_{\lambda}$, and $r-X_{\lambda}$ color combinations.
These intrinsic color-curve templates are  used to infer color excess values of reddened SE SNe  relying on the methodology described in Sect.~\ref{sec:reddening}.
 The results of this exercise are summarized in Table~\ref{tab:veryred} and Table~\ref{tab:EBMV}.

In the case of the seven most heavily reddened objects, best-fit $R_V^{host}$  and $A_V^{host}$ values are determined for each color combination. The best-fit results for these objects are listed in Table~\ref{tab:veryred}, as well as  the average $R_V^{host}$ and $A_V^{host}$ values obtained from each color combination.
The average values are characterized by relatively small standard deviations ranging between 0.1--1.1 and 0.03--0.15 for $R_V^{host}$ and $A_V^{host}$, respectively.

Reported in Table \ref{tab:EBMV} are the best-fit $A_V^{host}$ values for each color combination of the 14 SE SN suffering relatively low reddening.
 To obtain the best-fit value a universal $R_V^{host}$ value was assumed for each SE SN sub-type as determined from the most reddened objects (see Sect.~\ref{tab:veryred}).
In addition, the average of the best-fit $A_V^{host}$ values obtained from each of the four color combinations is also given.
 Comparing the $A_V^{host}$ values obtained from the four different color combinations reveals that for most of the objects the $A_V^{host}$ standard deviation is relatively small.
 We therefore  conclude   that when applied to different color combinations  the method gives consistent results.

\section{Discussion}
\label{sec:discussion}

The CSP-I  SE SN light-curve sample has enabled us to devise improved methods to quantify host-galaxy dust extinction, and with the addition of our NIR photometry the added ability to constrain $R_V^{host}$.
We now examine the consistency of the results obtained from the methods presented in Section~\ref{sec:colormini} and Section~\ref{sec:reddening}.
Plotted in Fig.~\ref{fig:checkEBMV} are  comparisons between the best-fit $A_V^{host}$ values obtained from fits of the F99 law to each object's set of color excesses (see the last column of Tables~\ref{tab:veryred} and \ref{tab:EBMV}) and the $A_V^{host}$ values inferred from individual
$E(V-X_{\lambda})_{host}$ color excess measurement and  the F99 law (see Table~\ref{tab:EVMX}).
In doing so these comparisons are made for each of the four color excess combinations:  $E(V-X_{\lambda})_{host}$ (top-left panel), $E(B-X_{\lambda})_{host}$ (top-right panel), $E(g-X_{\lambda})_{host}$ (top-left panel), and $E(r-X_{\lambda})_{host}$ (top-left panel).
Included in each panel's sub-plots is the root-mean-square (rms) difference between the two methods.
In general the  $A_V^{host}$ values obtained by the two methods are in good agreement, as indicated by the fact that the data follow the black lines in the plots.

We find that the the optical and NIR color combinations providing the smallest dispersion among the suite of $E(V-X_{\lambda})_{host}$, $E(B-X_{\lambda})_{host}$, $E(g-X_{\lambda})_{host}$, and $E(r-X_{\lambda})_{host}$ color combinations considered in this work are as follows.
 In the case of $E(V-X_{\lambda})_{host}$ and time relative to $V_{max}$ the $V-r$ and $V-H$ colors  exhibit minimal dispersion in $A_V^{host}$, while  for  $E(B-X_{\lambda})_{host}$ and time relative to $B_{max}$ the $B-i$ and $B-H$ colors exhibit minimal dispersion.
 Moreover, for $E(g-X_{\lambda})_{host}$ and time relative to $g_{max}$ the $g-r$ and $g-H$ colors exhibit minimal dispersion in $A_V^{host}$, while for $E(r-X_{\lambda})_{host}$ and time relative to $r_{max}$ the  $r-g/V$ and $r-H$ colors exhibit minimal dispersion.
 Discrepancies between the two methods are largely due to different spectral lines affecting the $A_V^{host}$ measurement from a single $V/B/g/r-X_{\lambda}$ color combination more than the global $V/B/g/r-X_{\lambda}$  fit value.

We now compare the $A_V^{host}$ value derived for each object in the sample to  their $EW_{\ion{Na}{i}~D}$ measurement obtained from visual-wavelength spectra \citep{holmbo17}.
Specifically, $A_V^{host}$ average values given in Table~\ref{tab:veryred} and Table~\ref{tab:EBMV} are plotted in Fig.~\ref{fig:NAID_EBMV} vs. $EW_{\ion{Na}{i}~D}$.
Over-plotted in the figure  is the relation (solid blue line) found  between these  parameters by \citet{poznanski12}, as well as a linear  fit (dashed line) to the data assuming $EW_{\ion{Na}{i}~D} = 0$ when $A_V^{host}=0$.  The latter of these is given by $A_V^{host}[{\rm mag}]~=~0.78(\pm0.15) \cdot EW_{\ion{Na}{i}~D}[{\rm \AA}]$.
As expected  there is significant scatter between the two quantities \citep[cf.][]{poznanski11}.

Returning to the $A_V^{host}$ CDFs shown in Fig.~\ref{fig:AV}, it is interesting to note that the SNe~Ic consistently exhibit higher extinction values compared to  SNe~IIb and SNe~Ib, and this  holds independent of the assume value of $R_V^{host}$.
This is consistent with the findings of \citet{anderson15} whom  used pixel statistics based on $H_{\alpha}$ emission  to study the environments of SNe by tracing on-going ($<$ 16 Myr old) star formation within SNe host galaxies.
\citeauthor{anderson15} found   SNe~Ic show a statistically significant preference  for being located in regions of on-going star formation, while  SNe~Ib are located in  regions  that do not closely follow ongoing star formation.
The implications are the progenitors of  SNe~Ic  have shorter lifetimes and are more massive compared to the progenitors of SNe~Ib.
Being preferentially located in on-going star forming regions could therefore explain  why SNe~Ic experience (on average) more dust extinction compared with SNe~Ib.
 As a caveat to the SNe~Ic, we stress this subtype contains a mixed bag of transient phenomena, and assuming intrinsic colors based on more normal examples may not be appropriate to apply to more extreme cases such as the broad-line SN~2009bb and SN~2009ca. Unfortunately, we are limited by small number statistics though this will hopefully change in the future.

The findings of Sect.~\ref{sec:extinction}  indicate  SE SNe appear to reside in host galaxies characterized by a range of  dust properties.
Although the sample size of heavily reddened objects is small,  hosts of SNe~Ic are found to be characterized by standard $R_V^{host}$ values, while SNe~IIb and SNe~Ib tend to favor  lower $R_V^{host}$ values.

Low $R_V^{host}$ values have been inferred from the study of some core-collapse SNe. For example, \citet{folatelli14} recently inferred reddening parameters of $R_V^{host} = 1.5$ and $E(B-V)_{host}  = 0.41$ mag for the flat-velocity Type~IIb SN~2010as. Additionally in the process of constructing a low red-shift Hubble diagram populated with SNe~IIP,  \citet{olivares10} obtained a global value of $R_V^{host} = 1.4$.

Low $R_V^{host}$ values are also obtained from the study of
 individual and samples of moderate-to-heavily reddened SNe~Ia, suggesting
 $R_V^{host}$ values ranging from 1.1 \citep{tripp98} to 2.2 \citep{kessler09,guy10,mandel11,burns14} to somewhere  between  $R_V^{host} = 1-2$ \citep{folatelli10}.
 Plotted in Fig.~\ref{RV_cdf} is the $R_V$ CDF  obtained for the most reddened SE SNe in our sample compared to the $R_V$ CDF   obtained from an extended sample of CSP-I SNe~Ia  characterized by $E(B-V)_{host} \geq 0.2$ mag \citep[see][]{burns14}.
 In general, compared to reddened SNe~Ia, SE SNe tend to exhibit larger $R_{V}^{host}$ values.

As noted in Sect.~\ref{sec:reddening} a handful of objects exhibit  $E(V-u)_{host}$ color excess values not reproduced by the fits obtained using the  F99 reddening law.
This motivated us to fit each set of color excess values of our sample with the CSM-motivated reddening law of G08
(see Fig.~\ref{fig:goobar}).
This model was originally introduced to explain the abnormally low $R_V^{host}$ values implied from the study of   some heavily reddened SN~Ia.
Although the G08 reddening law provides superior fits to the $u$ colors of some objects it does not necessarily imply this is due to dust in the SN's immediate circumstellar environment.
For example, \citet{burns14} found for several highly reddened SNe~Ia  with abnormal blue colors that they were   spectroscopically  peculiar objects exhibiting high-velocity features.
This implies the anomalous $V-u$ colors could be driven by the prominent \ion{Ca}{ii} H\&K feature rather than related to the circumstellar environment.
Unfortunately, the spectral coverage of our sample does not allow us to investigate the nature of this spectral feature given the wavelength cut-off off at the blue end of the spectrum for the majority of  the objects occurs at $\sim$ 4000 \AA\ \citep[see][]{holmbo17}.

In summary, we have identified a minimally-reddened sample of SE SNe enabling us to define intrinsic color-curve templates for the  main spectroscopically defined SE SN sub-types.
The minimally-reddened sample is used  to infer the color excess of the reddened objects, and in the case of heavily reddened objects, a constraint on $R_V^{host}$.
This work demonstrates the potential of securely estimating the host reddening parameters of SE SNe and highlights the  added benefits of NIR photometry when it comes to  determining the  host-galaxy reddening properties.
In the  future significant numbers of SE SNe will be discovered, and it is our hope that this work will provide a basis for further efforts  to improve upon our ability to understand their origins.

\begin{acknowledgements}

We thank R. Amanullah,  E. Baron, A. Goobar, and P. Mazzali for useful discussions.
M. D. Stritzinger, F. Taddia, E. Hsiao and C. Contreras gratefully acknowledge  support provided by the Danish Agency for Science and Technology and Innovation realized through a Sapere Aude Level 2 grant.
M. D. Stritzinger acknowledges funding by a research grant (13261) from the VILLUM FONDEN and  the Instrument Centre for Danish Astrophysics (IDA).
M. D.~Stritzinger conducted a portion of this research at the Aspen Center for Physics, which is supported by NSF grant PHY-1066293.
F. Taddia and J. Sollerman gratefully acknowledge the support from the Knut and Alice Wallenberg Foundation.
This material is also based upon work supported by NSF under
grants AST--0306969, AST--0607438, AST--0908886, AST--1008343, AST-1613426, AST--1613455, and AST--1613472.
This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.
\end{acknowledgements}

\bibliographystyle{aa}

redpajama_set_name: RedPajamaArXiv