\section{Introduction}
Developments in ultracold atomic gas experiments have boosted the study of many-body interaction effects in bosonic systems \cite{bloch_dalibar_08,1Dbosonrmp}. In particular, the ability to load Bose-Einstein condensates in optical lattices has provided a unique opportunity to control the effective dimensionality and interactions of Bose gases. A setup that is of interest to us here is that of a BEC loaded in a two-dimensional (2D) optical lattice, as a result of which it splits into an array of cigar-shaped Bose gases at very low temperatures \cite{1Dbosonrmp}. For a sufficiently deep 2D optical lattice and low energies, each cigar-shaped Bose gas can behave as an effective one-dimensional (1D) system \cite{greiner_bloch_01,moritz_stoferle_03} that is described by the Lieb-Liniger model \cite{liebI,liebII} (plus an additional term to account for trapping potentials). The effective 1D interaction strength is determined by the three-dimensional (3D) scattering length and the transverse confinement provided by the lattice \cite{olshaniiLL}. For sufficiently strong effective 1D interactions and sufficiently low 1D densities, these gases enter the so-called Tonks-Girardeau regime \cite{TGcont,kinoshita_wenger_05}, in which bosons behave as impenetrable particles (hard-core bosons) and can be mapped onto noninteracting fermions \cite{TGorigin}. Furthermore, if an additional (weaker, but not too weak) lattice is added along the longitudinal direction, the system can be described by the 1D version of the Bose-Hubbard model, which has been of much interest in condensed matter physics \cite{fisher_weichman_89}. Experiments in the presence of that additional lattice have allowed for the observation of the superfluid to Mott insulator transition in one dimension \cite{stoferle_moritz_04} and the lattice version of the Tonks-Girardeau regime \cite{TGlatt}.

The Lieb-Liniger model has attracted quite some attention since its introduction more than 50 years ago \cite{1Dbosonrmp,yangyang}. It is an integrable model whose exact solution, which can be obtained using the Bethe ansatz \cite{liebI}, provides insights into the universal behavior of 1D gapless systems. Its far-from-equilibrium dynamics has been recently scrutinized theoretically \cite{iyer_andrei_12,caux_konik_12,iyer_guan_13,kormos_shashi_13,nardis_wouters_14,zill_wright_15}, after experiments showed that 1D Bose gases taken far from equilibrium relax to states in which observables are not described by traditional ensembles of statistical mechanics \cite{kinoshita_wenger_06} and generalizations of the latter are needed \cite{rigol_dunjko_07,cassidy_clark_11}.

We should stress that despite the fact that the Lieb-Liniger model is Bethe ansatz solvable, the calculation of correlation functions is extremely challenging due to the complexity of the Bethe ansatz eigenfunctions. In the ground state, for example, the evaluation of one-particle correlations (of which the momentum distribution function is the Fourier transform) has only been possible numerically within algebraic Bethe ansatz \cite{caux_calabrese_07}, and through diffusion Monte Carlo simulations \cite{astrakharchikCF,astrakharchik_giorgini_06}. At finite temperatures, they have been calculated using simulations within the stochastic gauge method \cite{deuar_sykes_09}. Remarkably, the asymptotic behavior of the momentum distribution function [$m(k)$], for large values of the momentum $k$ (associated with short range correlations), has been determined analytically [$m(k)\propto k^{-4}$] \cite{LLtail}. In the Tonks-Girardeau regime, the simplification introduced by the mapping to noninteracting fermions has enabled the analytic calculation of one-particle correlations at zero and finite temperatures \cite{lenard_64,lenard_66,vaidya_tracy_79,jimbo_80}. On the other hand, the 1D Bose-Hubbard model {\it is not} integrable, so correlations can only be computed numerically \cite{1Dbosonrmp}. However, its Tonks-Girardeau limit is integrable and can be mapped onto the so-called $XX$ chain \cite{1Dbosonrmp}, which is the isotropic limit of the $XY$ chain introduced by Lieb, Shultz, and Mattis \cite{lieb_shultz_61}, and which can be mapped onto noninteracting fermions \cite{1Dbosonrmp}. This way, it is possible to calculate analytically one-particle correlations of the Tonks-Girardeau gas in the presence of a lattice \cite{mccoy_68,vaidya_tracy_78}.

In experiments with ultracold gases, a confining potential (that is, to a good approximation, harmonic) is needed to contain the gas. In the presence of such a confining potential, as in the Hubbard model, only the limit in which bosons are impenetrable (hence, mappable to noninteracting fermions) remains integrable. In this limit, the ground-state density and momentum distribution functions, as well as one-particle correlations, have been studied in the continuum \cite{girardeau_wright_01,minguzzi_vignolo_02,lapeyre_girardeau_02,papenbrock_03,forrester_03,gangardt_04} and in the lattice \cite{rigolUni,rigolGS}. In those studies, universal power laws with the same exponents as the ones known from homogeneous systems were found. At finite temperatures, a systematic study of one-particle observables has only been reported in the lattice \cite{rigolFT}, while, in the continuum, the momentum distribution function was recently calculated in Ref.~\cite{FTmdf} for trapped systems with five particles. For finite interaction strengths, for which trapped systems are not integrable anymore, quantum Monte Carlo simulations (with a lattice discretization) have been used to compute momentum distribution functions in the continuum and compare them with experimental results at relatively weak interactions, obtaining a good agreement between the two \cite{MDFboson}.

In this work, we are interested in systematically studying the density and momentum distribution functions of the Lieb-Liniger model in a harmonic trap at zero and finite temperature. We will focus on the scaling properties of those quantities and in the high momentum tails of the momentum distribution function. The latter have remained elusive to quantum Monte Carlo simulations so far. The weight of the high momentum tail, which is known as Tan's contact, is a quantity that plays a central role in a set of universal thermodynamic relations known as Tan's relations~\cite{tanI,tanII,tanIII}. The Tan's contact has been measured in experiments with 3D Bose-Einstein condensates \cite{wild_makotyn_12}.

In order to obtain very accurate results for all the quantities above, for arbitrary interaction strengths, here we use path integral quantum Monte Carlo \cite{pimc} for continuous systems with worm updates \cite{worm1,worm2}. We also obtain results in the Tonks-Girardeau limit using the Bose-Fermi mapping in the lattice and working at very low fillings \cite{rigolUni,rigolGS,rigolFT}. For completeness, in all cases, we also discuss the behavior of the quantities of interest in homogeneous systems, which helps motivate the scaling relations for the trapped systems.

The paper is organized as follows. In Sec.~\ref{sec:hamilQMC}, we introduce the Hamiltonian and numerical approaches used. We also discuss some of the checks done to gauge the accuracy of our calculations. In Secs.~\ref{sec:groundstate} and \ref{sec:finitetemperature}, we report a detailed study of density and momentum distribution functions, as well as one-particle correlations, of Lieb-Liniger systems at zero and finite temperatures, respectively. In both sections, we discuss results for homogeneous and harmonically trapped systems. Finally, our conclusions are presented in Sec.~\ref{sec:conclusion}.

\section{Hamiltonian and numerical approaches} \label{sec:hamilQMC}
We consider 1D bosons with repulsive contact interactions in the presence of an external harmonic trap. The Hamiltonian can be written as
\begin{equation}
 \mathcal{H}=\sum_j\left[-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x_j^2}+V(x_j)\right]+g\sum_{j<l}\delta(x_j-x_l),
\end{equation}
where $m$ is the mass of the bosons, $g$ is the strength of the contact interaction, $V(x_j)=m\omega^2 x_j^2/2$, and $\omega$ is the frequency of the harmonic trap. The strength of the effective 1D contact interaction is usually written as $g=-2\hbar^2/ma_\mathrm{1D}$, where $a_\mathrm{1D}=-a_\perp(a_\perp/a_s-C)$ is the 1D scattering length, $a_s$ is the (3D) $s$-wave scattering length, $a_\perp$ is the length of the transverse confinement, and $C=1.0326$~\cite{olshaniiLL}. Another parameter frequently used to describe the interaction strength is $c=mg/\hbar^2$. In the absence of the trap, this Hamiltonian reduces to the Lieb-Liniger model, which essentially has only one tunable parameter $\gamma=c/\rho$, where $\rho$ is the density.

Here, in order to simulate the Hamiltonian above in the continuum, we use a path integral quantum Monte Carlo method with worm updates (from now on referred to as the worm algorithm) \cite{worm1,worm2}. The configuration space in our simulation is built with discrete imaginary time world lines in a continuous position space. The worm algorithm operates both with diagonal and off-diagonal configurations by introducing one additional open world line in the configuration space. This makes it very efficient in simulating both diagonal observables (such as the energy, the density distribution, and the superfluid
fraction) and off-diagonal observables (such as two-point one-particle correlation functions). In our calculations, the action for contact interactions is approximated by the pair-product action, see Appendix~\ref{app:PPU}, and we always work in the grand canonical ensemble.

In Fig.~\ref{fig:liebenergy}, we compare results for the ground state kinetic [$K(\gamma)=\hbar^2 t(\gamma)/{2m}$] and interaction  [$U(\gamma)=\hbar^2 \nu(\gamma)/{2m}$] energies from the worm algorithm simulations (for a finite system with $N_b\approx20$ bosons) with the exact analytic results in the thermodynamic limit~\cite{liebI}. Average energies are obtained using the thermodynamic estimators reported in Appendix~\ref{app:E_est} (see also Ref.~\cite{pimc}). Figure~\ref{fig:liebenergy} shows that, already for systems with $N_b\approx20$, finite-size effects are negligible, and that the worm algorithm provides a very accurate estimation of observables in diagonal configuration space.

\begin{figure}[!t]
 \includegraphics[width=0.75\linewidth]{EGS_N24.pdf}
 \caption{(Color online) Worm algorithm results for the ground-state kinetic and interaction energies (see text) of the Lieb-Liniger model (circles) compared with the exact analytical results (lines) \cite{liebI}. The statistical errors (shown within the circles) can be seen to be much smaller than the circle sizes. The average number of particles in the system is $N_b\approx20$.}
 \label{fig:liebenergy}
\end{figure}

The observables on which we focus in this work, i.e., the density and momentum distribution functions, are calculated from the one-particle density matrix
\begin{equation}
 g(x,y)=\langle\Psi^\dag(x)\Psi(y)\rangle,
\end{equation}
where $\Psi^\dag(x)$ [$\Psi(x)$] is the bosonic field creation (annihilation) operator at position $x$. $g(x,y)$ is an example of an off-diagonal observable that can be efficiently calculated using the worm algorithm.

For improved statistics, in our calculations, we average $g(x,y)$ over all possible translations ($0<x,y<L$) so that it becomes a function of $r=|x-y|$. For $r>0$,
\begin{equation}
 g(r)=\frac1L\left[\int_0^Lg(x,x+r)dx+\int_0^Lg(y+r,y)dy\right],
\end{equation}
while, for $r=0$,
\begin{equation}
 g(0)=\frac1L\int_0^Lg(x,x)dx.
\end{equation}
The Fourier transform of $g(r)$ gives the momentum distribution function
\begin{eqnarray}
 m(k)&\equiv&\frac{1}{2\pi}\int_0^L dx\int_0^L dy \, e^{ik(x-y)}g(x,y)\nonumber\\&=&
 \frac{L}{2\pi}\int_0^L \cos(kr)g(r)dr.
 \label{eq:mdf_cont}
\end{eqnarray}

\begin{figure}[!t]
 \centering
 \includegraphics[width=0.85\linewidth]{LL_JCaux.pdf}
 \caption{(Color online) (a) Ground-state one-particle correlations and (b) momentum distribution functions of bosons for $\gamma=8.33$. Results are reported for systems with the same density but different average number of particles. Dashed lines are the results obtained using the worm algorithm and solid lines the results obtained in numerical calculations via algebraic Bethe ansatz~\cite{caux_calabrese_07}.}
 \label{fig:LL_JCaux}
\end{figure}

In Fig.~\ref{fig:LL_JCaux}, we show results for $g(r)$ [Fig.~\ref{fig:LL_JCaux}(a)] and $m(k)$ [Fig.~\ref{fig:LL_JCaux}(b)] for the ground state of systems with $\gamma=8.33$, the same density, and different average number of particles, and compare them with numerical results obtained via algebraic Bethe ansatz~\cite{caux_calabrese_07}. The results for $g(r)$ in both approaches are essentially indistinguishable in Fig.~\ref{fig:LL_JCaux}(a). The same is true for $m(k)$ in Fig.~\ref{fig:LL_JCaux}(b), except at the highest momenta. In the latter regime, $m(k)$ becomes very small and the statistical errors in the quantum Monte Carlo simulation become of the same order. This is why it is a challenge to study the asymptotic behavior of $m(k)$ at high momenta. Experiments will suffer from exactly the same limitation as $m(k)$ will at some point become of the same order of the experimental noise. The results in Fig.~\ref{fig:LL_JCaux} make apparent the accuracy of the worm algorithm for computing off-diagonal one-particle observables, which are the focus of this work. Once translational invariance is broken, no analytical solution is available. To gauge the accuracy of the worm algorithm in that case we focus on the Tonks-Girardeau limit.

In that limit, the one-particle density matrix (and, hence, the momentum distribution function) in a lattice can be calculated exactly via the Jordan-Wigner transformation and making use of properties of Slater determinants, both at zero \cite{rigolUni,rigolGS} and finite temperature \cite{rigolFT}. As we explain below, using this approach in the low-density limit in the lattice, one can efficiently compute one-particle properties in the continuum \cite{rigolGS}.

The Tonks-Girardeau Hamiltonian in the lattice can be written as
\begin{equation}
 \mathcal{H}_\mathrm{TG}=-t\sum_i(b_i^\dag b^{}_{i+1}+\text{H.c.})+V\sum_i i^2n_i,
 \label{Eq:H_TG}
\end{equation}
where $b_i^\dag$ ($b_i^{}$) is the creation (annihilation) operator of a hard-core boson at site $i$, $n_i=b_i^\dag b_i^{}$ is the site occupation operator, $t$ is the hopping amplitude, and $V$ sets the strength of the harmonic trap. This Hamiltonian is obtained from the Bose-Hubbard Hamiltonian by taking the limit in which the onsite repulsion $U\rightarrow\infty$ \cite{1Dbosonrmp}.

The hard-core boson creation and annihilation operators satisfy standard bosonic commutation relations with the constraint $(b_i^\dag)^2=b_i^2=0$. Mapping the Tonks-Girardeau Hamiltonian \eqref{Eq:H_TG} onto a spin-1/2 chain \cite{rey_satija_06,he_rigol_11}, and then the spin-1/2 chain onto noninteracting fermions in 1D via the Jordan-Wigner transition, one gets
\begin{equation}\label{Eq:H_F}
 \mathcal{H}_F=-t\sum_i(f_i^\dag f^{}_{i+1}+\text{H.c.})+V\sum_i i^2 n^f_i,
\end{equation}
where $f_i^\dag$ ($f_i^{}$) is the creation (annihilation) operator of a spinless fermion at site $i$, and $n^f_i=f_i^\dag f_i^{}$ is the fermionic site occupation operator. The fermionic Hamiltonian \eqref{Eq:H_F} can be straightforwardly diagonalized. One-particle bosonic correlations can then be obtained using properties of Slater determinants as discussed in Refs.~\cite{rigolUni,rigolGS,rigolFT}. We should stress that, within this lattice approach, the ground-state calculations are done in the canonical ensemble \cite{rigolUni,rigolGS}, while finite-temperature ones are done in the grand-canonical one \cite{rigolFT}.

To establish the relation between the parameters in the lattice Hamiltonian and in the continuum, we notice that the single-particle energy spectrum of the spinless fermion model without a trap is $\varepsilon_k=-2t\cos(ka)$, where $a$ is the lattice spacing. In the zero-density limit, the Fermi momentum $k_F\rightarrow 0$. As a result, the energy spectrum becomes quadratic in $k$ as, up to a constant, we can write $\varepsilon_k=ta^2k^2$. Thus, in the zero-density limit, the lattice model reduces to a continuum model with the effective mass given by $m=\hbar^2/2ta^2$. In this limit, according to Eq.~\eqref{Eq:H_TG}, the external trapping potential is
\begin{equation}
 V(x)=V(x/a)^2 .
\end{equation}
The trapping frequency is then given by $\omega=2\sqrt{tV}/\hbar$, and so the harmonic oscillator length $a_\mathrm{HO}=\sqrt{\hbar/m\omega}$ can be written as $a_\mathrm{HO}/a=(t/V)^{1/4}$. Finally, in finite-temperature calculations in the lattice, the lattice temperature $T_t$ is usually defined as $t\,T_t=k_BT$, i.e., it is given in units of $t$. The temperature in units of trapping frequency in the continuum is related to the lattice temperature by the expression $k_BT/\hbar\omega=T_t\sqrt{t/V}/2$.

\begin{figure}[!t]
 \centering
 \includegraphics[width=0.95\linewidth]{ContvsLatt_Nb5.pdf}
 \caption{(Color online) Density profiles (a) and momentum distribution functions (b) of the Tonks-Girardeau gas as calculated using the lattice approach (red solid lines) and in the continuum using the worm algorithm (dashed lines). The average number of particles is 5. Results are presented for different temperatures: $\tilde{T}\equiv k_BT/\hbar\omega$=0.1, 5, 10. Thin dotted lines indicate the $k^{-4}$ asymptotic behavior of the momentum tails. The statistical errors in the worm algorithm are of the order of the fluctuations seen in the results. They are not reported for clarity.}
 \label{fig:ContvsLatt_Nb5}
\end{figure}

In the lattice, the momentum distribution function is defined as the discrete Fourier transform of the one-particle density matrix
\begin{equation}
 m_k=\frac1N\sum_{jl}e^{ika(j-l)}b_i^\dag b^{}_{j},
 \label{eq:mdf_latt}
\end{equation}
where $N$ is the number of lattice sites. By comparing Eq.~\eqref{eq:mdf_cont} and Eq.~\eqref{eq:mdf_latt}, one can see that the momentum distribution function in the continuum and in the lattice are related by the expression $m_k=m(k)\delta k$, where $\delta k=2\pi/L$ (with $L=Na$) is the discretization of $k$ in the lattice model.

With the relations discussed so far at hand, we are ready to compare results for density and momentum distribution functions in the continuum and the lattice, which, as mentioned before, help us gauge the accuracy of the worm algorithm for diagonal and off-diagonal observables in the absence of translational invariance. In order to approach the Tonks-Girardeau limit in the continuum model, we choose a very small value of the 1D scattering length $a_\mathrm{1D}/a_\mathrm{HO}=5.66\times10^{-3}$.

\begin{figure}[!t]
 \centering
 \includegraphics[width=0.9\linewidth]{c_tgtemp.pdf}
 \caption{(Color online) Tan's contact for the trapped systems in Fig.~\ref{fig:ContvsLatt_Nb5} compared with results in Ref.~\cite{FTmdf}. The contact from the worm algorithm simulation is obtained through the interaction energy using the expression $\mathcal{C}a^3_\mathrm{HO}=-a^3_\mathrm{HO}N_b\nu(\gamma)/(\pi a_\mathrm{1D})$ [derived from Eq.~\eqref{eq:TanC}], where $\nu(\gamma)$ is the average interaction energy per particle [in units of $\hbar^2/(2m)$]. We calculated $\mathcal{C}a^3_\mathrm{HO}$ for two small values of $a_\text{1D}$ obtaining results that agree with the ones in Ref.~\cite{FTmdf}. The contact from the lattice calculation is obtained via a linear fit of the high momentum tails, and then extrapolated to the zero density limit via a quadratic fit to obtain the result for zero $a/a_\mathrm{HO}$ (see inset).}
 \label{fig:c_tgtemp}
\end{figure}

The scaled density [$n(x)=g(x,x)$] and momentum distribution functions obtained in the lattice and the continuum are shown in Fig.~\ref{fig:ContvsLatt_Nb5} at various temperatures ($\tilde{T}\equiv k_BT/\hbar\omega$). The agreement between the results of both approaches is almost perfect at the lowest temperatures $\tilde T=0.1$ and 5 (for $\tilde T=0.1$ the system is essentially in the ground state). As the temperature increases, lattice effects become important and lead to visible differences when the lattice results are compared with the results in the continuum [see the discussion after Eq.~\eqref{eq:TanC}]. In the momentum distribution function [Fig.~\ref{fig:ContvsLatt_Nb5}(b)], the worm algorithm at $\tilde T=0.1$ and 5 already captures part of the $k^{-4}$ momentum tails. For higher temperatures, such as $\tilde T=10$ in Fig.~\ref{fig:ContvsLatt_Nb5}(b), the $k^{-4}$ regime is not visible in the worm algorithm calculation. This is understandable because, as the lattice calculation clearly shows, the $k^{-4}$ tail starts at higher values of the scaled momentum [smaller values of the scaled $m(k)$] as the temperature increases.

We have studied the weight $\mathcal{C}$ (Tan's contact) of the $k^{-4}$ momentum tails [$m(k)=\mathcal{C}\,k^{-4}$] using both methods. In the lattice approach, $\mathcal C$ is computed directly via linear fits of the momentum tails. In the worm algorithm calculation, the contact $\mathcal{C}$ is calculated via the interaction energy estimator, using the relation \cite{FTmdf}
\begin{equation}
 \mathcal{C}=\frac{gm^2}{\pi\hbar^4}\langle\mathcal{H}_\mathrm{int}\rangle .
\label{eq:TanC}
\end{equation}
The results are shown in Fig.~\ref{fig:c_tgtemp}. The values of the contact obtained using the worm algorithm for two small values of the scattering length are consistent with the results in Ref.~\cite{FTmdf}. This indicates that such small values of $a_\text{1D}$ allow us to obtain results in the Tonks-Girardeau limit. Using the lattice approach, the contact in the continuum is calculated by taking the low density limit by means of a quadratic fit (see the inset in Fig.~\ref{fig:c_tgtemp}). The results of such a fitting procedure are in agreement with those obtained using the worm algorithm and Ref.~\cite{FTmdf}, as shown in the main panel in Fig.~\ref{fig:c_tgtemp}. Note that for the ground state calculation, lattice effects are negligible at the lowest densities studied.

\section{Ground-State} \label{sec:groundstate}

\subsection{Homogeneous Systems} \label{sec:homT0}

\subsubsection{Tonks-Girardeau limit}\label{sec:homT0TG}

In the Tonks-Girardeau limit, Tan's contact allows for a simple analytic expression that unveils how to scale the momentum distribution function (one-particle correlations) to achieve data collapse for high momenta (small distances). As indicated in Eq.~\eqref{eq:TanC}, $\mathcal{C}$ is linearly related to the interaction energy. In the language of Ref.~\cite{liebI},
\begin{equation}
 \mathcal{C}=\frac{N_b\rho\gamma\nu(\gamma)}{2\pi},
\label{eq:LLTanC}
\end{equation}
where $\nu(\gamma)=\rho^2\gamma de(\gamma)/d\gamma$, and $e(\gamma)$ is a monotonically increasing function of $\gamma$ that saturates at $\pi^2/3$ when $\gamma\rightarrow\infty$ (for its definition, see Ref.~\cite{liebI}). For large values of $\gamma$, $\nu$ has the following asymptotic form: $\nu(\gamma)=4\rho^2e(\gamma)/(\gamma+2)$. Thus, in the limit $\gamma\rightarrow\infty$, Tan's contact reduces to $\mathcal{C}=2\pi L\rho^4/3$. Hence, $\tilde{\mathcal{C}}\equiv\mathcal{C}/L\rho^4=2\pi/3$. This means that if one plots $m(k)/L$ vs $k/\rho$, no matter the number of particles and the system size, the curves must collapse for high values of $k$. From Eq.~\eqref{eq:mdf_cont}, it then follows that the universal behavior of $m(k)$ for high values of $k$ implies that, if one plots $g(r)/\rho$ vs $\rho r$, the curves must collapse for small values of $r$. Results for $m(k)$ and $g(r)$, obtained using the {\it lattice approach} discussed in the previous section, are reported in Fig.~\ref{fig:mdist_TG}

\begin{figure}[!t]
 \centering
 \includegraphics[width=1.0\linewidth]{mdist_TG.pdf}
 \caption{(Color online) Scaled momentum distribution functions $m(k)/L$ in the Tonks-Girardeau limit as a function of the scaled momentum $k/\rho$. (a) Systems with different number of particles and densities, (b) same as (a) but in log-log scale, (c) systems with the same number of particles and different densities, (d) systems with different number of particles and the same density. Insets: (a) scaled average one-particle density matrix, (c) and (d) same as main panels but in log-log scale. Thin dotted lines indicate the asymptotic behavior discussed in the text. Here, and throughout this work, $\rho$ in the lattice is given in units of the inverse lattice spacing $a$.}
 \label{fig:mdist_TG}
\end{figure}

In Fig.~\ref{fig:mdist_TG}(a), we plot $m(k)/L$ vs $k/\rho$ for systems with different number of particles and densities. Figure~\ref{fig:mdist_TG}(b) reports the same results but in a log-log scale. Both panels make apparent that, for intermediate and high momentum, all curves collapse onto a universal result. In the $k\rightarrow 0$ limit, $m(k=0)/L=m_{k=0}\propto\sqrt{N_b}$ because of the existence of quasi-long-range order \cite{1Dbosonrmp} and, as such, only curves with the same number of particles can collapse for low values of $k$. This can be better seen in Fig.~\ref{fig:mdist_TG}(c), where, for the same number of particles, all momentum distribution curves collapse for all momenta. Differences can only be seen at the highest momenta because of finite-density effects in the lattice. On the other hand, as shown in Fig.~\ref{fig:mdist_TG}(d), plotting results for the same density and different number of particles, again produce curves that collapse at intermediate and high momenta but differ for the lowest values of $k$. For systems with larger number of particles, the universal behavior can be seen to start at lower values of $k$. Consequently, irrespectively of the density in thermodynamically large system sizes, one expects the momentum distribution functions to become universal starting at infinitesimal values of $k$.

In the inset in Fig.~\ref{fig:mdist_TG}(a), we plot results for $g(r)/\rho$ vs $\rho r$ for systems with different number of particles and densities. The plots can be seen collapse onto universal results, in which, at long distance, $g(r)\propto1/\sqrt{r}$ \cite{1Dbosonrmp}. The latter behavior can be understood within the bosonization approach \cite{bosonizing}, which predicts that $g(r)\propto1/r^{1/(2K)}$ at long distances, where $K$ is the so called Tomonaga-Luttinger parameter. In the Tonks-Girardeau limit, $K=1$.

\subsubsection{Finite interaction strength}\label{sec:homT0LL}

\begin{figure}[!b]
 \centering
 \includegraphics[width=1.0\linewidth]{mdist_Gama.pdf}
 \caption{(Color online) Scaled momentum distribution functions $m(k)/L$ vs the scaled momentum $k/\rho$ for finite values of $\gamma$, plotted in [(a) and (c)] a linear scale and in [(b) and (d)] a log-log scale. Results are reported for [(a) and (b)] $\gamma=8.33$ and [(c) and (d)] $\gamma=1.0$. Insets: $g(r)/\rho$ vs $\rho r$ for the same systems as in the main panels. In all plots, thin dotted lines indicate the asymptotic behavior discussed in the text. A linear fit of the long distance behavior of the curves in the insets indicate that the Luttinger parameter is $K\approx1.52$ for $\gamma=8.33$ and $K\approx3.63$ for $\gamma=1.0$. Here, and throughout this work, $\rho$ in the continuum is reported in an arbitrary unit.}
 \label{fig:mdist_Gama}
\end{figure}

For finite interaction strengths, particles can penetrate through one another and the systems tend to have lower kinetic energies (less particles in the high momentum tails) than in the Tonks-Girardeau limit. As a result, Tan's contact decreases with decreasing $\gamma$. From Eq.~\eqref{eq:LLTanC}, $\mathcal{C}$ can be written as $\mathcal{C}=L\gamma^2e'(\gamma)\rho^4/(2\pi)$. The scaled contact, $\tilde{\mathcal{C}}\equiv\mathcal{C}/L\rho^4=\gamma^2e'(\gamma)/(2\pi)$, is now {\it only} a function of $\gamma$. From this result it follows that, so long as $\gamma$ is kept constant (no matter the system size and the number of particles in the system), the curves for $m(k)/L$ vs $k/\rho$ will be universal for large values of $k$. This is the same scaling discussed for the Tonks-Girardeau regime, in which $\gamma$ is constant ($\gamma=\infty$). Similarly, the curves for $g(r)/\rho$ vs $\rho r$ must be universal at short distances. Results for $m(k)$ and $g(r)$, obtained using the {\it worm algorithm}, are reported in Fig.~\ref{fig:mdist_Gama}.

In the main panels in Fig.~\ref{fig:mdist_Gama}, we plot the momentum distribution for $\gamma=8.3$ [Fig.~\ref{fig:mdist_Gama}(a) and \ref{fig:mdist_Gama}(b)] and for $\gamma=1.0$ [Fig.~\ref{fig:mdist_Gama}(c) and \ref{fig:mdist_Gama}(d)], in systems with different densities and average number of particles. Consistent with the previous discussion, the momentum distributions collapse for intermediate and high momentum. More importantly for experiments, with decreasing $\gamma$, the onset of the $k^{-4}$ behavior becomes less sharp (see also Fig.~\ref{fig:mdist_TG}, and Fig.~5 in Ref.~\cite{caux_calabrese_07}). This makes the experimental identification of the $k^{-4}$ tails increasingly difficult as one departs from the Tonks-Girardeau limit.

As $k\rightarrow 0$, the collapse obtained for intermediate and high values of $k$ breaks down because of quasi-long-rage order. For finite values of $\gamma$, more particles have low momenta than in the Tonks-Girardeau limit. In particular, within the bosonization approach, the zero momentum occupation is predicted to be $m(k=0)/L\propto N_b^{1-1/(2K)}$, with $K\geq1$ \cite{bosonizing}. The insets in Figs.~\ref{fig:mdist_TG}(a) and \ref{fig:mdist_TG}(c) depict $g(r)/\rho$ vs $\rho r$ for the same systems for which momentum distributions are shown in the main panels. The curves can be seen to collapse at short distances as expected. At long distances, $g(r)\propto1/r^{1/(2K)}$ as predicted by the bosonization approach, with $K\approx1.52$ for $\gamma=8.3$ and $K\approx3.63$ for $\gamma=1.0$.

\subsection{Trapped Systems}

\subsubsection{Tonks-Girardeau limit}

In order to generalize the scaling discussed in Sec.~\ref{sec:homT0TG} to harmonically trapped systems, we need to find the replacement for $\rho$ and $L$ in the presence of a trap. Let us start by discussing how to scale density profiles in a harmonic trap. Since the density profiles of impenetrable bosons are the same as those of the non-interacting spinless fermions to which they can be mapped \cite{1Dbosonrmp}, we can focus on the noninteracting spinless fermions. The density profiles of the latter systems have been studied in detail in the past \cite{kolomeisky_Newman_00,vignolo_minguzzi_00,brack_zyl_01}. For a homogeneous Fermi system, the density is a linear function of the Fermi momentum $k_F$, $\rho=k_F/\pi$. The chemical potential $\mu$ in such systems is the Fermi energy $E_F=\hbar^2k_F^2/(2m)$, so $\rho=\sqrt{2m\mu}/(\hbar\pi)$.

In a harmonic trap, the local chemical potential changes according to the relation
\begin{equation}
 \mu(x)=\mu_0-V(x),
 \label{eq:eff_CP}
\end{equation}
where $\mu_0$ is the chemical potential at the center of the trap. Within the local density approximation (LDA), the density at each position in the trap is solely determined by the local chemical potential, through the relation obtained for homogeneous systems. Hence
\begin{equation}
 n(x)=\frac{1}{\pi}\left(\tilde{\mu}-\frac{x^2}{a^4_\mathrm{HO}}\right)^{\frac12},
 \label{eq:LDA_TG}
\end{equation}
where $\tilde{\mu}=2m\mu_0/\hbar^2$. The integral of Eq.~\eqref{eq:LDA_TG} over the entire space is the number of particles $N_b$. It provides the following relation between $\tilde{\mu}$ and $N_b$, $\tilde{\mu}=2N_b/a^2_\mathrm{HO}$. The density in the center of the trap is then $\tilde{\rho}=\sqrt{2N_b}/\pi a_\mathrm{HO}$, and the position in the trap at which the density vanishes is $\tilde{L}=\sqrt{2N_b}a_\mathrm{HO}$. Using these quantities, the density profile in the trap can be written as
\begin{equation}
 n(x)=\tilde{\rho}\left[1-\left(\frac{x}{\tilde{L}}\right)^2\right]^{\frac12} .
\end{equation}
By defining the scaled density $\tilde{n}(x)=n(x)/\tilde{\rho}$ and the scaled position $\tilde{x}=x/\tilde{L}$, one obtains
\begin{equation}\label{eq:scalHCB2}
 \tilde{n}(\tilde{x})=(1-\tilde{x}^2)^{\frac12} .
\end{equation}

Since the momentum distribution function at high momenta is determined by the average short-distance one-particle correlations, from the LDA and Eq.~\eqref{eq:scalHCB2} one can advance that, using the scaling relations for homogeneous systems with $\rho\rightarrow\tilde{\rho}$ and $L\rightarrow\tilde{L}$, there will be data collapse for the momentum distribution function at high momenta.

\begin{figure}[!t]
 \centering
 \includegraphics[width=1.0\linewidth]{mdistH_TG.pdf}
 \caption{(Color online) [(a) and (c)] Universal scaling of density and [(b) and (d)] momentum distribution functions in the Tonks-Girardeau limit. We report results from lattice calculations for [(a) and (b)] up to 25 bosons, and from the worm algorithm for [(c) and (d)] up to 100 bosons. In (a) and (c), we also report the LDA prediction for the density profiles. In (b) and (d), thin dotted line depict $k^{-4}$ behavior.}
 \label{fig:mdistH_TG}
\end{figure}

In Fig.~\ref{fig:mdistH_TG}, we report results for density and momentum distribution functions in the Tonks-Girardeau limit obtained using the lattice approach and the worm algorithm. Given the constraint of very low site occupancies in the lattice, our lattice calculations can only be done for significantly smaller numbers of particles than within the worm algorithm. The normalized density profiles in Figs.~\ref{fig:mdistH_TG}(a) and \ref{fig:mdistH_TG}(c) exhibit the expected collapse, even for as few as nine bosons [Fig.~\ref{fig:mdistH_TG}(a)]. Data collapse can also be seen in the momentum distribution functions, except at very low momenta. The $k^{-4}$ high momentum tails are apparent both in the lattice and worm algorithm results, though they are significantly better seen in the lattice calculation where much smaller values of $m(k)$, for higher values of $k$, can be computed.

We note that, using the relations we established in Sec.~\ref{sec:hamilQMC} between quantities in the continuum and in the lattice, one finds that the scaling properties of trapped systems in the continuum discussed so far are consistent with those discussed for lattice systems in Ref.~\cite{rigolUni,rigolGS}.

\subsubsection{Finite interaction strength}

For finite interaction strengths, the relation between the chemical potential and the density is not as simple as in the Tonks-Girardeau limit, where $\mu=\hbar^2\pi^2\rho^2/(2m)$. Instead, $\mu=\hbar^2\rho^2f(\gamma)/(2m)$, where $f(\gamma)=3e(\gamma)-\gamma e'(\gamma)$ \cite{liebI}. This means that, using Eq.~\eqref{eq:eff_CP} and within LDA,
\begin{equation}
 n^2(x)f[\gamma(x)]=n^2(0)f[\gamma(0)]-\frac{x^2}{a^4_\mathrm{HO}} .
\end{equation}
The position in the trap at which the density vanishes is $L_\text{LL}=\sqrt{f[\gamma(0)]}n(0)a_\mathrm{HO}^2$. Defining, as for the Tonks-Girardeau limit, $\tilde{n}(x)=n(x)/n(0)$ and $\tilde{x}=x/L_\text{LL}$, and noticing that $\gamma(x)=mg/[\hbar^2n(x)]=\gamma(0)/\tilde{n}(x)$, we obtain
\begin{equation}
 \tilde{n}^2(\tilde{x})f\left[\frac{\gamma(0)}{\tilde{n}(\tilde{x})}\right]=f[\gamma(0)](1-\tilde{x}^2) .
 \label{eq:LDA_Gama}
\end{equation}

Equation~\eqref{eq:LDA_Gama} makes apparent that the scaled density profiles of systems with the same value of $\gamma$ in the center of the trap must be identical. As in the Tonks-Girardeau limit, this result and the LDA allow one to advance a scaling collapse in the momentum distribution of trapped systems at high momenta \cite{LLtail}.

However, the scaling above requires the knowledge of $n(0)$ and $L_\text{LL}$, which can only be obtained numerically (or measured experimentally). We can predict another universal way to scale densities and positions by noting that, from Eq.~\eqref{eq:LDA_Gama}, it follows that
\begin{equation}\label{eq:ldallm}
 \tilde{n}(\tilde{x})=h_0[1-\tilde{x}^2,\gamma(0)],
\end{equation}
where $h_0(x,y)$ is an unknown function that can, in principle, be computed numerically (from now on we will call dummy functions of this kind as $h_\alpha$, with $\alpha=0,1,\ldots$). Since $N_b=2\int_0^{L_\text{LL}}n(x)dx=2n^2(0)a^2_\text{HO}\sqrt{f[\gamma(0)]}\int_0^1h_0[1-\tilde{x}^2,\gamma(0)]d\tilde{x}$, we see that $n(0)=h_1[\gamma(0)]\sqrt{N_b}/a_\text{HO}=h_2[\gamma(0)]\tilde{\rho}$, and $\tilde{\rho}$ is the density in the center of the trap in the Tonks-Girardeau limit. From this result for $n(0)$ it follows that $L_\text{LL}=h_3[\gamma(0)]\tilde{L}$, where $\tilde{L}$ is the point at which the density vanishes in the Tonks-Girardeau limit. Hence, after fixing $\gamma(0)$, one can use exactly the same scaling for finite interaction strengths as the one used in the Tonks-Girardeau limit. We should stress that this is unique to bosons in the continuum. In the Bose-Hubbard model, one needs two parameters (the so-called characteristic density and the on-site interaction strength) to be kept fixed in order to be able to scale density profiles \cite{rigol_batrouni_09}.

\begin{figure}[!t]
 \centering
 \includegraphics[width=1.0\linewidth]{mdistH_Gama.pdf}
 \caption{(Color online) [(a) and (c)] Universal behavior of scaled density profiles and [(b) and (d)] momentum distribution functions for harmonically trapped systems with $\gamma(0)\approx9.3$ in (a) and (b), and $\gamma(0)\approx0.5$ in (c) and (d). In (b) and (d), thin dotted lines depict $k^{-4}$ behavior.}
 \label{fig:mdistH_Gama}
\end{figure}

In Fig.~\ref{fig:mdistH_Gama}, we show the scaled density [Figs.~\ref{fig:mdistH_Gama}(a) and \ref{fig:mdistH_Gama}(c)] and momentum [Figs.~\ref{fig:mdistH_Gama}(b) and \ref{fig:mdistH_Gama}(d)] profiles for systems with different numbers of particles, trapping frequencies, and 1D scattering lengths. Results are presented for $\gamma(0)\approx9.3$ [Figs.~\ref{fig:mdistH_Gama}(a) and \ref{fig:mdistH_Gama}(b)] and $\gamma(0)\approx2.9$ [Figs.~\ref{fig:mdistH_Gama}(c) and \ref{fig:mdistH_Gama}(d)]. The plots exhibit an excellent data collapse in the density profiles, and in the momentum distribution functions at intermediate and high momentum, as predicted by the analysis above. The high momentum $k^{-4}$ tails are barely apparent in Figs.~\ref{fig:mdistH_Gama}(b) and \ref{fig:mdistH_Gama}(d). By comparing the results for $\gamma(0)\approx9.3$ to those in the Tonks-Girardeau limit in Fig.~\ref{fig:mdistH_TG}, one can see that, even for that large value of $\gamma(0)$, the range and sharpness of the $k^{-4}$ tails have decreased. This worsens as $\gamma(0)$ is further decreased and, as in the homogeneous case, poses an increasing experimental challenge for detecting the $k^{-4}$ tails as one departs from the Tonks-Girardeau regime.

\section{Finite Temperature}\label{sec:finitetemperature}

\subsection{Homogeneous Systems}

\subsubsection{Tonks-Girardeau limit}

At finite temperature, the long-distance behavior of one-particle correlation functions is qualitatively different from the one in the ground state. This because quasi-long-range order is destroyed by thermal fluctuations, i.e., the correlation functions decay exponentially with increasing distance \cite{1Dbosonrmp}. As a result, the population of all momentum modes $m(k)/L$ becomes intensive. This means that, if the density and the temperature are kept constant, the momentum distribution function of systems with different number of particles should collapse (provided the system sizes are much larger than the correlation length). The curves with the same value of $\rho$ in Figs.~\ref{fig:mdistT_TG}(a) and \ref{fig:mdistT_TG}(b) show that this is indeed the case.

\begin{figure}[!b]
 \centering
 \includegraphics[width=1.0\linewidth]{mdistT_TG.pdf}
 \caption{(Color online) Finite temperature results for the scaled momentum distribution function of the homogeneous Tonks-Girardeau gas in systems with different average number of particles and densities. The scaled temperature in the plots is: (a) $\tilde T=2.51$ and (b) $\tilde T=15.06$. The black solid line in both panels was obtained using the lattice approach, while all other curves were obtained using the worm algorithm. Thin dotted lines depict $k^{-4}$ behavior.}
 \label{fig:mdistT_TG}
\end{figure}

The question we would like to address here is how to scale the momentum profiles when the density and the temperature are changed in these systems. For that, as we did in Sec.~\ref{sec:homT0TG}, we look into Tan's contact in the limit $\gamma\rightarrow\infty$. As follows from the derivation in Appendix~\ref{app:YY_bosons}, the Tan's contact takes the form~\cite{LLreview,guanpolylog,guancontact}
\begin{equation}
 \mathcal{C}=\frac{(2mk_BT)^2L}{\hbar^42\pi}f_{1/2}\left(\frac{\mu}{k_BT}\right)f_{3/2}\left(\frac{\mu}{k_BT}\right) ,
\end{equation}
where $f_{\nu}(\cdot)$ is the Fermi-Dirac function (see Appendix~\ref{app:YY_bosons} for its definition). From Eq.~\eqref{eq:fermi}, one can see that $\mu/(k_BT)$ is a function of $2mk_BT/(\hbar^2\rho^2)$. Thus, one can define the scaled (dimensionless) temperature $\tilde{T}=2mk_BT/(\hbar^2\rho^2)$. The scaled contact $\tilde{\mathcal C}=C/(L\rho^4)$ then has the form
\begin{equation}
 \tilde{\mathcal C}=\frac{\tilde{T}^2}{2\pi}f_{1/2}\left(\frac{\mu}{k_BT}\right)f_{3/2}\left(\frac{\mu}{k_BT}\right) ,
\end{equation}
which is only a function of $\tilde T$. This means that, once $\tilde{T}$ is fixed, $m(k)/L$ vs $k/\rho$ is universal. The results in Figs.~\ref{fig:mdistT_TG}(a) and \ref{fig:mdistT_TG}(b), in which we plot systems with different densities and temperatures but the same values of $\tilde{T}$, show that this is indeed the correct scaling. The $k^{-4}$ momentum tails are clearly visible at high momenta. Two important features of those tails are apparent when comparing Figs.~\ref{fig:mdistT_TG}(a) and \ref{fig:mdistT_TG}(b). The first one is that the weight of the tails (the contact) increases with increasing $\tilde{T}$, as discussed in Ref.~\cite{FTmdf}. The second one, probably more important for experiments, is that with increasing $\tilde{T}$ the values of $m(k)$ at which the tails start to develop become smaller. This means that the $k^{-4}$ momentum tails are more difficult to observe experimentally as the temperature increases, despite the fact that Tan's contact becomes larger.

\subsubsection{Finite interaction strength}

\begin{figure}[!b]
 \centering
 \includegraphics[width=1.0\linewidth]{mdistT_Gama.pdf}
 \caption{(Color online) Finite temperature results for the scaled momentum distribution function of the homogeneous Lieb-Liniger gas with different average number of particles and densities at finite $\gamma=8.33$. The scaled temperature in the plots is: (a) $\tilde{T}=2.51$ and (b) $\tilde T=15.06$. All results were obtained using the worm algorithm. Thin dotted lines depict $k^{-4}$ behavior.}
 \label{fig:mdistT_Gama}
\end{figure}

For systems with a finite value of $\gamma$ and finite temperatures one could advance, based on the results for the ground state (in which scaling collapse was found so long as $\gamma$ was kept fixed) and on the results for finite temperature in the Tonks-Girardeau limit (in which scaling collapse was found so long as $\tilde{T}$ was kept fixed), that scaling collapse requires us to keep $\gamma$ and $\tilde{T}$ constant. The results reported in Fig.~\ref{fig:mdistT_Gama} for a finite value of $\gamma=8.33$, which parallel the results in the Tonks-Girardeau limit reported in Fig.~\ref{fig:mdistT_TG}, show that this is indeed the case. The behavior of the high momentum tails with increasing $\tilde{T}$ in Fig.~\ref{fig:mdistT_Gama} is qualitatively similar to that discussed in the Tonks-Girardeau limit.

While a close expression for Tan's contact for finite values of $\gamma$ and finite temperatures is not available,
a calculation up to order $1/\gamma$ (see Appendix~\ref{app:YY_bosons}) reveals that $\tilde{\mathcal C}=C/(L\rho^4)$ can be written as~\cite{LLreview,guanpolylog,guancontact}
\begin{equation}
 \tilde{\mathcal C}=\frac{\tilde{T}^2}{2\pi}f_{3/2}\times\left[f_{1/2}+\frac{\tilde{T}^{1/2}}{\sqrt\pi\gamma}(2f^2_{1/2}+f_{-1/2}f_{3/2})\right] ,
\end{equation}
where by $f_{\nu}$ it is meant $f_{\nu}(\mu/k_BT)$. This further supports the correctness of the scaling used in Fig.~\ref{fig:mdistT_Gama}.

\subsection{Trapped Systems}

\subsubsection{Tonks-Girardeau limit}

\begin{figure}[!b]
 \centering
 \includegraphics[width=1.0\linewidth]{mdistHT_TG.pdf}
 \caption{(Color online) Finite temperature results for [(a) and (c)] the scaled density profiles and [(b) and (d)] momentum distribution functions of harmonically trapped Tonks-Girardeau gases with [(a) and (b)] $\tilde{T}(0)\approx2.63$ and [(c) and (d)] $\tilde{T}(0)\approx37.48$. In (a) and (c), black solid lines show the LDA predictions. In (b) and (d), the black solid lines show results of the lattice calculations. All other curves were obtained using the worm algorithm. In (b) and (d), thin dotted lines depict $k^{-4}$ behavior.}
 \label{fig:mdistHT_TG}
\end{figure}

As for the ground-state case, let us discuss how to scale density profiles in the Tonks-Girardeau limit at finite temperature. Dividing Eq.~\eqref{eq:eff_CP} by $k_BT$, we obtain
\begin{equation}
 \frac{\mu(x)}{k_BT}=\frac{\mu(0)}{k_BT}-\frac{m\omega^2x^2}{2k_BT},
 \label{eq:eff_CPT}
\end{equation}
Using the LDA and the fact that, as discussed in the previous section, in homogeneous systems $\mu/(k_BT)=h_4(2mk_BT/\hbar^2\rho^2)$, we can rewrite Eq.~\eqref{eq:eff_CPT} as
\begin{equation}
 h_4\left[\frac{\tilde{T}(0)n^2(0)}{n^2(x)}\right]=h_4[\tilde{T}(0)]-\frac{x^2}{a_{HO}^4\tilde{T}(0)n^2(0)},
 \label{eq:eff_CPT1}
\end{equation}
where $\tilde{T}(0)=2mk_BT/[\hbar^2n^2(0)]$. Now, let us find the position $L_\text{FT}$ in the trap at which the chemical potential vanishes. At zero temperature, the position at which the chemical potential vanishes is also the position at which the density vanishes. However, this is not the case anymore at finite temperature. For $\mu(x)=0$, $h_4[2mk_BT/\hbar^2n^2(x)]=0$ and, using Eq.~\eqref{eq:eff_CPT1}, $L_\text{FT}=\sqrt{\tilde{T}(0)h_4[\tilde{T}(0)]}n(0)a_{HO}^2$. Defining $\tilde{x}=x/L_\text{FT}$, Eq.~\eqref{eq:eff_CPT1} can be rewritten as
\begin{equation}
 h_4\left[\frac{\tilde{T}(0)n(0)^2}{n(x)^2}\right]=h_4[\tilde{T}(0)](1-\tilde{x}^2).
 \label{eq:eff_CPT2}
\end{equation}
Defining $\tilde{n}(x)=n(x)/n(0)$, we see that Eq.~\eqref{eq:eff_CPT2} implies that
\begin{equation}
 \tilde{n}(x)=h_5[1-\tilde{x}^2,\tilde{T}(0)].
 \label{eq:eff_CPT3}
\end{equation}
At this point, we note that Eq.~\eqref{eq:eff_CPT3} has exactly the same form as Eq.~\eqref{eq:ldallm}, with $h_0\rightarrow h_5$ and $\gamma(0)\rightarrow\tilde{T}(0)$. Using that $N_b=2\int_0^{\infty}n(x)dx$, we get that $n(0)=h_6[\tilde{T}(0)]\tilde{\rho}$ and $L_\text{FT}=h_7[\tilde{T}(0)]\tilde{L}$. This means that, once $\tilde{T}(0)$ has been fixed, a scaling as the one used for the ground state case will produce universal density profiles at finite temperature. In other words, for a fixed $\tilde{T}(0)$, $\tilde{\rho}$ and $\tilde{L}$ determine the average density in the system and an effective system size. Hence, given the fact that LDA is expected to work for sufficiently large system sizes at any finite temperature (because of finite correlation lengths), plots of $m(k)/\tilde{L}$ vs $k/\tilde{\rho}$ are expected to exhibit data collapse.

In Figs.~\ref{fig:mdistHT_TG}(a) and \ref{fig:mdistHT_TG}(c), we show plots of $n(x)/\tilde{\rho}$ vs $x/\tilde{L}$ for systems with different average number of particles, densities in the center of the trap, and temperatures, but all with the same value of $\tilde{T}(0)\approx2.63$ [Fig.~\ref{fig:mdistHT_TG}(a)] and $\tilde{T}(0)\approx37.48$ [Fig.~\ref{fig:mdistHT_TG}(c)]. All the scaled density profiles exhibit the predicted data collapse and are in excellent agreement with the LDA prediction. In Figs.~\ref{fig:mdistHT_TG}(b) and \ref{fig:mdistHT_TG}(d), we show the corresponding scaled momentum distribution functions, which also exhibit an almost perfect data collapse. The results for the momentum distribution functions exhibit an excellent agreement between the worm algorithm calculations and the lattice approach results. The latter allow us to reach smaller values of $m(k)/\tilde{L}$ (larger values of $k/\tilde{\rho}$) so that the $k^{-4}$ momentum tails can be better resolved.

\subsubsection{Finite interaction strength}

\begin{figure}[!t]
 \centering
 \includegraphics[width=1.0\linewidth]{mdistHT_Gama.pdf}
 \caption{(Color online) Finite temperature results for [(a) and (c)] the scaled density profiles and [(b) and (d)] momentum distribution functions of harmonically trapped Lieb-Liniger gases with [(a) and (b)] $\gamma(0)\approx9.3$, and $\tilde{T}(0)\approx1.56$ and [(c) and (d)] $\tilde{T}(0)\approx21.38$. All curves were obtained using the worm algorithm. In (b) and (d), thin dotted lines depict $k^{-4}$.}
 \label{fig:mdistHT_Gama}
\end{figure}

The scaling of trapped systems with finite interaction strengths at finite temperature is the most challenging one. From the ground state analysis in the presence of a trap we know that $\gamma(0)$ needs to be kept fixed. From the finite temperature analysis of homogeneous systems we learned that $\gamma$ and $\tilde{T}$ need to be fixed, while the finite temperature analysis of the trapped Tonks-Girardeau gas revealed that $\tilde{T}(0)$ needs to be fixed. Putting all this together we can advance that, for finite interaction strength, trapped systems at finite temperature will only exhibit scaling collapse if $\gamma(0)$ and $\tilde{T}(0)$ are fixed. We can also advance, from the analysis of trapped systems with finite interaction strength in the ground state and in the Tonks-Girardeau limit at finite temperature that using $\tilde{\rho}$ and $\tilde{L}$ will allow one to achieve the expected collapse.

In Fig.~\ref{fig:mdistHT_Gama} we show results for the scaled momentum distribution function of trapped systems with the same value of $\gamma(0)$ and two values of $\tilde{T}(0)$. They all exhibit a nearly perfect collapse. The features observed in the momentum tails of the homogeneous case are also apparent here. Namely, Tan's contact increases with increasing temperature but at the same time the $k^{-4}$ tails develop starting at higher values of $k$  and smaller values of $m(k)$.

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

We carried out an unbiased study of one-dimensional trapped bosons in the ground state and at finite temperature using the worm algorithm, and the Bose-Fermi mapping in the presence of a lattice at low fillings. Whenever possible, we compared our results to the predictions of the Bethe ansatz obtaining an excellent agreement.

Our study focused in the behavior of the density and momentum distribution functions, as well as one-particle correlations. We discussed in detail how to scale density and momentum profiles to observe universal behavior, which was demonstrated by our numerical results. For trapped systems in the ground state, we showed that fixing $\gamma$ in the center of the trap is all that is required to obtain a universal scaling of density and momentum profiles. This is to be contrasted with the Bose-Hubbard model, for which one needs to fix the so-called characteristic density and the on-site interaction strength \cite{rigol_batrouni_09}. At finite temperature in systems in the continuum, in addition to the condition for the ground state, one needs to fix the scaled temperature in the center of the trap. Those two parameters fully characterize the densities and momentum distribution functions.

In our study of the momentum distribution function, we payed special attention to the $k^{-4}$ asymptotic behavior at high momenta. As shown by our numerical calculations, those momentum tails become increasingly challenging to resolve as the temperature increases and the interactions become weaker, making the Tonks-Girardeau limit at low temperatures the ideal regime for their experimental observation.

\begin{acknowledgments}
This work was supported by the National Science Foundation Grant No.~PHY13-18303 and by the Office of Naval Research. The computations were performed in the Institute for CyberScience at Penn State and the Center for High-Performance Computing at the University of Southern California. We thank Nikolay Prokof'ev for providing us with the worm algorithm code used in the calculations, Jean-S\'ebastien Caux for providing us with the Bethe ansatz results reported in Fig.~\ref{fig:LL_JCaux}, and Juan Carrasquilla and David Weiss for insightful discussions.
\end{acknowledgments}

\newpage\


\section{Introduction}
An \textit{$(m,n)$-colored mixed graph} $G = (V, A \cup E)$ is a graph $G$ with set of vertices $V$, set of arcs $A$ and set of edges $E$ where each arc is
 colored by one of the $m$ colors $\alpha_1, \alpha_2, ..., \alpha_m$ and each edge is  colored by one of the $n$ colors $\beta_1, \beta_2, ..., \beta_n$. We denote  the number of vertices and the number of edges of the underlying graph of $G$ by $v_{G}$  and $e_{G}$, respectively.
Also, we will consider only those
$(m,n)$-colored mixed graphs for which the underlying undirected graph is simple.
Ne\v{s}et\v{r}il and Raspaud~\cite{raspaud_and_nesetril} generalized the notion of vertex coloring and chromatic number
 for $(m,n)$-colored mixed graphs
 by definining colored homomorphism.

Let $G = (V_1, A_1 \cup E_1)$ and $H = (V_2, A_2 \cup E_2)$
be two $(m,n)$-colored mixed graphs.
A colored homomorphism of $G$ to $H$ is a function $f : V_1 \rightarrow V_2$ satisfying
$$uv \in A_1 \Rightarrow f(u)f(v) \in A_2,$$
$$uv \in E_1 \Rightarrow f(u)f(v) \in E_2,$$
and the color of the arc or  edge linking $f(u)$ and $f(v)$ is the same as the color of the arc or the edge linking $u$ and $v$~\cite{raspaud_and_nesetril}.
  We write $ G \rightarrow  H$ whenever there exists a
 homomorphism of $ G$ to $ H$.

Given an $(m,n)$-colored mixed graph $G$ let $H$ be an $(m,n)$-colored mixed graph with minimum \textit{order} (number of vertices) such that $G \rightarrow H$.
Then the order of $H$ is the \textit{$(m,n)$-colored mixed chromatic number} $\chi_{(m,n)}(G)$ of $G$.
For an undirected simple graph $G$, the maximum $(m,n)$-colored mixed chromatic number taken over all $(m, n)$-colored mixed
graphs having underlying undirected simple graph $G$ is denoted by  $\chi_{(m,n)}(G)$.
Let $\mathcal{F}$ be a family of undirected simple graphs.
Then  $\chi_{(m,n)}(\mathcal{F})$ is the maximum of $\chi_{(m,n)}(G)$
taken over all $G \in \mathcal{F}$.

Note that a $(0,1)$-colored mixed graph $G$  is nothing but an undirected simple graph while
$\chi_{(0,1)}(G)$ is the ordinary chromatic number.
Similarly, the study of $\chi_{(1,0)}(G)$
is  the study of oriented chromatic number which is considered by several researchers in the last two decades (for details please check the recent updated survey~\cite{sopena_updated_survey}).
Alon and Marshall~\cite{Marshall-edgecoloring}  studied the homomorphism of $(0,n)$-colored mixed graphs with a particular focus on $n=2$.

A simple graph $G$ is \textit{$k$-acyclic colorable} if we can color its vertices with $k$ colors such that each color class induces an independent set and any two color class induces a forest.
The \textit{acyclic chromatic number} $\chi_a(G)$ of a simple graph $G$ is the minimum  $k$
such that $G$ is $k$-acyclic colorable.
Ne\v{s}et\v{r}il and Raspaud~\cite{raspaud_and_nesetril} showed that $\chi_{(m,n)}(G) \leq k(2m+n)^{k-1}$ where $G$ is a $k$-acyclic colorable graph.
As planar graphs are $5$-acyclic colorable due to Borodin~\cite{Borodinacyclic}, the same authors implied
 $\chi_{(m,n)}(\mathcal{P}) \leq 5(2m+n)^4$ for the family $\mathcal{P}$ of planar graphs as a corollary.
This result, in particular, implies $\chi_{(1,0)}(\mathcal{P}) \leq 80$ and $\chi_{(0,2)}(\mathcal{P}) \leq 80$ (independently proved
before in~\cite{planar80} and~\cite{Marshall-edgecoloring}, respectively).

 Let $\mathcal{A}_k$ be the family of graphs with acyclic chromatic number at most $k$. Ochem~\cite{Ochem_negativeresults} showed that the upper bound  $\chi_{(1,0)}(\mathcal{A}_k) \leq 80$   is tight.
We generalize it for all $(m,n) \neq (0,1)$ to show that the upper bound $\chi_{(m,n)}(\mathcal{A}_k) \leq k(2m+n)^{k-1}$
 obtained by Ne\v{s}et\v{r}il and Raspaud~\cite{raspaud_and_nesetril}
is tight.
This implies that the upper bound $\chi_{(m,n)}(\mathcal{P}) \leq 5(2m+n)^4$ cannot be improved using
the upper bound of $\chi_{(m,n)}(\mathcal{A}_5)$.

The arboricity $arb(G)$ of a graph $G$ is the minimum $k$ such that the edges of $G$ can be decomposed into $k$ forests.
Kostochka, Sopena and Zhu~\cite{Kostochka97acyclicand} showed that given a simple graph $G$, the acyclic chromatic number $\chi_a(G)$ of $G$ is also bounded by a function of  $\chi_{(1,0)}(G)$.
We generalize this result for all $(m,n) \neq (0,1)$ by showing that for a graph $G$ with
$\chi_{(m,n)}(G) \leq k$ we have
$\chi_a(G) \leq k^2 + k^{2 + \lceil log_2 log_p k \rceil}$ where $p = 2m+n$. Our bound slightly improves the bound obtained by  Kostochka, Sopena and Zhu~\cite{Kostochka97acyclicand} for $(m,n) = (1,0)$.
For achieving this result we first establish some relations among arboricity of a graph, $(m,n)$-colored mixed chromatic number and
acyclic chromatic number.

 Let $\mathcal{G}_{\Delta}$ be the family of graphs with maximum degree  $\Delta$.
Kostochka, Sopena and Zhu~\cite{Kostochka97acyclicand} proved that $2^{\Delta/2}  \chi_{(1,0)}(\mathcal{G}_{\Delta}) \leq 2\Delta^2 2^{\Delta}$.
We improve this result in a generalized setting by proving $p^{\Delta/2}   \leq  \chi_{(m,n)}(\mathcal{G}_{\Delta}) \leq 2(\Delta-1)^p p^{\Delta-1}$ for all  $(m,n) \neq (0,1)$ where
$p = 2m+n$.

\section{Preliminaries}

A \textit{special 2-path} $uvw$ of an $(m,n)$-colored mixed graph $G$ is a 2-path satisfying one of the following conditions:

\begin{itemize}
\item[(i)] $uv$ and $vw$ are edges of different colors,

\item[(ii)] $uv$ and $vw$ are arcs (possibly of the same color),

\item[(iii)] $uv$ and $wv$ are arcs of different colors,

\item[(iv)] $vu$ and $vw$ are arcs of different colors,

\item[(v)] exactly one of $uv$ and $vw$ is an edge and the other is an arc.
\end{itemize}

\begin{observation}\label{special}
The endpoints of a special 2-path must have different image under any homomorphism of $G$.
\end{observation}

\begin{proof}
Let $uvw$ be a special 2-path in an $(m,n)$-colored mixed graph $G$. Let $f: G \rightarrow H$ be a colored homomorphism of $G$ to an
$(m,n)$-colored mixed graph $H$ such that $f(u) = f(w)$. Then $f(u)f(v)$ and $f(w)f(v)$  will induce parallel edges in the underlying graph of $H$.
But as we are dealing with $(m,n)$-colored mixed graphs with underlying simple graphs, this is not possible.
\end{proof}

Let $G = (V, A \cup E)$ be an $(m,n)$-colored mixed graph. Let $uv$ be an arc of $G$ with color $\alpha_i$ for some $i \in \{1,2,...,m\}$.
Then $u$ is a
\textit{$-\alpha_i$-neighbor} of $v$ and $v$ is a \textit{$+\alpha_i$-neighbor} of $u$. The set of all $+\alpha_i$-neighbors
and $-\alpha_i$-neighbors of $v$ is denoted by $N^{+\alpha_i}(v)$ and $N^{-\alpha_i}(v)$, respectively. Similarly,
let $uv$ be an edge of $G$ with color $\beta_i$ for some $i \in \{1,2,...,n\}$. Then $u$ is a
\textit{$\beta_i$-neighbor} of $v$ and the set of all $\beta_i$-neighbors of $v$ is denoted by $N^{\beta_i}(v)$.
Let $\vec{a} = (a_1, a_2, ..., a_j)$ be a \textit{$j$-vector} such that
$a_i \in \{\pm \alpha_1, \pm \alpha_2, ...,  \pm \alpha_m,  \beta_1, \beta_2, ..., \beta_n\}$ where $i \in \{1,2,...,j\}$.
Let $J = (v_1, v_2, ..., v_j)$ be a \textit{$j$-tuple} (without repetition) of  vertices  from $G$. Then we define the set
$N^{\vec{a}}(J) = \{v \in V | v \in N^{a_i}(v_i) \text{ \textit{for all} } 1 \leq i \leq j\} $.
Finally, we say that $G$ has property $Q^{t,j}_{g(j)}$ if for each $j$-vector $\vec{a}$ and each $j$-tuple $J$ we have
$|N^{\vec{a}}(J)| \geq g(j)$ where $j \in \{0,1,...,t\}$ and $g: \{0,1,...,t\} \rightarrow \{0,1,... \infty\}$ is an integral function.

  \section{On graphs with bounded acyclic chromatic number}\label{sec acyclic}
  First we will construct examples of $(m,n)$-colored mixed graphs $H_k^{(m,n)}$ with acyclic chromatic number at most $k$ and
  $\chi_{(m,n)}(H_k^{(m,n)}) = k(2m+n)^{k-1}$ for all $k \geq 3$ and for all $(m,n) \neq (0,1)$. This, along with the upper bound established by Ne\v{s}et\v{r}il and Raspaud~\cite{raspaud_and_nesetril}, will imply the following result:

  \begin{theorem}\label{acyclic-chromatic}
  Let $\mathcal{A}_k$ be the family of graphs with acyclic chromatic number at most $k$.
  Then $\chi_{(m,n)}(\mathcal{A}_k) = k(2m+n)^{k-1}$ for all $k \geq 3$ and for all $(m,n) \neq (0,1)$.
  \end{theorem}

\begin{proof}
First we will construct an $(m,n)$-colored mixed graph $H_k^{(m,n)}$, where $p = 2m+n \geq 2$, as follows.
Let $A_{k-1}$ be the set of all $(k-1)$-vectors.
Thus, $|A_{k-1}| = p^{k-1}$.

Define $B_i$ as a set of $(k-1)$ vertices $B_i = \{b^i_1, b^i_2, ...,b^i_{k-1}\}$ for all $i \in \{1, 2, ..., k\}$ such that $B_r \cap B_s = \emptyset$ when $r \neq s$. The vertices of $B_i$'s are called \textit{bottom} vertices for each $i \in \{1, 2, ..., k\}$.
Furthermore, let $TB_i = (b^i_1, b^i_2, ...,b^i_{k-1})$ be a $(k-1)$-tuple.

After that define the set of  vertices $T_i = \{t^i_{\vec{a}}| t^i_{\vec{a}} \in N^{\vec{a}}(TB_i) \text{ \textit{for all} } \vec{a} \in A_{k-1}\}$
for all $i \in \{1, 2, ..., k\}$.
The vertices of $T_i$'s are called \textit{top} vertices for each $i \in \{1, 2, ..., k\}$.
Observe that there are $p^{k-1}$ vertices in $T_i$ for each $i  \in \{1, 2, ..., k\}$.

Note that the definition of $T_i$ already implies some colored arcs and edges between the set of vertices
$B_i$ and $T_i$ for all   $i \in \{1, 2, ..., k\}$.

As $p \geq 2$ it is possible to construct a special 2-path.
Now for each  pair of vertices $u \in T_i$ and $v \in T_j$ ($i \neq j$), construct a special 2-path $uw_{uv}v$ and  call these new vertices $w_{uv}$ as
\textit{internal} vertices for all $i, j \in \{1, 2, ..., k\}$.
This so obtained graph is $H_k^{(m,n)}$.

Now we will show that $\chi_{(m,n)}(H_k^{(m,n)}) \geq k(2m+n)^{k-1}$.
Let $\vec{a} \neq \vec{a'}$ be two distinct
$(k-1)$-vectors.
Assume that the $j^th$ co-ordinate of $\vec{a}$ and  $\vec{a'}$ is different.
Then note that $t^i_{\vec{a}}b^i_jt^i_{\vec{a'}}$ is a special 2-path.
Therefore, $t^i_{\vec{a}}$ and $t^i_{\vec{a}}$ must have different homomorphic image under any homomorphism.
Thus, all the vertices in $T_i$ must have distinct homomorphic image under any homomorphism.
Moreover, as a vertex of $T_i$ is connected by a special 2-path with a vertex of $T_j$ for all $i \neq j$, all the top vertices must have distinct
homomorphic image under any homomorphism. It is easy to see that $|T_i| =  p^{k-1}$ for all   $i \in \{1, 2, ..., k\}$.
Hence $\chi_{(m,n)}(H_k^{(m,n)}) \geq \sum_{i=1}^k |T_i| = k(2m+n)^{k-1}$.

 Then we will show that $\chi_{a}(H_k^{(m,n)}) \leq k$. From now on, by $H_k^{(m,n)}$, we mean the underlying undirected simple graph of the
 $(m,n)$-colored mixed graph $H_k^{(m,n)}$. We will provide an acyclic coloring of this graph with $\{1,2,...,k\}$.
 Color all the vertices of $T_i$ with $i$ for all $i \in \{1,2,...,k\}$. Then color all the vertices of $B_i$ with distinct $(k-1)$ colors
 from the set $\{1,2,...,k\} \setminus \{i\}$ of colors for all $i \in \{1,2,...,k\}$.
Note that each internal vertex have exactly two neighbors. Color each internal vertex with a color different from its neighbors.
It is easy to check that this is an acyclic coloring.

Therefore, we showed that $\chi_{(m,n)}(\mathcal{A}_k) \geq k(2m+n)^{k-1}$
while, on the other hand, Ne\v{s}et\v{r}il and Raspaud~\cite{raspaud_and_nesetril} showed that $\chi_{(m,n)}(\mathcal{A}_k) \leq k(2m+n)^{k-1}$ for all $k \geq 3$ and for all $(m,n) \neq (0,1)$.
  \end{proof}

Consider a complete graph  $K_t$. Replace all its edges  by a 2-path to obtain the graph $S$.
For all $(m,n) \neq (0,1)$, it is possible to assign colored edges/arcs to
the edges of $S$ such that it becomes an $(m,n)$-colored mixed graph with $t$ vertices that are pairwise connected by
a special 2-path. Therefore, by Observation~\ref{special} we know that $\chi_{(m,n)}(S) \geq t$ whereas, it is easy to note that $S$ has arboricity 2.
Thus, the $(m,n)$-colored mixed chromatic number  is not bounded by any function of arboricity.
Though the reverse type of bound exists. Kostochka, Sopena and Zhu~\cite{Kostochka97acyclicand}  proved such a bound for $(m,n) =(1,0)$. We generalize their result for all $(m,n) \neq (0,1)$.

\begin{theorem}\label{chromatic-arboricity}
Let $G$ be an $(m,n)$-colored mixed graph with $\chi_{(m,n)}(G) = k$ where $p = 2m+n \geq 2$. Then $arb(G) \leq \lceil log_{p}k +k/2 \rceil $.
\end{theorem}

\begin{proof}
Let $G'$ be an arbitrary labeled subgraph of $G$ consisting $v_{G'}$ vertices and $e_{G'}$ edges. We know from
Nash-Williams' Theorem~\cite{nash1page}  that the arboricity $arb(G)$ of any graph $G$ is equal to the maximum of
$\lceil e_{G'}/(v_{G'} - 1) \rceil$ over all subgraphs $G'$ of $G$. So it is sufficient to prove that for any subgraph $G'$ of $G$,
$e_{G'}/(v_{G'} - 1) \leq log_p k + k/2$.
As $G'$ is a labeled graph, so there are $p^{e_{G'}}$ different $(m,n)$-colored mixed graphs with underlying graph $G'$. As $\chi_{(m,n)}(G) = k$, there exits a homomorphism from $G'$ to a $(m, n)$-colored mixed graph $G_k$ which has the complete graph on $k$ vertices as its underlying graph.
Note that the number of possible homomorphisms of $G'$ to $G_k$ is at most $k^{v_{G'}}$.
For each such homomorphism of  $G'$ to $G_k$ there are at most $p^{k \choose 2}$ different $(m,n)$-colored mixed graphs with underlying labeled graph  $G'$ as there are $p^{k \choose 2}$ choices of $G_k$.
Therefore,
\begin{equation}\label{eqn wolog}
p^{k \choose 2}. k^{v_{G'}} \geq p^{e_{G'}}
\end{equation}
  which implies
  \begin{equation}\label{eqn wlog}
   log_p k \geq (e_{G'}/v_{G'}) - {k \choose 2} / v_{G'}.
\end{equation}

If $v_{G'} \leq k$, then $e_{G'}/(v_{G'} - 1) \leq v_{G'}/2 \leq k/2$. Now let $v_{G'} > k$. We know that $\chi_{(m,n)}(G') \leq \chi_{(m,n)}(G) = k$. So

\begin{equation*} \label{.}
\begin{split}
log_p k & \geq \frac{e_{G'}}{v_{G'}} - \frac{k(k - 1)}{2 v_{G'}} \\
 & \geq \frac{e_{G'}}{(v_{G'} -1)} - \frac{e_{G'}}{v_{G'}(v_{G'} - 1)} - \frac{k - 1}{2} \\
 & \geq \frac{e_{G'}}{(v_{G'} -1)} - 1/2 - k/2 + 1/2 \\
 & \geq \frac{e_{G'}}{(v_{G'} -1)} - k/2.
\end{split}
\end{equation*}

Therefore, $\frac{e_{G'}}{(v_{G'} -1)} \leq log_{p}k +k/2$.
\end{proof}

We have seen that the $(m,n)$-colored mixed chromatic number of a graph $G$ is bounded by a function of the acyclic chromatic number of $G$. Here we show that
 it is possible to bound the acyclic chromatic number of a graph in terms of its $(m,n)$-colored mixed chromatic number and arboricity.
Our result is a generalization of a similar result proved for $(m,n) =(1,0)$ by  Kostochka, Sopena and Zhu~\cite{Kostochka97acyclicand}.

\begin{theorem}\label{arboricity.chromatic-acyclic}
Let $G$ be an $(m,n)$-colored mixed graph with $arb(G) = r$ and $\chi_{(m,n)}(G) = k$ where $p = 2m+n \geq 2$.  Then
$\chi_a(G) \leq k^{\lceil log_pr \rceil +1}$.
\end{theorem}

\begin{proof}
First we rename the following symbols:
$\alpha_{1} = a_0, -\alpha_{1} = a_{1}, \alpha_2 = a_2, -\alpha_2 = a_3, ..., \alpha_m = a_{2m-2}, -\alpha_m = a_{2m-1}, \beta_1 = a_{2m}, \beta_2 = a_{2m+1},..., \beta_n = a_{2m+n-1}$.

Let $G$ be a graph with $\chi_{(m,n)}(G) = k$ where $2m+n = p$.
Let $v_1, v_2, ..., v_t$ be some ordering of the vertices of $G$.
Now consider the $(m,n)$-colored mixed graph $G_0$ with underlying graph $G$ such that for any $i < j$ we have
$v_j \in N^{a_0}(v_i)$
whenever $v_iv_j$ is
an edge of $G$.

Note that the edges of  $G$ can be covered by $r$ edge disjoint forests $F_1, F_2, ..., F_r$ as $arb(G) = r$.
Let $s_i$ be the number $i$ expressed with base $p$
for all $i \in \{1,2,...,r\}$. Note that $s_i$ can have at most $s = \lceil log_pr \rceil$ digits.

  Now we will construct a sequence of $(m,n)$-colored mixed graphs $G_1, G_2, ..., G_s$ each having underlying graph $G$.
  For a fixed $l \in \{1,2,...,s\}$ we will describe the construction of $G_l$.
  Let $i <j$ and $v_iv_j$ is an edge of $G$.
  Suppose $v_iv_j$ is an edge of the forest $F_{l'}$ for some $l' \in \{1,2,...,r\}$.
  Let the $l^{th}$ digit of $s_{l'}$  be $s_{l'}(l)$. Then  $G_l$ is constructed in a way such that
  we have $v_j \in N^{a_{s_{l'}(l)}}(v_i)$ in $G_l$.

  Note that there is a homomorphism  $f_l: G_l \rightarrow H_l$ for each $l \in \{1,2,...,s\}$ such that $H_l$ is
  an $(m,n)$-colored mixed graph on $k$ vertices.
  Now we claim that $f(v) = (f_0(v), f_1(v), ..., f_s(v))$ for each $v \in V(G)$ is an acyclic coloring of $G$.

  For adjacent vertices $u,v$ in $G$ clearly we have $f(v) \neq f(u)$ as $f_0(v) \neq f_0(u)$.
  Let $C$ be a cycle in $G$. We have to show that at least 3 colors have been used to color this cycle with respect to the coloring given by $f$.
  Note that in $C$ there must be two incident edges $uv$ and $vw$ such that they belong to different forests,
  say, $F_i$ and $F_{i'}$, respectively.
 Now suppose that $C$ received two colors with respect to $f$. Then we must have $f(u) = f(w) \neq f(v)$. In particular we must have
 $f_0(u) = f_0(w) \neq f_0(v)$.
 To have that we must also have $u,w \in N^{a_i}(v)$ for some $i \in \{0,1,...,p-1\}$ in $G_0$.
 Let  $s_i$ and $s_{i'}$ differ in their $j^{th}$ digit. Then in $G_j$ we have $u \in N^{a_i'}(v)$ and $w \in N^{a_i''}(v)$
 for some $i' \neq i''$. Then we must have $f_j(u) \neq f_j(w)$.  Therefore, we also have $f(u) \neq f(w)$. Thus, the cycle $C$ cannot be colored with two colors under the coloring $f$. So $f$ is indeed an acyclic coloring of $G$.
    \end{proof}

Thus, combining   Theorem~\ref{chromatic-arboricity} and~\ref{arboricity.chromatic-acyclic} we have
$\chi_a(G) \leq k^{\lceil log_p \lceil log_{p}k +k/2 \rceil \rceil +1}$ for $\chi_{(m,n)}(G) = k$ where $p = 2m+n \geq 2$.
However,  we managed to obtain the following better bound.

\begin{theorem}\label{chromatic-acyclic}
Let $G$ be an $(m,n)$-colored mixed graph with  $\chi_{(m,n)}(G) = k \geq 4$ where $p = 2m+n \geq 2$.  Then
$\chi_a(G) \leq k^2 +k^{2+ \lceil log_2log_pk \rceil}$.
\end{theorem}

\begin{proof}
Let $t$ be the maximum real number such that there exists a subgraph $G'$ of $G$ with $v_{G'} \geq k^2$ and $e_{G'} \geq t.v_{G'}$.
Let $G''$ be the biggest subgraph of $G$ with $e_{G''} > t.v_{G''}$. Thus, by maximality of $t$, $v_{G''} < k^2$.

Let $G_0 = G - G''$. Hence $\chi_a(G) \leq \chi_a(G_0) + k^2$.
By maximality of $G''$, for each subgraph $H$ of $G_0$, we have $e_{H} \leq t.v_{H}$.

If $t \leq \frac{v_{H} - 1}{2}$, then $e_{H} \leq (t + 1/2)(v_{H} -1)$. If $t > \frac{v_{H} - 1}{2}$, then $\frac{v_{H}}{2} < t + 1/2$. So $e_{H} \leq \frac{(v_{H} -1).v_{H}}{2} \leq (t + 1/2)(v_{H} -1)$. Therefore, $e_{H} \leq (t + 1/2)(v_{H} -1)$ for each subgraph $H$ of $G_0$.

By Nash-Williams' Theorem~\cite{nash1page}, there exists $r = \lceil t + 1/2 \rceil$ forests $F_1, F_2, \cdots, F_r$ which covers all the edges of $G_0$. We know from Theorem~\ref{arboricity.chromatic-acyclic} $\chi_a(G_0) \leq k^{s +1}$ where $s = \lceil log_pr \rceil$.

Using  inequality~(\ref{eqn wlog}) we get $log_p k \geq t - 1/2$. Therfore

$$s = \lceil log_p(\lceil t + 1/2 \rceil)\rceil \leq \lceil log_p(1 + \lceil log_pk \rceil)\rceil \leq 1 + \lceil log_plog_pk \rceil.$$

Hence $\chi_a(G) \leq k^2 +k^{2+ \lceil log_plog_pk \rceil}$.
\end{proof}

Our bound, when restricted to the case of $(m,n) = (1,0)$, slightly improves the existing bound~\cite{Kostochka97acyclicand}.

\section{On graphs with bounded maximum degree}
Recall that  $\mathcal{G}_{\Delta}$ is the family of graphs with maximum degree  $\Delta$.
It is known that $\chi_{(1,0)}(\mathcal{G}_{\Delta}) \leq 2 \Delta^2 2^\Delta$~\cite{Kostochka97acyclicand}.
Here we prove that $\chi_{(m,n)}(\mathcal{G}_{\Delta}) \leq 2 (\Delta-1)^p .p^{(\Delta-1)} + 2 $ for all $p = 2m+n \geq 2$ and $\Delta \geq 5$.
Our result, restricted
to the case $(m,n) = (1,0)$, slightly improves the upper bound of Kostochka, Sopena and Zhu~\cite{Kostochka97acyclicand}.

\begin{theorem}\label{chromatic-degree}
For the family $\mathcal{G}_{\Delta}$  of graphs with maximum degree
$\Delta$ we have  $p^{\Delta/2} \leq \chi_{(m,n)}(\mathcal{G}_{\Delta}) \leq 2 (\Delta-1)^p. p^{(\Delta-1)} +2$ for all $p = 2m+n \geq 2$ and for all $\Delta \geq 5$.
\end{theorem}

If every subgraph of a graph $G$ have at least one vertex with degree at most $d$, then $G$ is \textit{$d$-degenerated}.
Minimum such $d$ is the \textit{degeneracy} of $G$. To prove the above theorem we need the following result.

\begin{theorem}\label{chromatic-degree.degeneracy}
Let $\mathcal{G}'_{\Delta}$  be the family of graphs with maximum degree  $\Delta$ and degeneracy $(\Delta - 1)$. Then
  $\chi_{(m,n)}(\mathcal{G}'_{\Delta}) \leq 2 (\Delta-1)^p .p^{(\Delta-1)}$ for all $p = 2m+n \geq 2$ and for all $\Delta \geq 5$.
\end{theorem}

To prove the above theorem we need the following lemma.

\begin{lemma}\label{key-lemma}
There exists an  $(m,n)$-colored complete mixed graph  with property $Q^{t-1,j}_{1+(t-j)(t-2)}$ on $c = 2 (t-1)^p .p^{(t-1)}$ vertices
where $p = 2m+n \geq 2$ and $t \geq 5$.
\end{lemma}

\begin{proof}
Let $C$ be a random $(m,n)$-colored  mixed graph with underlying complete graph. Let $u,v$ be two vertices of $C$ and the events $u \in N^{a}(v)$
for $a \in \{\pm \alpha_1, \pm \alpha_2, ...,  \pm \alpha_m,  \beta_1, \beta_2, ..., \beta_n\}$ are equiprobable and independent  with probability $\frac{1}{2m+n} = \frac{1}{p}$.
We will show that the probability of $C$ not having property $Q^{t-1,j}_{1+(t-j)(t-2)}$ is strictly less than 1 when
$|C| = c = 2 (t-1)^p .p^{(t-1)}$. Let $P(J,\vec{a})$ denote the probability of the event $|N^{\vec{a}}(J)| < 1+(t-j)(t-2)$ where $J$ is a $j$-tuple of $C$ and $\vec{a}$ is a $j$-vector for some $j \in \{0,1,...,t-1\}$. Call such an event a \textit{bad event}. Thus,

\begin{equation} \label{eq1}
\begin{split}
P(J, \vec{a}) & = \sum\limits_{i=0}^{(t - j)(t - 2)} {c-j \choose i} p^{-ij} (1 - p^{-j})^{c - i - j} \\
 & < (1 - p^{-j})^c \sum\limits_{i=0}^{(t - j)(t - 2)} \frac{c^i}{i!} (1 - p^{-j})^{- i - j} p^{-ij} \\
 & < 2 e^{-cp^{-j}} \sum\limits_{i=0}^{(t - j)(t - 2)} c^i \\
 & < e^{-cp^{-j}} c^{(t - j)(t - 2) + 1}.
\end{split}
\end{equation}

Let $P(B)$ denote the probability of the occurrence of at least one bad event.
To prove this lemma it is enough to  show that $P(B) < 1$. Let $T^j$ denote the set of all $j$-tuples and $W^j$ denote the set of all $j$-vectors.  Then

\begin{equation} \label{eq2}
\begin{split}
P(B) = \sum_{j=0}^{t - 1} \sum_{J \in T^j} \sum_{\vec{a} \in W^j} P(J, \vec{a}) & < \sum\limits_{j=0}^{t - 1} {c \choose j} p^j e^{-cp^{-j}} c^{(t - j)(t - 2) + 1} \\
 & < \sum\limits_{j=0}^{t - 1} \frac{c^j}{j!} p^j e^{-cp^{-j}} c^{(t - j)(t - 2) + 1} \\
 & = 2 \sum\limits_{j=0}^{t - 1} \frac{p^j}{2^j} \frac{2^{j-1}}{j!}  c^j e^{-cp^{-j}} c^{(t - j)(t - 2) + 1} \\
 & < 2 \sum\limits_{j=0}^{t - 1} \frac{p^j}{2^j} e^{-cp^{-j}} c^{(t - j)(t - 2) + 1 + j}.
\end{split}
\end{equation}

Consider the function $f(j) = 2 (p/2)^j e^{-cp^{-j}} c^{(t - j)(t - 2) + 1 + j}$. Observe that $f(j)$ is the $j^{th}$
summand  of the last sum from equation~(\ref{eq2}).
Now
\begin{equation} \label{eq3}
\begin{split}
 \frac{f(j + 1)}{f(j)} & = \frac{p}{2} \frac{e^{(p-1)cp^{-j-1}}}{c^{t-3}} \\
 & >  \frac{p}{2} \frac{e^{(p-1)cp^{-(t-1)}}}{c^{t-3}} \\
 & >  \frac{p}{2} \left(\frac{e^{2 (p-1) (t-1)^{p-1}}}{c}\right)^{t-3} \\
\end{split}
\end{equation}

As $\frac{p - 1}{p} > \frac{1}{2}$,  $$\frac{(k-1)^{p-1}}{2} > \ln (k-1) \implies (p-1)(k-1)^{p-1} > \ln (k-1)^p.$$

Furthermore,  $$\frac{(p-1)}{\ln p} (k-1)^{p-1} > \frac{\ln 2}{\ln p} + (k-1) \implies (p-1) (k-1)^{p-1} > \ln(2 p^{k-1}).$$

 Adding the above two inequalities we get
 $$e^{2 (p-1) (t-1)^{p-1}} > 2 (t-1)^p p^{t-1} = c.$$

Hence $\frac{f(j + 1)}{f(j)} > \frac{p}{2}$. Thus, using inequality~(\ref{eq2}) we get
$P(B)  < \sum\limits_{j=0}^{t - 1} f(j)$.  This implies

\medskip

\[P(B) < \begin{cases}
\frac{(p/2)^t - 1}{(p/2) - 1}f(0), & \text{  if } p > 2 \\
t f(0), & \text{  if } p = 2
\end{cases}
\]

\medskip

\textbf{Case.1:} $p > 2$.

\begin{equation} \label{eq4}
\begin{split}
 P(B) & < 2. \frac{(p/2)^t - 1}{(p/2) - 1}. \frac{c^{(t - 1)^2}}{e^{2(t - 1)^p p^{t - 1}}} \\
 & < 4. \frac{(p/2)^t - 1}{p - 2}. \left(\frac{c}{e^{2 p^{t - 1}}}\right)^{(t - 1)^p} \\
 & <  4. (p/2)^t. \left(\frac{c}{e^{2 p^{t - 1}}}\right)^{(t - 1)^p} \\
 & < \left(\frac{pc}{e^{2 p^{t - 1}}}\right)^{(t - 1)^p}
\end{split}
\end{equation}

Now, we observe that

\begin{equation*} \label{eq5}
\begin{split}
 \ln (pc) & < \ln p + \ln 2 + p\ln (t - 1) + (t - 1)\ln p \\
 & = t \ln p + p\ln (t - 1) + \ln 2 \\
 & <  tp + p(t - 1) + 2 \\
 & < 2tp < 2p^{t - 1}
\end{split}
\end{equation*}

So from the inequality (\ref{eq4}), we can say that $P(B) < 1$ for $p >2$.

\medskip

\textbf{Case.2:} $p = 2$.

\begin{equation} \label{eq6}
\begin{split}
 P(B) & < 2t. \frac{c^{(t - 1)^2}}{e^{(t - 1)^2 2^t}} \\
 & = 2t. \left(\frac{c}{e^{2^t}}\right)^{(t - 1)^2} \\
 & < \left(\frac{2tc}{e^{2^t}}\right)^{(t - 1)^2}
\end{split}
\end{equation}

Observe that, $\ln c = 2 \ln (t-1) + t \ln 2 < 2(t - 1) + 2t = 4t - 2$.

Now, we see that
\[
 \ln (2tc) < 4t - 2 + 2t < 6t < 2^t \\
 \implies 2tc < e^{2^t} \\
 \implies \frac{2tc}{e^{2^t}} < 1
\]

So from the inequality (\ref{eq6}), we can say that $P(B) < 1$ for $p = 2$.
\end{proof}

Now we are ready to prove Theorem~\ref{chromatic-degree.degeneracy}.

\medskip

\noindent \textit{Proof of Theorem~\ref{chromatic-degree.degeneracy}.} Suppose that $G$ is an $(m,n)$-colored mixed graph with maximum degree $\Delta$ and degeneracy
$(\Delta-1)$.
By Lemma~\ref{key-lemma} we know that there exists an $(m,n)$-colored mixed graph $C$ with property
$Q^{\Delta-1,j}_{1+(\Delta-j)(\Delta-2)}$ on $2 (\Delta-1)^p .p^{(\Delta-1)}$ vertices
where $p = 2m+n \geq 2$ and $\Delta \geq 5$. We will show that $G$ admits a homomorphism to $C$.

As $G$ has degeneracy $(\Delta-1)$, we can provide an ordering $v_1, v_2, ..., v_k$ of the vertices of $G$ in such a way that each vertex $v_j$ has at most $(\Delta-1)$ neighbors with lower indices.
Let $G_l$ be the  $(m,n)$-colored mixed graph induced by the vertices $v_1, v_2, ..., v_l$ from $G$ for $l \in \{1,2,...,k\}$.
Now we will recursively construct a  homomorphism $f: G \rightarrow C$ with the following properties:

\begin{itemize}
\item[$(i)$] The partial mapping $f(v_1), f(v_2), ..., f(v_l)$ is a homomorphism of  $G_l$ to $C$ for all $l \in \{1,2,...,k\}$.

\item[$(ii)$]  For each $i > l$, all the neighbors of $v_i$ with indices less than or equal to $l$ has different images with respect to the mapping $f$.
\end{itemize}

Note that the base case is trivial, that is, any partial mapping $f(v_1)$ is enough.
Suppose that the function $f$ satisfies the above properties for all $j \leq t$ where $t \in \{1,2,...,k-1\}$ is  fixed.
Now assume that $v_{t+1}$ has $s$ neighbors with indices greater than $t+1$.
Then $v_{t+1}$ has at most $(\Delta - s)$ neighbors with  indices less than $t+1$.
Let $A$ be the set of neighbors of $v_{t+1}$   with  indices greater than $t+1$.
 Let $B$ be the set of vertices   with indices at most $t$ and with at least one neighbor in $A$.
 Note that as each vertex of $A$ is a neighbor
 of $v_{t+1}$ and has at most $\Delta-1$ neighbors with lesser indices, $|B| = (\Delta-2)|A| = s(\Delta-2)$.
 Let $D$ be the set of possible  options for $f(v_{t+1})$ such that the partial mapping
 is a homomorphism of $G_{t+1}$ to $C$.
 As $C$ has property $Q^{\Delta-1,j}_{1+(\Delta-j)(\Delta-2)}$ we have $|C| \geq 1+s(\Delta-1)$.
 So the set $D \setminus B$ is non-empty.
 Thus, choose any vertex from $D \setminus B$ as the image $f(v_{t+1})$.
 Note that this partial mapping satisfies the required conditions.  \hfill $ \square $

\medskip

Finally, we are ready to prove Theorem~\ref{chromatic-degree}.

\medskip

\noindent \textit{Proof of Theorem~\ref{chromatic-degree}.}  First we will prove the lower bound.
Let $G_t$ be a $\Delta$ regular graph on $t$ vertices. Thus, $G_t$ has $\frac{t \Delta}{2}$  edges. Then
we have $$k_t = \chi_{(m,n)}(G_t)  \geq \frac{p^{\Delta /2}}{p^{\left. {k_t \choose 2} \middle/  t \right.}}$$
using inequality~(\ref{eqn wolog}) (see Section~\ref{sec acyclic}). If $\chi_{(m,n)}(G_t) \geq p^{\Delta /2} $ for some $t$, then we are done.
Otherwise, $\chi_{(m,n)}(G_t) = k_t$ is bounded. In that case, if $t$ is sufficiently large, then $\chi_{(m,n)}(G_t) \geq p^{\Delta /2}$ as
$\chi_{(m,n)}(G_t)$ is a positive integer.

\medskip

Let $G = (V, A \cup E)$ be a connected $(m,n)$-colored mixed graph with maximum degree $\Delta \geq 5$ and $p = 2m+n \geq 2$.
If $G$ has a vertex of degree at most $(\Delta-1)$ then it has degeneracy at most $(\Delta-1)$. In that case  by Theorem~\ref{chromatic-degree}
we are done.

Otherwise, $G$ is $\Delta$ regular. In that case, remove an edge $uv$ of $G$ to obtain the graph $G'$. Note that $G'$ has
maximum degree at most $\Delta$  and has degeneracy at most $(\Delta-1)$. Therefore, by Theorem~\ref{chromatic-degree} there exists an
$(m,n)$-colored complete mixed graph $C$ on
$2 (\Delta-1)^p .p^{(\Delta-1)}$ vertices to which $G'$ admits a $f$ homomorphism to.
Let $G''$ be the graph obtained by deleting the vertices $u$ and $v$ of $G'$. Note that the homomorphism $f$ restricted to $G''$ is
 a homomorphism $f_{res}$ of $G''$ to $C$. Now include two new vertices $u'$ and $v'$ to $C$ and obtain a new graph $C'$.
Color the edges or arcs between the vertices of $C$ and $\{u',v'\}$ in such a way so that we can extend the homomorphism $f_{res}$ to a homomorphism
$f_{ext}$
of $G$ to $C'$  where $f_{ext}(u) = u'$, $f_{ext}(v) = v'$ and $f_{ext}(x) = f_{res}(x)$ for all $x \in V(G) \setminus \{u,v\}$.
It is easy to note that the above mentioned process is possible.

Thus, every connected $(m,n)$-colored mixed graph with maximum degree $\Delta$ admits a homomorphism to $C'$.
\hfill $ \square$

\bibliographystyle{abbrv}