\section{Introduction}

About thirty years ago, the Mannheim-Kazanas (MK) metric has been obtained as an exact exterior solution to conformal Weyl gravity
associated with a static, spherically symmetric gravitational source (Mannheim, 1989). The Schwarzschild term $2M/r$ is accompanied by a linear potential, $\gamma r$, and a quadratic contribution, of the form $- \lambda r^2$.
The linear term is generated
through the effect of cosmology on individual galaxies, while the quadratic one is induced by inhomogeneities in the cosmic background. These inhomogeneities are associated with distances larger than 1 Mpc.
Since the two contributions have opposed signs, the
galaxies are able to support bound orbits.

Even though the MK metric has not been obtained in Einstein's General Relativity, it has the advantage of solving the problem of the shape of the galactic rotation curves (Rubin, 1978), without postulating the dark matter existence.

As it is known, in real galaxies, at large values of $r$, the rotational velocity remains at nearly the same level or
is increasing.
Thus, in order to maintain the high
velocity, an additional contribution to the Newtonian one is required.

In this respect, the MK metric agrees with the empirical expression of the orbiting velocities of the visible matter around the center of galaxy
\[
v^2 = \frac{a}{R} + b R \, ,
\]
where $R$ is the distance to the center of the galaxy, while the parameters $a$ and $b$ contain the number of stars in the galaxy.

In the last years, the MK  metric was tested on a broad sample of galaxies whose rotation curves extend
well beyond the galactic optical disks. The universal quantity $\gamma / \lambda$ gives a natural limit on the size of galaxies (Mannheim, 2011).

The present work is devoted to particles evolving on the Mannheim-Kazanas spacetime.
The effective potential and the conditions to have closed trajectories are discussed.
Working within a tetradic approach, the solutions to the Gordon equation are obtained, for different ranges of the radial variable.

The study of the motions of particles in this metric is very important since it provided a way to detect the presence of a global de Sitter-like
component and to found a specific value for its strength (Mannheim, 2011).

\section{The Mannheim-Kazanas metric and the Galaxy rotation curves}

Let us start with the original expression of the Mannheim-Kazanas (MK) vacuum solution (Mannheim, 1989)
\begin{equation}
g_{00} = 1-3 \beta \gamma - \frac{\beta (2-3 \beta \gamma)}{r} + \gamma r - \lambda r^2\, ,
\end{equation}
where $\gamma$ and $\lambda$ are universal constants, the same for all galaxies,
$\gamma = \gamma_0 = 3 \times 10^{-28} m^{-1}$ and $\lambda = 9.54 \times 10^{-50} m^{-2}$.

Also, $\beta =M$ can be associated to the mass
at the galaxy center, and one may work in the approximation $3 \gamma M \ll 1$.

The metric (1) can be written in the more transparent physical expression
\begin{equation}
g_{00} = 1 -\frac{2M}{r} + \gamma r - \lambda r^2
\end{equation}
and the spherically symmetric line element
\begin{equation}
ds^2= g_{11} (dr)^2 + r^2 \left[  (d \theta )^2 + \sin^2 \theta (d \varphi)^2 \right] - g_{00} (dt)^2 \; ,
\end{equation}
for $g_{11}= g_{00}^{-1}$, is given by
\begin{equation}
ds^2= \frac{(dr)^2}{g_{00}} + r^2 \left[  (d \theta )^2 + \sin^2 \theta (d \varphi)^2 \right] - g_{00} (dt)^2 = -d \tau^2 \; ,
\end{equation}
where $\tau$ is the proper time.

Using the Lagrangian
\[
- L = \frac{1}{2} \left[ \frac{(\dot{r})^2}{g_{00}} + r^2 \dot{\theta}^2 + \sin^2 \theta \dot{\varphi}^2- g_{00} \dot{t}^2 \right]  ,
\]
one can derive the conserved energy and angular momentum
\begin{equation}
g_{00} \dot{t} = E \; , \;
r^2 \sin^2 \theta \dot {\varphi} = - \, K \, ,
\end{equation}
where {\it dot} means the derivatives with respect to $\tau$.

If the motion of the particle with zero angular momentum is on the equatorial plane, i.e. $\dot{\theta} = \dot{ \varphi} =0$, the relations (4) and (5) lead to the important constraint
\[
g_{00} \dot{t}^2 - \frac{\dot{r}^2}{g_{00}} =1 \, ,
\]
i.e.
\begin{equation}
\dot{r}^2 = E^2 - g_{00}  \, .
\end{equation}
In the above expression, one may identify the effective potential
\begin{equation}
V_{eff} = g_{00} = 1 - \frac{2M}{r} + \gamma r - \lambda r^2 \; ,
\end{equation}
which is represented in the figure 1. The potential goes to $- \infty$, for $r \to \pm \infty$ and has a singularity in $r=0$.

\begin{figure}
  \centering
  \includegraphics[width=0.45\textwidth]{Fig1.pdf}
  \caption{The effective potential (7), for $\Delta >0$.}
\end{figure}

The force acting on a particle,
\[
F = - V_{eff}^{\prime} = - \frac{2M}{r^2} - \gamma + 2 \lambda r \, ,
\]
contains besides the attractive contributions, the repulsive term $\lambda r^2$ which becomes important at cosmological distances.

For the metric (2), the horizons are given by the equation $g_{00}=0$, namely they are the solutions of the cubic equation
\begin{equation}
- \lambda  r^3 + \gamma r^2 + r - 2M  =0 \; .
\end{equation}
This has three real roots if
the discriminant $\Delta = 18abcd-4b^3d+b^2c^2-4ac^3-27a^2d^2$ is a positive quantity (Shelbey, 1969),
i.e.
\[
\Delta = -108 M^2 \lambda^2  + 4(1+9 \gamma M) \lambda + \gamma^2(1+8 \gamma M) > 0 \,  .
\]
To first order in $\gamma M \ll 1$ and $\lambda M^2 \ll1$, the three roots of (8) are given by the simple expressions: .
\begin{eqnarray}
r_1 \approx - \frac{1}{\gamma} \left[ 1 - \frac{\lambda}{\gamma^2} \right] ; \; \;
r_2 \approx 2 M \left[ 1 + 4 \lambda M^2 - 2 \gamma M \right] ; \; \; r_3 = \frac{\gamma}{\lambda}  + \frac{1}{\gamma}  \; .
\end{eqnarray}
Thus, there are two positive roots of the equation (8), corresponding to the physical horizons.
The first one is the black hole horizon, $r_h = r_2$, situated in the Schwarzschild region, while the second one corresponds to the cosmological horizon, $r_{\lambda} = r_3$.

The existence of a circular orbit of radius $R_c$ imposes the conditions
$\dot{r}=0$ and $\dot{r}^{\prime} =0$, at $r= R_c = const.$ The corresponding equation
\[
V_{eff}^{\prime} ( r= R_c ) =0 \,  ,
\]
i.e.
\begin{equation}
-2 \lambda R_c^3 + \gamma R_c^2 + 2M =0
\end{equation}
has a negative discriminant $\Delta = -8M (\gamma^3 + 54 \lambda^2 M)$ and therefore we have one real solution
\[
R_c = \frac{1}{6 \lambda} \left[ \gamma + \frac{\gamma^2}{K^{1/3}} + K^{1/3} \right] \; , \; \;
K = \gamma^3 + 108 \lambda^2 M - \left[ 216 \lambda^2 M \left( \gamma^3 + 54 \lambda^2 M \right) \right]^{1/2} .
\]
To first order in $M \lambda^2 / \gamma^3$, the circular orbit's radius is
\begin{equation}
R_c \approx \frac{\gamma}{2 \lambda} + \frac{4M \lambda}{\gamma^2} \, ,
\end{equation}
and it depends, besides the universal constants $\gamma$ and $\lambda$, on the mass $M$.

Thus, a particle with the energy
\[
E^2  = V_{eff} (r= R_c )  = V_{max}
\]
will follow a circular orbit of radius $R_c$, with $R_c$ in between the two horizons.
This circular orbit is unstable because the potential is concave down.

For particles with the energy $E^2 < V_{max}$, there may be two turning points (see figure 2). If one imposes $\dot{r}=0$ in the relation (6), he obtains the equation
\[
E^2 = 1 - \frac{2M}{r} + \gamma r - \lambda r^2 \, .
\]
\begin{figure}
  \centering
  \includegraphics[width=0.45\textwidth]{Fig2.pdf}
  \caption{The particle with the energy $E^2 < V_{max}$ is represented by the dashed line. The two turning point are given by the intersections with the effective potential.}
\end{figure}
The corresponding cubic equation
\begin{equation}
\lambda r^3 - \gamma r^2 + \varepsilon r + 2M =0 \, ,
\end{equation}
with $\varepsilon = E^2 -1$, has three real solution if $\Delta >0$. To first order in $\varepsilon$, this means the energy range
\[
0 \leq \varepsilon < \frac{2 \gamma^3 - 27 \lambda^2 M}{9 \gamma \lambda} \; ,
\]
for $M< 2 \gamma^3 /(27 \lambda^2) \approx 2 \times 10^{14} m \approx 10^{11} M_S$.

The two positive roots of (12), denoted by $r_1$ and $r_2$, are situated in between the two horizons. In the figure 2, these are given by the intersections of the dashed line with the potential. When the energy is increasing, the interval between $r_1$ and $r_2$ becomes smaller and the elliptic trajectory turns into a circle for $E^2 = V_{max}$.
The relations between the important values of $r$ discussed up to now are:
$r_h < r_1 < R_c < r_2 < r_{\lambda}$.

Let us mention that a particle coming from large $r$ values, on her way to the horizon $r_h$, is encountering a Minkowskian region where the terms $M/r$ and $\lambda r^2$ are compensated by $\gamma r$ and therefore $g_{00} \approx 1$. The positive real solutions of the equation $g_{00}=1$, i.e.
\[
- \lambda r^3 + \gamma r^2 - 2M =0  \; ,
\]
whose $\Delta$ is positive for $2 \gamma^3 > 27 \lambda^2 M$,
can be approximated to
\begin{equation}
R_{1} \approx \frac{M \lambda}{\gamma^2} + \sqrt{\frac{2M}{\gamma}}  \, , \; \;
R_{2} \approx \frac{\gamma}{\lambda} - \frac{2M \lambda}{\gamma^2} \, .
\end{equation}
One can check that $R_1 \ll R_2$ and both $R_1$ and $R_2$ are in between the two horizons in (9). Also, for large $r$ values, where the $M/r$ contribution can be neglected, one has only one universal Minkowskian region, at $R = \gamma / \lambda$, which gives the size of galaxies.

Finally, let us discuss the intriguing subject of
orbiting velocities of matter around the galaxy center.

We are using a free of coordinates approach and introduce, for the line element (3), the pseudo-orthonormal frame ${E_a}_{(a=\overline{1,4})}$
\begin{equation}
E_1=  \sqrt{g_{00}} \, \partial_r \; , \; \;
E_2= \frac{1}{r}  \, \partial_{\theta} \; ,
\; \;
E_3=\frac{1}{r \sin \theta}\, \partial_{\varphi} \; , \; \;
E_4= \frac{1}{\sqrt{g_{00}}} \, \partial_t \; ,
\end{equation}
whose corresponding dual base is
\begin{equation}
\omega^1=  \frac{dr}{\sqrt{g_{00}}} \; , \; \;
\omega^2=  r  \,d\theta \; ,
\; \;
\omega^3= r \sin \theta \, d\varphi \; , \; \;
\omega^4= \sqrt{g_{00}} \,dt \; ,
\end{equation}
so that $ds^2 = \eta_{ab} \omega^a \omega^b$, with $\eta_{ab} = {\rm diag} [ 1 , 1 , 1 , -1]$.

Using the first Cartan's equation,
\begin{eqnarray}
& &
d\omega^a=\Gamma^a_{.[bc]}\,\omega^b\wedge \omega^c \;  ,
\end{eqnarray}
with $1 \leq b<c \leq4$ and $\Gamma^a_{.[bc]} = \Gamma^a_{.bc} - \Gamma^a_{. cb}$,
we obtain the connection coefficients
\begin{eqnarray}
& &
\Gamma_{122} =  \Gamma_{133} = - \frac{ \sqrt{g_{00}}}{r} = - \Gamma_{212}  = - \Gamma_{313}
\nonumber \\*
& &
\Gamma_{233} = - \Gamma_{323} = - \, \frac{\cot \theta}{r}  \; ,
\nonumber \\*
& &
\Gamma_{144} = - \Gamma_{414} = \frac{g_{00}^{\prime}}{2 \sqrt{g_{00}}} \; .
\end{eqnarray}

The equations for the timelike geodesics
\[
\frac{du^a}{d \tau} + \Gamma^a_{\; bc} u^b u^c = 0 \;  ,
\]
where the four velocity components,
\[
u^a = \frac{\omega^a}{d \tau} \, ,
\]
are given by:
\[
u^1 = \frac{\dot{r}}{\sqrt{g_{00}}} \; , \; \; u^2 = r \dot{\theta} \; , \; \; u^3 = r \sin \theta \dot{\varphi} \; , \; \;
u^4 = \sqrt{g_{00}} \dot{t} \; ,
\]
have the explicit form
\begin{eqnarray}
& &
\ddot{r} - \frac{g_{00}^{\prime}}{2 g_{00}} \dot{r}^2 - r g_{00} \left[ \dot{\theta}^2 + \sin^2 \theta \dot{\varphi}^2 \right] + \frac{g_{00} g_{00}^{\prime}}{2} \, \dot{t}^2 = 0  \; ,
\nonumber \\*
& & r \ddot{\theta} + 2 \dot{r} \dot { \theta} - r \sin \theta \cos \theta \dot{\varphi}^2 = 0 \; , \nonumber \\*
& & r \sin \theta \ddot{\varphi} + 2 \sin \theta \dot{r} \dot{ \varphi} + 2 r \cos \theta \dot{\theta} \dot{ \varphi} = 0 \, ,
\nonumber \\*
& & \ddot{t} + \frac{g_{00}^{\prime}}{g_{00}} \, \dot{r} \dot{t} = 0  \;  .
\end{eqnarray}
For the circular orbit in the equatorial plane, i.e.
$\dot{r}=0$, $\theta = \pi/2$, $\dot{\theta} =0$, the first relation in (18) becomes
\[
\frac{d \varphi}{dt} = \sqrt{\frac{g_{00}^{\prime}}{2r}} \,  ,
\]
leading to the velocity
\begin{equation}
v = r \frac{d \varphi}{dt} = \sqrt{\frac{r g_{00}^{\prime}}{2}} \,  .
\end{equation}
Obviously, in the Schwarzschild case, one obtains the well known Newtonian expression
\[
v = \sqrt{\frac{M}{r}} \, ,
\]
while for the metric (2), the velocity is given by
\begin{equation}
v = \sqrt{\frac{1}{2} \left[ \frac{2M}{r} + \gamma r - 2 \lambda r^2 \right]} .
\end{equation}
For $2 \lambda r \ll \gamma$, when only the first two terms in the parenthesis are taken into account, one may notice that the velocity is significantly increased, compared to the Schwarzschild  expression.

The relation $v^2 \geq 0$ is satisfied by $r \leq R_c$, where $R_c$ is the circular radius (11).

\section{The Gordon equation and its Heun solutions}

Working in a free of coordinates approach (Dariescu, 2017), based on Cartan's equations, the
Einstein tensor components in the tetradic frame (14) are:
\begin{eqnarray}
& &
G_{11} = \frac{2 \gamma}{r} - 3 \lambda \; , \nonumber \\*
& &
G_{22} = G_{33} = \frac{\gamma}{r} - 3 \lambda \;   , \nonumber  \\*
& & G_{44} = 3 \lambda - \frac{2\gamma}{r} \; .
\end{eqnarray}
One may notice that, from the point of view of Einstein's General Relativity, the proper energy density $\rho = T_{44} =G_{44}$ is negative for $r< ( 2 \gamma)/(3 \lambda)$, leading to the conclusion that, at values of $r$ smaller than $2 \gamma / ( 3 \lambda)$, the astrophysical object described by the MK metric is surrounded by exotic matter.

The scalar curvature,
\[
R = 12 \lambda - \, \frac{6\gamma}{r} \, ,
\]
has a singularity in $r=0$ and is positive for $r> \gamma/(2 \lambda)$.

In the general expression of the Klein-Gordon equation (Dariescu, 2017)
\begin{eqnarray}
& &
\frac{1}{r^2} \frac{\partial \;}{\partial r} \left[ r^2 g_{00} \frac{\partial \Phi}{\partial r} \right]
+ \frac{1}{r^2 \sin \theta} \frac{\partial \;}{\partial \theta} \left[ \sin \theta \frac{\partial \Phi}{\partial \theta} \right]
+ \frac{1}{r^2 \sin^2 \theta} \frac{\partial^2 \Phi}{\partial \varphi^2}
\nonumber \\*
& & - \;
\frac{1}{g_{00}} \frac{\partial^2 \Phi}{\partial t^2} - \mu^2 \Phi \, = 0 \; ,
\end{eqnarray}
where $g_{00}$ given in (2) is depending only on $r$, one can perform the variables separation
\begin{equation}
\Phi = F(r) Y_{\ell}^m ( \theta , \varphi ) \, e^{-i \omega t} \, ,
\end{equation}
where $Y_{\ell}^m$ are the spherical functions. The corresponding radial equation,
\begin{eqnarray}
\frac{1}{r^2} \frac{d \;}{dr} \left[ r^2 g_{00} \frac{dF}{dr} \right]
+ \left[ \frac{\omega^2}{g_{00}} - \frac{\ell ( \ell+1) }{r^2} - \mu^2 \right] F = 0 \, ,
\end{eqnarray}
can be analytically solved only in particular cases.

Let us consider the very large distances, typical to clusters and superclusters, for which the metric function can be approximated to
\begin{equation}
g_{00} \approx 1 +\gamma r - \lambda r^2  = - \lambda \left[ r + \frac{b-\gamma}{2\lambda} \right] \left[ r - \frac{b+ \gamma}{2 \lambda} \right] ,
\end{equation}
where $b = \sqrt{\gamma^2 + 4 \lambda}$.
There is an unique horizon, which is the positive root of the equation $g_{00}=0$, i.e.
\begin{equation}
r_h = \frac{b+\gamma}{2 \lambda}  \; ,
\end{equation}
which can be understood as the boundary of this spacetime.
Using the numerical values of $\gamma = \gamma_0$ and $\lambda$, we get $r_h \approx 100 Mpc$.

The radial function, solution to (24) with $\mu=0$, is expressed in terms of Heun general functions (Ronveaux, 1995; Slavyanov, 2000), as
\begin{eqnarray}
F(r)  & = & \left[ r -  \frac{b + \gamma}{2 \lambda} \right]^{-4i \omega \lambda/b}  \left[ r +  \frac{b - \gamma}{2 \lambda} \right]^{4i \omega \lambda/b}  r^{- \frac{1 \pm \sqrt{1+16 \lambda \ell ( \ell +1)}}{2}}  \nonumber \\* & \times &
HeunG \left[ {\rm a}  \, , q , \, \alpha , \, \beta , \, \gamma , \, \delta ,  - ( \lambda r_h ) r  \right] .
\end{eqnarray}

As it is known from the theory (Slavyanov, 2000), the Heun general function of variable $x$ has four singular points and one may check that $r= r_h$ is one of them.

As we turn to regions well inside a galactic halo boundary, the $\lambda$ contribution can be neglected and the term $2M/r$ becomes important.
The metric (2) can be approximated to
\begin{equation}
g_{00} = 1- \frac{2M}{r} + \gamma r = \frac{\gamma}{r} \left[ r +  \frac{a+1}{2\gamma} \right]  \left[ r -  \frac{a-1}{2\gamma} \right] ,
\end{equation}
where
\[
a = \sqrt{1+8 \gamma M}   \, .
\]
This case has been investigated, in detail, in (Dariescu, in press).
Similarly to the previous case, there is one horizon, close to the Schwarzschild radius,
\begin{equation}
r_h = \frac{a-1}{2 \gamma} \approx 2M ( 1-2 \gamma M) \; .
\end{equation}

In the massless case, the radial function is expressed in terms of Heun general functions (Ronveaux, 1995; Slavyanov, 2000), as
\begin{eqnarray}
F(r)  =  \left[ r -  \frac{a-1}{2 \gamma} \right]^{\frac{ i \omega (a-1)}{2 \gamma a}}\left[ r +  \frac{a+1}{2 \gamma} \right]^{- \, \frac{ i \omega (a+1)}{2 \gamma a}}
HeunG \left[ {\rm a} \, , q , \, \alpha , \, \beta , \, \gamma , \, \delta , 1 + \frac{2 \gamma r}{a+1} \right] .
\end{eqnarray}

The function (30) is much more complex than the one obtained for the Schwarzschild metric
\[
g_{00} = 1 - \frac{2M}{r} \, ,
\]
where the amplitude function expressed in terms of the Heun confluent functions as
\begin{equation}
F_S = e^{i \omega r} ( r-2M)^{\mp 2i \omega M} HeunC \left[ - 4i \omega M , \mp 4 i \omega M, 0 , - 8 \omega^2 M^2 , 8 \omega^2 M^2 - \ell (\ell+1), 1 - \frac{r}{2M} \right] .
\end{equation}

Compared to the Heun general functions which have four singular points, the Heun confluent functions have two regular and one irregular singularities and these can be obtained from the Heun general functions by a confluence process, when two of the singularities coalesce.
The Heun confluent function can be computed as a power series expansion around the origin $z=0$, i.e. $r=2M$, and the series converges for $|z| <1$,
 (Ronveaux, 1995; Slavyanov, 2000).

The horizon defined in (29) is closer to the $r=0$ singularity compared to the Schwarszchild horizon, $r_S =2M$.
On her way to the horizon $r_h$, the particle is crossing the Minkowskian region $r=R = \sqrt{2M/\gamma}$, where $g_{00} \approx 1$.

One may notice that, when $r \to 0$, the amplitude function (30) goes to infinity since the Heun general function, of variable $x$, has a singularity in $x=1$.
However, we should not worry about that because the metric (2) is valid only on the external region of an astrophysical object of radius $R$.

\section{Conclusions}

In this work, we have considered a spacetime described by the metric proposed by Mannheim and Kazanas (MK metric) (Mannheim, 1989), which contains two parameters, $\gamma$ and $\lambda$. We have worked with the universal values: $\gamma = \gamma_0 = 3 \times 10^{-28} m^{-1}$ and $\lambda = 9.54 \times 10^{-50} m^{-2}$, so that $\gamma / \lambda \approx 100 kpc$ (Mannheim, 2011). Obviously, this value of $\gamma$ is independent of the galactic mass.

At a first site, the MK metric can be seen as an extension of the Schwarzschild solution. However, one can notice significant differences between these two.
Thus, for the metric (2), the horizons are given by the roots of the cubic equation (8) and have the approximate expressions given in (9). One can identify a Minkowskian region, in between the two horizons, where $g_{00} \approx 1$, which has no analogue to the Schwarzschild case. The constant potential is leading to approximately constant rotational velocities which agree with the observed flat galactic rotation curves (Rubin, 1978).

For a typical galaxy with $M = 10^{11} M_s$, so that
\[
M = \frac{10^{11} M_S G}{c^2} \approx 1.5 \times 10^{14} m \, ,
\]
the terms $2M/r$ and $\gamma r$ are almost equal for $r_* = \sqrt{2M / \gamma } \approx 30 kpc$, which corresponds to the Minkowskian region, once we neglect the $\lambda r^2$ contribution.

Also, a particle moving in the potential (7), with suitable energy values, can follow an elliptic or a circular orbit. The radius (11) has been calculated in the approximation $M \lambda^2 / \gamma^3 \approx 3M \times 10^{-16} \ll 1$. The universal quantity $\gamma / \lambda$ has been seen as a natural geometric limit on the size of galaxies, in the sense that beyond this distance there could no longer
be any bound galactic orbits (Mannheim, 2011). In addition, a very important point has been made in (Nandi, 2012), with respect to the actual size of the galaxies, which must be related to the radius of the last circular stable orbit.

The velocity expression (20), derived from the geodesics equation, agrees with the one obtained by O'Brien and Mannheim in the framework of conformal gravity theory (O'Brien, 2012).
As it is known, the Newtonian velocity falls below data as the orbital radius increases. The presence of the linear term is leading to an excess in the observed velocities.
With no free parameters
other than the galactic mass, O'Brien and Mannheim proved that their result can fit an impressive number of galaxies.

However, for fitting the observed orbital velocities for more than 100 galaxies, $\gamma$ has been taken as $\gamma = \gamma_0 + \gamma_G$, where $\gamma_G = N \gamma^*$, with $N = M/M_S$ and $\gamma^* = 5.42 \times 10^{-39} m^{-1}$ (Mannheim, 2012).
Thus, besides the global cosmological linear term, there is a second one, which is depending on the matter within the galaxy.
Such an analysis is beyond the purpose of our paper and for a detailed discussion on this subject, we recommend (Nesbet, 2018).

In the final part of the paper, we have worked out the Gordon equation for massless bosons evolving in the MK spacetime. The radial equation is considered for different regions of the variable $r$, starting with distance scales between $1 Mpc$ and
$100 Mpc$, where $g_{00}$ can be taken as in (25), down to regions where the $2M/r$ term has a significant contribution and the term $\lambda r^2$ can be neglected.

The whole analysis is much more involved compared to the Schwarzschild exterior metric, the solutions to the Gordon equation being given by the Heun general functions, (Slavyanov, 2000).


\section{Introduction}
Some years ago, Tan derived a set of exact relations that link the short distance large-momentum correlations to the bulk thermodynamic properties of a fermionic system with short-range inter-particle interactions \cite{Tan1, Tan2, Tan3}. They are connected by a single coefficient $C$, referred to as the integrated contact intensity or ``contact". Further, Tan's ideas were developed \cite{Br8, Br9} and extended to Bose systems also
\cite{Br10,Lang11,Werner1,Combescot13}. It has been established that
Tan's contact measures the density of pairs at short distances and determines the exact large-momentum
or high frequency behavior of various physical observables. It serves as an important
quantity to characterize strongly interacting many-body systems \cite{Pitbook14}.

Experimentally, Tan's relations have been confirmed both for fermions (especially for quasi-1D systems) \cite{kuhle2011, sagi2012, stewart2010} and Bose gases \cite{chang2016, Fletcher2016, makotyu2014, wild2012, zou2021}.

Although, there are several theoretical studies of contact of low dimensional Fermi systems even at
finite temperatures \cite{Hoffman,yuPRA80,Matveeva,Br8,brateen2013,vigndofermi,Patu}, the temperature dependence of bosonic contact in three dimensions is mostly unknown.

Studying the temperature dependence of contact parameter gives an opportunity to know its critical behavior near phase transitions. For example, Chen {\it et al}. \cite{chencrit} have shown {that contact and its derivatives are uniquely determined by the universality class of the phase transition for fermions in the critical region.}
In a system of bosons, critical behavior of thermodynamic quantities have some specifics due to Bose-Einstein condensation (BEC) at low temperatures. For instance, the specific heat is discontinuous at the critical temperature $T_c$, and the Gr{\"u}neisen parameter changes its sign \cite{garst,ourMCE}.

The main goal of the present work is to study the temperature dependence of Tan's contact for bosons, as well as its derivatives, to draw conclusions about their critical behavior near $T_c$  and close to the quantum critical
point (QCP) of the quantum phase transition (QPT).

It is well known that physical systems at equilibrium are studied in statistical mechanics through statistical
ensembles. For example, a system {that} the canonical ensemble exchanges only energy with the surroundings
is called $NVT$, where $N$ is the number of atoms of a gas in the volume $V$ at temperature $T$.
On the other hand,
a grand canonical ensemble {also makes it possible for particles to exchange},
with the parameter $\mu$, which may be
interpreted as the chemical potential, and is called $\mu VT$. In the former case, $N$ or $\rho=N/V$ is fixed
and the chemical potential is determined by the equation of state, while in the latter case $\mu$ is fixed as an
input parameter which defines number of particles through a thermodynamic equation \cite{Marzolino}.
Tan's relations including the contact have been thoroughly studied for $NVT$ ensembles, {however, the investigations of contact for a $\mu VT$ ensemble is still missing. In this work, we choose
a system} of bosonic quasiparticles introduced in bond operator formalism \cite{Sachdev,OurAnn2} to describe
the singlet-triplet excitations in spin gapped magnetic materials \cite{Zapf}. It is well established that low temperature
properties of such a class of quantum magnets could be described within the paradigm of Bose-Einstein
condensation of these quasiparticles {called} triplons. Therefore,  one may conclude that,  at low
temperature, thermodynamic properties of such materials are {mainly determined} (but not only) by the
condensation (and depletion) of triplons \cite{Zapf,ourMCE,ourJT,ourCharac}. Triplon concept has found
application in other spin-gapped models as well \cite{kumar}.

There is one more reason of our choice of spin gapped quantum magnets to study Tan's contact.
In such systems, the number of particles $N$ and the condensate fraction $N_0/N$ can be directly evaluated by measuring the uniform magnetization $M$ and staggered magnetization $M_{\perp}$ respectively,
as $N=M/g\mu_B$, $N_0^2=2M_{\perp}^2/(g\mu_B)^2$
\cite{delamorre}, where $g$ is the electron Land{\'e} factor and $\mu_B$ is the  Bohr magneton.
Evidently, the final analytic expression for $C$ will include $M$ and $M_{\perp}$, making the contact parameter
to be easily measured.

In practice, Tan's contact for bosons with zero range interaction may be evaluated by any of following Tan's relations \cite{Pitbook14}. {For example,}
\begin{enumerate}
\item By the asymptotic behavior of momentum distribution $n_k$,
\begin{equation}
C=\lim_{k\rightarrow \infty}k^4 n_k\, .
\end{equation}\normalsize
\item By Tan's sweep theorem as
\begin{equation}
C=\frac{8\pi ma^2}{V} \
\begin{cases}
\left(\displaystyle\frac{\partial F}{\partial a}\right)_{N,T} \ \ \ & \text{for $NVT$ ensembles}\, , \\
\left(\displaystyle\frac{\partial\Omega}{\partial a}\right)_{\mu,T} \ \ \ & \text{for ${\mu}VT $ensembles}\, ,
\end{cases}
\label{Cdef2}
\end{equation}\normalsize
where  $F$ and $\Omega$ are the free energy and the grand canonical potential, respectively \cite{bouchul21}; $a$ is the scattering length and $m$ is the effective mass of the particle.
\item By using Hellman-Feynman theorem \cite{Br8}
\begin{equation}
C=\frac{16\pi^2 a^2}{V}\int d\vec{r}\langle\psi^{\dagger}(\vec{r})\psi^{\dagger}(\vec{r})\psi(\vec{r})\psi(\vec{r})\rangle\, .
\end{equation}\normalsize
\end{enumerate}
In our previous paper \cite{ourTan1}, it has been shown that if one uses mean field theory (MFT), evaluation of $C$
from Eq.\,(\ref{Cdef2}) will be the most convenient and reliable. {Consequently, we study $C$}
of a triplon gas at finite temperature in the framework of MFT, namely in the Hartree-Fock-Bogoliubov (HFB)
approximation \cite{ourTan13}.
{We stress here that below we shall discuss $C$ only for homogeneous systems. For inhomogeneous systems,
especially at very low temperature, one may also use the Gross-Pitaevskii equation to study the
dynamics of the condensate \cite{wangpra81,wangpra84} as well as  that of the Tan's contact.}

{For convenience, we adopt units such that $\hbar=1$, $k_B=1$ and $V=1$ in the following text}. In these units
$C$ is dimensionless.
In natural units $C$ may be obtained further by dividing it by ${\bar a}^4$, getting it in
(length)$^{-4}$ where ${\bar a}$
is the average lattice parameter, ${\bar a}=V^{1/3}$, and $V$ is the unit cell volume of the crystal.

{This article is organized as follows}. In Section II, we outline main equations for a triplon gas at finite temperature in the HFB approximation, which will be used to evaluate the Tan's contact and its derivatives
in Section III. Our discussions and conclusions will be presented in Section IV.

\section{Triplon density in the HFB approximation}

We start with the effective Hamiltonian of triplons as a non-ideal Bose gas with contact repulsive interaction
\begin{equation}
{\ds\hat H}=\int d\vec{r}\left[\psi^{\dagger}(\vec{r})(\hat{K}-\mu)\psi(\vec{r})+\frac{U}{2}(\psi^{\dagger}(\vec{r}) \psi(\vec{r}))^2 \right]
\label{Ham}
\end{equation}\normalsize
where $\psi(\vec{r})$ is the bosonic field operator, $U=4\pi a/m$ is the interaction strength,
and $\hat{K}$ is the kinetic energy operator which defines the bare triplon dispersion $\varepsilon_{\vec{k}}$
in momentum space. Since the triplon BEC occurs in solids, the integration is performed over the unit cell of the
crystal with the corresponding momenta defined in the first Brillouin zone \cite{ouraniz1}. The parameter $\mu$ characterizes an additional direct contribution to the triplon energy due to the external field
\begin{equation}
\mu=g \mu_B H-\Delta
\label{chempot}
\end{equation}\normalsize
and it can be interpreted as the chemical potential of the $S_z=-1$ triplons
\cite{Giamarchi}.
In Eq.\,(\ref{chempot}),  $\Delta=g\mu_B H_c$,
is the spin gap separating the ground-state from the lowest energy triplet excitation
and $H_c$ is the critical external magnetic field \cite{grundmann}. {The effective
Hamiltonian} ${\hat H}$ has the following
free parameters characterizing a given material: $g$, $H_c$, $U$ and possibly, an effective triplon mass $m$
for the case of a simple bare dispersion given as $\varepsilon _{\vec k}={\vec k}^2/2m$.\footnote{
In general, a realistic  $\varepsilon _{\vec k}$ includes several other parameters \cite{misguich}.}

In general the Hamiltonian in Eq.\,(\ref{Ham}) is invariant under global $U(1)$ gauge transformation $\psi(\vec{r})\rightarrow e^{i\alpha}\psi(\vec{r})$ with {a real number $\alpha$}. However, this symmetry is broken in the condensate phase, where $T<T_c$, and it is restored in the normal phase, $T\geqslant T_c$.

This transition temperature $T_c$, which corresponds to the vanishing of the condensate density
$\rho_0(T_c)=0$ may be calculated from the following equation
\begin{equation}
\sum_k \frac{1}{e^{\varepsilon_k/T_c}-1}=\frac{\mu}{2U}\equiv \rho_c\, .
\label{Tc}
\end{equation}\normalsize
Note that for existing quantum magnets with a spin gap $T_c$ is rather large, $T_c\approx 2$\,K, in contrast to the critical temperature of BEC in atomic gases, where the number of particles is fixed
and the critical temperature is of the order of nanokelvins. However, in the $\mu VT$ ensemble
under discussion, the particle density depends on the external magnetic field through the chemical potential $\mu$
as
\begin{equation}
\begin{array}{l}
\rho(T \ge T_c)=\ds\sum_k \frac{1}{e^{\beta\omega_k}-1}, \quad \quad
\omega_k=\varepsilon_k-\mu+2U\rho\, ,
\label{rhobig}
\end{array}
\end{equation}\normalsize
in the normal phase. In the BEC phase, an explicit expression for $\rho(T<T_c)$ strongly depends on the chosen
version of an approximation. For example, in the HFB approximation, which is employed here, the particle
density and condensate fraction are given by the following set of equations \cite{ourAnn}
\begin{eqnarray}
\begin{aligned}
 \rho(T < T_c)=\rho_0+\rho_1=\frac{\mu+\mu_{eff}}{2U}\, , \\
 \mu_{eff}=\mu+2U(\sigma-\rho_1)\, , \\
 \rho_1=\frac{1}{2}\sum_k\left[\frac{\coth(\beta E_k/2)(\varepsilon_k+\mu_{eff})}
{E_k}-1\right]\, , \\
 \sigma=-\frac{\mu_{eff}}{2}\sum_k\frac{\coth(\beta E_k/2)}{E_k}\, ,
\label{MainEqs}
\end{aligned}
\end{eqnarray}
where $\rho_1$ and $\sigma$ correspond to the density of non-condensed particles and anomalous
density, respectively.
The hyperbolic function in Eq.\,\re{MainEqs} can be represented also as  $\coth(\beta x/2)=1/2+f_B(x)$,
where $f_B(x)=1/(\exp(\beta x)-1)$ is the Bose distribution function.

Note that neglecting $\sigma$ leads to an unexpected cusp in magnetization \cite{ourPRB81}.
In Eq.\,(\ref{MainEqs}), $E_k$ is the dispersion of quasiparticles (bogolons) given as $E_k=\sqrt{\varepsilon_k}
\sqrt{\varepsilon_k+2\mu_{eff}}$ with the speed of sound $c=\sqrt{\mu_{eff}/m}$, where $m$ has the meaning of the triplon effective mass, corresponding to  the limit of small momenta $\varepsilon_k \approx \vec{k}^2/2m$. The typical value of $m$, used in the literature,
\cite{Zapf} is  rather small, $m\approx 0.02$\,K.

\section{Tan's contact for a triplon gas}

As it is pointed out in the Introduction, for the $\mu VT$ ensemble it is convenient to calculate $C$ using
\begin{equation}
C=2U^{2}m^{2}\left(\frac{\partial\Omega}{\partial U}\right)_{T,\mu}\, ,
\label{CDeRho}
\end{equation}\normalsize
where the grand thermodynamic potential $\Omega$ has the following total derivative \cite{Pitbook14,ourJT,ourCharac}
\begin{equation}
d\Omega=-SdT-PdV-Nd\mu-MdH+\frac{C da}{8\pi m a^2}\, ,
\label{domega}
\end{equation}
in which $S$ is the entropy and $P$ is the pressure.
Note that the relation between the magnetization $M$ and the triplon density may be directly
obtained from Eq.\,(\ref{domega}) as $M=g\mu_B\rho$.
Moreover, using Eq.\,(\ref{domega}) one may find useful expressions for the derivatives of $C$ as
\begin{eqnarray}
\begin{aligned}
  \left(\frac{\partial C}{\partial T}\right)=-2U^2m^2\left(\frac{\partial S}{\partial U}\right)\, , \\
  \left (\frac{\partial C}{\partial \mu}\right )=
\displaystyle\frac{1}{g\mu_B}\left(\frac{\partial C}{\partial H}\right) =
-2U^2m^2\left(\frac{\partial \rho}{\partial U}\right)\, .
\label{cdef}
\end{aligned}
\end{eqnarray}

It is well known that {from BEC to a normal phase transition} is a second order phase transition,
in which the entropy $S(T,H,a)$ is continuous across the critical temperature i.e.,
\begin{equation}
dS(T=T_c^-)=dS(T=T_c^+)=\left ( \displaystyle\frac{\partial S}{\partial T}\right )_{H,a}dT+
\left ( \displaystyle\frac{\partial S}{\partial H}\right )_{T,a}dH+\left ( \displaystyle\frac{\partial S}{\partial a}\right )_{H,T}da\, .
\lab{ds}
\end{equation}\normalsize
Therefore, using Eqs.\,\re{domega}-\re{ds} leads to an extended Ehrenfest relation:
\begin{equation}
T_{c}^{-1}\Delta  [C_H]=-\Delta\left[\left ( \displaystyle\frac{\partial M}{\partial T}\right )_{H,U}\right]\left (\displaystyle\frac {dH}{dT}
\right )+
\displaystyle\frac{1}{2U^2m^2}\Delta\left[\left ( \displaystyle\frac{\partial C}{\partial T}\right )_{H,U}\right]\left (\displaystyle\frac{dU}{dT}\right )
\lab{eren}
\end{equation}\normalsize
where $C_H=T(\partial S/\partial T)$ is the heat capacity and $\Delta [f]\equiv f(T=T_c^-)-f(T=T_c^+)$ is the
jump in the function $f$ at $T=T_c$.\footnote{A similar relation for $NVT$ ensembles will be
$T_{c}^{-1}\Delta [C_P]=-\Delta\left[\left ( \displaystyle\frac{\partial P}{\partial T}\right )_{V,U}\right]
\left (\displaystyle\frac{dV}{dT}\right )+
\displaystyle\frac{1}{2U^2m^2}\Delta\left[\left ( \displaystyle\frac{\partial C}{\partial T}\right )_{V,U}\right]\left (\displaystyle\frac{dU}{dT}\right )$
when the entropy is considered as a function of temperature, volume and the strength of the contact interaction.}

In the following we discuss the normal and BEC phases, separately.

\subsection{Normal phase ($T\ge T_c$)}

In the normal phase $\Omega$ is given by
\begin{equation}
\Omega(T \ge T_c)=-U\rho^2+T\sum_k \ln(1-e^{\beta\omega_k}),
\label{omegaBig}
\end{equation}\normalsize
where $\rho$ and $\omega_k$ are defined in Eq.\,(\ref{rhobig}).
Now, taking the derivative of $\Omega$ and using Eq.\,(\ref{rhobig}) one may easily find
\begin{equation}
\left (\frac{\partial\Omega}{\partial U} \right )_{T \ge T_c}=-\rho^2-2U\rho\rho'+(2\rho+2U\rho')\sum_k\frac{1}{e^{\beta\omega_k}-1}=\rho^2\, ,
\label{difomUBg}
\end{equation}\normalsize
and hence
\begin{equation}
C(T \ge T_c)=2U^{2}m^{2}\rho^{2}=\frac{2U^{2}m^{2}M^{2}}{(g\mu_{B})^{2}}\, .
\label{Cbig}
\end{equation}\normalsize
The derivatives of $C$ can be evaluated using Eq.\,(\ref{cdef}) or directly from Eq.\,\re{Cbig} as
\begin{equation}
\begin{aligned}
\left (\frac{\partial C}{\partial T}\right )_{T\ge T_c}=4U^2m^2\rho\left ( \frac{\partial\rho}{\partial T}
\right )
=-2m^2U^2\left(\frac{\partial S}{\partial U} \right)\, , \\
\left (\frac{\partial C}{\partial\mu}\right )_{T \ge T_c}=4U^2m^2\rho\left ( \frac{\partial\rho}{\partial\mu}
\right ) =-2m^2U^2\left(\frac{\partial\rho}{\partial U} \right)\, ,
\label{difCTCmu}
\end{aligned}
\end{equation}\normalsize
where \cite{ourMCE}

\begin{equation}
\begin{aligned}
\left ( \frac{\partial S}{\partial U}\right )_{T \ge  T_c}=\frac{2\rho\beta I_1}{1-2I_2}, \ \ \
\left ( \frac{\partial\rho}{\partial U}\right )_{T\ge T_c}=\frac{2\rho I_2}{U[1-2I_2]}, \\
 \\
\label{drodSUBg}
\end{aligned}
\end{equation}\normalsize
and  the integrals $I_1$ are $I_2$ {are} given in the Appendix.
Note that Eqs.\,(\ref{difCTCmu}) lead to the following relations for the susceptibility, $\chi=(\partial M/\partial H)$
in the normal phase
\begin{equation}
\begin{aligned}
\left ( \frac{\partial C}{\partial H}\right )_{T \ge T_c}=\frac{4U^2m^2}{(g\mu_B)^2} \left ( M \chi \right )_{T \ge T_c}\, , \\
\chi ({T \ge T_c})=-\frac{(g\mu_B)^2}{2M}\left(\frac{\partial M}{\partial U}\right) _{T \ge T_c}.
\label{chi1}
\end{aligned}
\end{equation}\normalsize

\subsection{Condensed phase ($T<T_c$)}
	In the condensed phase the thermodynamic potential is given by \cite{ourJT,ourCharac,OurAnn2}
\begin{equation}
\Omega(T<T_c)=\frac{U\rho_0^2}{2}-\mu\rho_0-\frac{U}{2}(2\rho_1^2+\sigma^2)+\frac{1}{2}\sum_k(E_k-\varepsilon_k-\mu_{eff})+T\sum_k \ln(1-e^{-\beta E_k})
\label{omegBEC}
\end{equation}\normalsize
with $\mu_{eff}$, $\rho_1$ and $\rho$ are given by Eqs.\,(\ref{MainEqs}), $\rho_0=\rho-\rho_1=\mu_{eff}/U-\sigma$ and $E_k=\sqrt{\varepsilon_k}\sqrt{\varepsilon_k+2\mu_{eff}}$. For simplicity, one may start with the well-known Bogoliubov approximation, where the depletion and anomalous density are completely neglected. For atomic gases, this leads to the Gross-Pitaevskii equation which is successfully used in a lot of systems \cite{Pitbook14}. In the present case Bogoliubov approximation corresponds to $\rho_0\approx\rho$, $\mu_{eff}\approx\mu=U\rho$ and hence
\begin{equation}
\Omega_{\rm Bogoliubov}(T<T_c)=-\frac{\mu^2}{2U}+\frac{1}{2}\sum_k[E_k-\varepsilon_k-\mu]+T\sum_k \ln(1-e^{-E_k\beta})\, .
\label{Ombog}
\end{equation}\normalsize
Since in $\mu VT$ ensembles the chemical potential does not depend on the inter-particle interaction, i.e.,
$(\partial\mu/ \partial U)=0$, one immediately obtains from Eq.\,(\ref{Ombog}) Tan's contact in the Bogoliubov
approximation as\footnote{In natural units $C_{0}=m^2\mu^2/\bar{a}^4$, where
 the average lattice parameter $\bar{a}$ can be taken from Table 1.}
\begin{equation}
C_{0}=m^2\mu^2\, .
\label{Cbog}
\end{equation}\normalsize
From this equation, it is seen that in this simple approximation Tan's contact does not depend on temperature,
which means that in $\mu VT$ ensembles Bogoliubov approximation is not applicable to finite temperatures.
The quality of this approximation at $T=0$ is tested below.
In Figs.\,1, Tan's contact as a function of the magnetic field is presented for two different compounds
TlCuCl$_3$ (Fig.\,1a) and Ba$_3$Cr$_2$O$_8$ (Fig.\,1b) in natural units, i.e., in \AA$^{-4}$. It is seen that for both compounds, the difference between the exact HFB (solid lines) and simple Bogoliubov approximation
(dashed lines) is rather small, especially for small chemical potentials. Remarkably, the contact given
by Eq.\,(\ref{Cbog}) does not include the scattering length of triplon interaction but includes only the effective
mass as an internal microscopic parameter of the system. This is why the contact for TlCuCl$_3$ at the
same magnetic field e.g., at $H=13.5$\,T, is much smaller than that for Ba$_3$Cr$_2$O$_8$. In fact,
as it is seen from Table 1 the ratio of effective masses is of the order of
$m$({Ba$_3$Cr$_2$O$_8$})/$m$({TlCuCl$_3$})$\approx 10$.

\begin{figure}[h!]
\begin{minipage}[H]{0.49\linewidth}
\center{\includegraphics[width=1.25\linewidth]{Fig1a.pdf} \\ a)}
\end{minipage}
\hfill
\begin{minipage}[H]{0.49\linewidth}
\center{\includegraphics[width=1.25\linewidth]{Fig1b.pdf} \\ b)}
\end{minipage}
\caption{Tan's contact at zero temperature for TlCuCl$_3$ (a) and Ba$_3$Cr$_2$O$_8$ (b) in \AA$^{-4}$
in HFB (solid lines) and Bogoliubov  approximations  (dashed lines). The input parameters are given in Table 1.}
\label{Fig1}
\end{figure}

\begin{table}[h!
\caption {Material parameters used in our numerical calculations. The parameters $g$, $H_{c}$, ${\bar a}$ and $\Delta$ are taken directly from the experimental measurements, while $m$ and $U$ are optimized to fit experimental magnetizations (see Ref. [\onlinecite{ourMCE}] for details)}.
\begin{tabular}{ p{3cm} c c c c c c}
 \hline
 \hline
    {}                    &   $g$   & $H_{c} [\rm T]$
&$m [\rm K^{-1}]$&  $U [\rm K]$ &
$\Delta [\rm K]$ & ${\bar a }$ $ [\AA]$ \\
 \hline
 Ba$_{3}$Cr$_{2}$O$_{8}$  &  1.95 &  12.10     & 0.2
&  20     &  15.85 & 3.97 \\
 Sr$_{3}$Cr$_{2}$O$_{8}$  &  1.95 &  30.40       &   0.06     & 51.2      &   39.8 & 3.82             \\
  TlCuCl$_{3}$  &  2.06 &  5.1       &   0.02   & 315
&   7.1 & 7.93             \\
 \hline
 \hline
   \end{tabular}

 \label{tab1}
\end{table}

Although $C_{0}$ in Eq.\,(\ref{Cbog}) is not good at finite temperatures, it may serve as a characteristic scale
for the contact $C$. Thus, differentiating $\Omega$ in Eq.\,\re{omegBEC} with respect to
$U$  we obtain the following explicit expression for the contact in HFB approximation
\begin{equation}
C(T<T_{c})=C_{0}\left\lbrace \widetilde{\rho} - \frac{\widetilde{\rho}_{1}^{2}+\widetilde{\rho^2}-\widetilde{\sigma}^2}{4}+2\rho_c \mu^{\prime}_{\rm eff}\left[ I_3-1-U\widetilde{\sigma}{\sigma}^{\prime}+\frac{\widetilde{\rho}_{0}}{2}+U{\rho}^{\prime}_1 (2-\widetilde{\rho}_1-\widetilde{\rho})\right] \right\rbrace
\label{CSM}
\end{equation}\normalsize
where $\widetilde{\rho}=\rho/\rho_{c}$, $\widetilde{\rho}_{1}=\rho_{1}/\rho_{c}$, $\widetilde{\rho}_{0}=\rho_{0}/\rho_{c}$, $\widetilde{\sigma}=\sigma/\rho_{c}$, ${\rho}^{\prime}_1=(\partial\rho_1/\partial\mu_{\rm eff})$, ${\sigma}^{\prime}=(\partial\sigma/\partial\mu_{\rm eff})$,
 ${\mu}^{\prime}_{\rm eff}=(\partial\mu_{\rm eff}/\partial U)=(\widetilde{\sigma}-\widetilde{\rho}_1)
 /[1+2U({\rho}^{\prime}_1-{\sigma}^{\prime})]$ and $\rho_c=\mu /2U$. The explicit expressions for $I_3$,
 ${\rho}^{\prime}_1$ and ${\sigma}^{\prime}$ are given in the Appendix. Note that, due to the factoring out
 of $C_0$, the quantity in curly brackets is dimensionless.

Now, using Eqs.\,\re{cdef} and \re{omegBEC}, yields
\begin{equation}
\begin{array}{l}
S(T<T_c)=-\ds\sum_k \ln(1-e^{-\beta E_k})+\beta \ds\sum_k E_k f_B(E_k), \\
\\
\left( \displaystyle\frac{\partial C}{\partial T}\right) _{T<T_c}=
-2m^2U^2\left(\displaystyle\frac{\partial S}{\partial U}\right)=
8\beta(Um\rho_c)^2\rho_c I_4 {\mu}^{\prime}_{\rm eff}\, , \\
\\
 \left(  \displaystyle\frac{\partial C}{\partial H}\right)  _{T<T_c}= \displaystyle\frac{1}{g\mu_B}\left( \displaystyle\frac{\partial C}{\partial\mu}\right) =-\displaystyle\frac{2m^2U^2}{g\mu_B}\left( \displaystyle\frac{\partial\rho}{\partial U}\right) =\displaystyle\frac{m^2\mu(\widetilde{\rho}-2\rho_c{\mu}^{\prime}_{eff})}{g\mu_B} = -\displaystyle\frac{2m^2U^2}{(g\mu_B)^2}
\left( \displaystyle\frac{\partial M}{\partial U}\right)\, .
\end{array}
\label{dCSM}
\end{equation}\normalsize
Unlike in the normal phase, for $T<T_c$ we failed to find a simple relation
between susceptibility and $({\partial C}/{\partial\mu})$ similar to Eq.\,(\ref{chi1}).

In the previous subsection we have shown that at zero temperature the contact may be conveniently approximated
as $C(T=0)\approx C_0= m^2\mu^2$. We now analyze the temperature dependence of $C(T<T_c)$ given in
Eq.\,(\ref{CSM}) which requires solving the set of Eqs.\,(\ref{MainEqs}).
In Fig.\,2 {a)}, we present $C(T)$ for three compounds Ba$_3$Cr$_2$O$_8$, Sr$_{3}$Cr$_{2}$O$_{8}$, and
TlCuCl$_{3}$ calculated in the HFB approximation (solid curves) using Eqs.\,(\ref{Cbig}) and
(\ref{CSM}).
It is seen that contact at finite temperature can be at least qualitatively approximated as
\begin{equation}
C(T)\approx 2U^2m^2M^2/(g\mu_B)^2\, ,
\label{C0}
\end{equation}\normalsize (dashed lines) which coincides with the exact one in the normal phase.
In {Fig.\,2 b)}
the magnetizations for the same compounds as a function of temperature are also presented for completeness.

\begin{figure}[h!]
\begin{minipage}[H]{0.49\linewidth}
\center{\includegraphics[width=1.25\linewidth]{Fig2a.pdf} \\ a)}
\end{minipage}
\hfill
\begin{minipage}[H]{0.49\linewidth}
\center{\includegraphics[width=1.25\linewidth]{Fig2b.pdf} \\ b)}
\end{minipage}
\caption{(a): Contact for spin gapped magnets in dimensionless units. The solid lines are for
HFB approximation. The dashed lines correspond to the classic approximation given in Eq.\,\re{C0};
(b): the magnetizations $M(T)$ per dimer.}
\label{Fig2}
\end{figure}

\begin{figure}[h!]
\begin{minipage}[H]{0.49\linewidth}
\center{\includegraphics[width=1.25\linewidth]{Fig3a.pdf} \\ a)}
\end{minipage}
\hfill
\begin{minipage}[H]{0.49\linewidth}
\center{\includegraphics[width=1.25\linewidth]{Fig3b.pdf} \\ b)}
\end{minipage}
\caption{Derivatives of contact with respect to magnetic field (a) and to temperature
	for TlCuCl$_{3}$. Here we use $V=1$ units, so $(\partial C/ \partial H)$ and $(\partial C/ \partial T)$
	are given in T$^{-1}$ and K$^{-1}$, respectively. The curves are obtained from Eqs.\,\re{dCSM}.
	}
\label{Fig3}
\end{figure}

\subsection{Contact and its derivatives near critical temperature}

The crude approximation in Eq.\,(\ref{C0}) is {also useful to explain} the temperature dependence of $(\partial C/\partial H)$ and $(\partial C/\partial T)$ presented in Fig.\,3a and Fig.\,3b, respectively. In fact, it follows that
\begin{eqnarray}
\begin{aligned}
\left( \frac{\partial C}{\partial H}\right )\propto
 \left(\frac{\partial M^2}{\partial H}\right )\propto M\chi\, , \\
\left(\frac{\partial C}{\partial T}\right )\propto \left(\frac{\partial M^2}{\partial T}\right )=2 M\left(\frac{\partial M}{\partial T}\right )\, .
\label{propto}
\end{aligned}
\end{eqnarray}
The experimental susceptibilities of quantum magnets with a spin gap show a well-pronounced maximum
at $T=T_c$,  and they decrease exponentially at low temperatures indicating the existence of a gap $\Delta$
between the ground-state and first excited triplet states \cite{delamorre,shermansound,tanaka2020}, while the magnetization
has a local minimum i.e., $(\partial M/\partial T) |_{T=T_c}=0$ at the critical point (see Fig.\,2b).
This is in good agreement with our predicted temperature dependence of derivatives of the contact:
$(\partial C/\partial H)$ has a maximum and $(\partial C/\partial T)$ changes its sign near the phase transition from BEC to a normal phase.

On the other hand, it is well established that the heat capacity $C_V$ has a cusp near the transition point
\cite{ourImanPRE,ourImanIJP,ourAniz2part2,Huangbook}. Therefore, we ask whether the quantities under
consideration have a jump near $T_c$.
To study this question we consider Eqs.\,(\ref{Cbig}) and (\ref{CSM}). From Eq. (\ref{Cbig}) one may see
that $C(T\rightarrow T_c^+)=2U^2m^2\rho_{c}^{2}$ . On the other side,  setting in
(\ref{CSM}) $\tilde{\rho}_1=\tilde{\rho}=1$, $\sigma=0$, $\mu_{\rm eff}=0$ and $E_k=\varepsilon _k$,
one obtains $C(T\rightarrow T_c^-)=m^2\mu^2[1/2+(I_3-1)2\rho_c{\mu}^{\prime}_{\rm eff}]=2U^2m^2\rho_{c}^{2}$,
where we used $I_3(T=T_c)=\sum_k f_B(\varepsilon _k)=\rho_c$  and  $\mu=2U\rho_c$. Therefore, Tan's contact is continuous at the critical temperature: $C(T\rightarrow T_c^-)=C(T\rightarrow T_c^+)=2(Um\rho_c)^2$.
As to its derivatives, using explicit expression for $(\partial C/\partial \mu)$
and $(\partial C/\partial T)$ given by Eqs.\,(\ref{difCTCmu}) and (\ref{dCSM}), it can be demonstrated that
they have a cusp at $T=T_c$. Presenting the details in the Appendix, we bring here the final expression for the discontinuities  in dimensionless units
\begin{equation}
\begin{array}{l}
R_{H}\equiv\frac{H_c}{C_{0}}\left\lbrace \left( \frac{\partial C}{\partial H}\right) _{T=T_c-0}-
\left( \frac{\partial C}{\partial H}\right) _{T=T_c+0} \right\rbrace=\frac{4T_c}{U}\left\lbrace \sum_k\left[1+\frac{4f_B(\varepsilon _k)T_c}{\varepsilon _k}-4f_{B}^{2}(\varepsilon _k) \right] \right\rbrace^{-1}\, , \\
\\
R_T=\frac{T_c}{C_{0}}\left\lbrace \left( \frac{\partial C}{\partial T}\right) _{T=T_c-0}-
\left( \frac{\partial C}{\partial T}\right) _{T=T_c+0}\right\rbrace=-\displaystyle\frac{R_H}{\rho_c T_c}\sum_k\varepsilon _k e^{\beta\varepsilon _k}f_B(\varepsilon _k)\, .
\label{RHURT}
\end{array}
\end{equation}\normalsize
The numerical values of discontinuities $R_H$ and $R_T$ for TlCuCl$_3$ for some values of the magnetic field are presented in Table 2.
\begin{table}[h!
\caption {Some critical parameters of TlCuCl$_{3}$.}
\begin{tabular}{p{0cm} c c c c c c }
 \hline
 \hline
    & $\ \ \ \ \ H [\rm T]\ \ \ \ $  & $\ \ \ \ \ T_{c} [\rm K] \ \ \ \ \ $ & $\ \ \ \ \ \rho_{c}\ \ \ \ \ $ & $\ \ \ \ \ C (T=T_c)[\rm nm^{-4}] \ \ \ \ \ $ & $\ \ \ \ \ R_{T}\ \ \ \ \ $ & $\ \ \ \ \ R_{\mu}\ \ \ \ \ $ \\ [0.5ex]
 \hline
  & 6.0 & 2.15 & 0.00120 & 0.000384 &  -0.0245 & 0.0134 \\
  & 6.5 & 3.06 & 0.00230 & 0.00140 & -0.0365 & 0.0190 \\
  & 7.0 & 3.73 & 0.00340 & 0.00307 & -0.0458 & 0.0230 \\
  & 7.5 & 4.28 & 0.00450 & 0.00537 & -0.0540 & 0.0263 \\
  & 8.0 & 4.75 & 0.00559 & 0.00831 & -0.0609 & 0.0290\\
  & 8.5 & 5.16 & 0.00669 & 0.0119 & -0.0671 & 0.0314 \\
  & 9.0 & 5.54 & 0.00779 & 0.0161 & -0.0727 & 0.0336 \\
  & 9.5 & 5.88 & 0.00889 & 0.0210 & -0.0779 & 0.0355 \\
  & 10.0 & 6.20 & 0.00999 & 0.0265 & -0.0826 & 0.0373 \\
  & 10.5 & 6.49 & 0.0111 & 0.0326 & -0.0870 & 0.0390 \\
  & 11.0 & 6.77 & 0.0122 & 0.0394 & -0.0911 & 0.0405 \\
  & 11.5 & 7.04 & 0.0133 & 0.0468 & -0.0950 & 0.0420 \\
  & 12.0 & 7.28 & 0.0144 & 0.0549 & -0.0986 & 0.0420 \\
 \hline
 \hline
 \end{tabular}
    \label{tab2}
\end{table}
From Table 2 it may be concluded that the jump in $(\partial C/\partial T)\propto R_T$
is negative ($R_T<0$), in contrast to the jump in $(\partial C/\partial\mu)\propto R_\mu>0$,
and in the heat capacity as well as in the Gr{\"u}neisen parameter \cite{ourMCE}.
Thus, from the relation $R_T\approx-2m^2U^2(\partial S/\partial U)\vert_{T=T_c}$  one may conclude that,
in $\mu VT$ ensembles entropy increases with increasing repulsive interaction strength,
i.e., $(\partial S/\partial U)>0$. Next, we discuss another critical phenomenon, which is
related to the low temperature critical behavior.

\subsection{Low temperature expansion and QPT}

In a $\mu VT$ ensemble, studied in the present work,  a quantum phase transition occurs at zero temperature
upon tuning the external magnetic field or equivalently the chemical potential. At $T=0$ the distance to the
quantum critical point  is determined by a control parameter, $r(H)$. Near the QCP the control parameter
can be linearized around the critical magnetic field $H_c$ as $r(H)=(H-H_c)/H_c$, which corresponds to
$\mu=r\Delta$ (where $\Delta=g\mu_B H_c$ is the spin gap).
It is expected that some thermodynamic quantities will diverge at QPT. For example, the Gr{\"u}neisen parameter
which plays an important role in magnetocaloric effect diverges as $\left.\Gamma_{H}\right|_{r\rightarrow 0}\sim 1/r$
\cite{garst,Zhu}. In the following we address the critical behavior of Tan's contact and its derivatives close to QCP.

The case of $C$ and $\left (\partial C/\partial\mu\right )$ is rather straightforward.
Actually, in Section\,II, we have shown that at $T=0$ the contact can be simply approximated as
$\left. C \right|_{T=0}=m^{2}\mu^{2}$, especially at weak magnetic fields, i.e., as $r\rightarrow 0$ (see Fig.\,1).
Thus, one directly obtains
\begin{equation}
\begin{array}{l}
 C(T=0) \vert_{ \ r\rightarrow 0}=m^{2}r^{2}\Delta^{2}+O(r^{5/2})\, , \\
\left(  \displaystyle\frac{\partial C}{\partial\mu} \right) \vert_{T=0, \ r\rightarrow 0}=2 m^{2}r\Delta+O(r^{3/2})\, ,
\label{QPT1}
\end{array}
\end{equation}\normalsize
to show that in contrast to  $\Gamma_{H}$, contact and its derivative  $(\partial C/ \partial H)$ are regular near QCP.

The case of $(\partial C/ \partial T)$ is a little more complicated. Here, one has to find the low temperature expansion
for $(\partial S/ \partial U)$, as well as small $r$ expansion for the effective chemical potential $\mu_{\rm eff}(r)$.
This leads to following expressions \cite{ourMCE}
\begin{equation}
\begin{array}{l}
\left( \displaystyle\frac{\partial S}{\partial U}\right) =-\mu_{\rm eff}^{\prime} \beta^{2}\sum_k \varepsilon_k e^{\beta E_k}f_{B}^2(E_k)=-\displaystyle\frac{\mu_{\rm eff}^{\prime}\pi^2 T^{3}}{15 m s_{0}^5}+O((Tm)^5)\, , \\
\\
\mu^{\prime}_{\rm eff}=\left( \displaystyle\frac{\partial\mu}{\partial U}\right) =-\displaystyle\frac{rQ_{0}m\Delta_{1}}{\pi\Delta}+O(r^{3/2})\, , \\
\mu_{\rm eff}=r\Delta_1+r^{3/2}\Delta_{32}\, ,
\label{Qpir}
\end{array}
\end{equation}\normalsize
where $s_0=\sqrt{\mu_{\rm eff}/m}$ is the sound velocity at $T=0$, $\Delta_1=\Delta\pi/(\pi+Q_{0}Um)$, $\Delta_{32}=4(\Delta m)^{3/2}\sqrt{\pi}U/3(Q_{0}Um+\pi)^{5/2}$ and $Q_{0}=(6/\pi)^{1/3}$ is the Debye momentum
\cite{ourMCE}.
Then, the low temperature expansion for $(\partial C/\partial T)$ near QCP can be directly obtained as
\begin{equation}
\left( \frac{\partial C}{\partial T}\right) =-2m^{2}U^{2}\left(\frac{\partial S}{\partial U} \right)=
-\alpha\frac{T^{3}}{r^{3/2}}+O((Tm)^{5})\, ,
\label{dCdTQCP}
\end{equation}\normalsize
where $\alpha={m^{9/2}U^{2}\pi Q_{0}}/{15\Delta\sqrt{\Delta_1}}$.
It is seen that at low but finite temperatures $(\partial C/\partial T)$ diverges as $1/r^{3/2}$, but there
is no divergence at exactly $T=0$ , where QPT occurs.
Numerical estimations show that $\alpha$ is positive and rather small and hence, at low temperatures the contact decreases very slowly as seen in {Fig.\,2 a).}

\section{Discussion and Conclusion}

We have studied the Tan's contact for spin-gapped quantum magnets, assuming that their low temperature properties are related to that of a triplon gas.  To the best of our knowledge, this is the first time the parameter $C$ and its dependence on temperature is investigated for a solid.
Naturally,  it will be interesting to compare our results with existing ones in the literature.
However, there is neither a theoretical nor an experimental study of Tan's contact for $\mu VT$ ensembles at finite temperatures and most of the literature on $C(T)$ concern 1D gases in $NVT$ ensembles.

The temperature dependence of $C$ of harmonically trapped 1D Leib-Liniger Bose gas has been studied by Yao \textit{et al.} \cite{YaoPRE121}.
They have found that contact increases at low temperatures, reaches a maximum at $T=T^{*}$ and then starts to
decrease, i.e., $\left.(\partial C/\partial T)\right|_{T<T^{*}}>0$, $\left.(\partial C/\partial T)\right|_{T=T^{*}}=0$ and
$\left.(\partial C/\partial T)\right|_{T>T^{*}}<0$.
From Figs.\,1 and 2, one may note that for a 3D system of bosons the situation is quite the opposite;
$\left.(\partial C/\partial T)\right|_{T<T_c}<0$ and $\left.(\partial C/\partial T)\right|_{T>T_c}>0$.
But at the critical temperature $T_c$, in both cases $C(T)$ exhibits an extremum as a function of temperature,
maximum for 1D bosons and minimum for triplons, i.e., in both cases $(\partial C/\partial T)$, changes its sign
near the critical temperature.
Bearing in mind the relation $(\partial C/\partial T)\sim (\partial S/\partial U)$, one may conclude that at $T=T_c$ the interaction dependence of the entropy displays a maximum.
Note that in the case of 1D bosons this maximum provides a signature of the crossover to the fermionized regime, while in the case of quantum magnets it corresponds just to the point of finite temperature phase transition, where accumulation of entropy occurs \cite{garst}.

In the present work the temperature dependence of $C(T)$ at large temperatures $T>T_c$, is found to be
almost linear. This is to be compared with the results by Vignolo and Minguzzi \cite{Vignolo}
who obtained the large temperature behavior as $C\sim\sqrt{T}$ for a 1D Bose gas in the Tonks-Girardeau limit.
It would be interesting to study the effect of dimensionality on $C(T)$ in $\mu VT$ ensembles  also
\cite{kosterlic}.

We have  studied the critical behavior of $C$ at the quantum phase transition and found that $C$ and
$(\partial C/\partial\mu)$ are regular at QCP which is in good agreement with predictions made by
Chen \textit{et al.} \cite{chencrit}.

In conclusion it should be underlined that our analytical expressions for Tan's contact of spin gapped
compounds are expressed in terms of magnetizations which are directly observable.
For example, at very low temperatures, $T\ll T_c$, one evaluates $C$ simply from the expression
$C=m^2\mu^2/{\bar a}^{4}$, where the only parameter, effective mass should be estimated.
Note that, in this region contact of a $\mu VT$ ensemble is almost not sensitive to the strength of
the boson-boson interaction $U$, while the dependence on $U$ of the contact of an $NVT$ ensemble
is rather strong even in the classical approximation, when $C$ is  given by
$C(NVT)\approx 16\pi^2a^2\rho^2=U^2 m^2 \rho^2$ \cite{ourTan1,Pitbook14}.

Finally, we venture to suggest a way to measure $C$ in quantum magnets,
{which can be performed in a similar way as it was done in Ref.\,[\onlinecite{wild2012}] (see Appendix B for
details).}

We hope our work will stimulate more studies, especially experimental ones exploring the universality of
$\mu VT$ ensembles in quantum magnets. {Results presented here will deepen our understanding of the
connection between short-ranged two-body correlations and magnetic phase transitions in
antiferromagnetic materials, which
have the potential to be used in developing the next generation of spintronic devices.}

\acknowledgements
This work is supported by the Ministry of Innovative Development of the Republic of Uzbekistan and
the Scientific and Technological Research Council of Turkey (TUBITAK) under Grant No.\,119N689. AR
acknowledges partial support from the Academy of Sciences of the Republic of Uzbekistan. BT acknowledges
support from the Turkish Academy of Sciences (TUBA).

\section*{Appendix A}
\defB.\arabic{equation}{A.\arabic{equation}}
\setcounter{equation}{0}

In this appendix we provide explicit expressions for the integrals introduced in
Eqs.\,(\ref{drodSUBg}), (\ref{CSM}) and (\ref{dCSM}):
\begin{eqnarray}
\begin{aligned}
 &I_1=-\beta\sum_k\omega_k f_{B}^{2}(\omega_k)e^{\beta\omega_k}, \ \ \ \ \ \
I_2=-U\beta\sum_k f_{B}^{2}(\omega_k)e^{\beta\omega_k}, \\
 &I_3=\rho_{c}^{-1}\sum_k\left(\frac{\varepsilon_k W_k}{E_k}-\frac{1}{2} \right), \ \ \ \ \ \
I_4=\beta\sum_k f_B^{2}(\omega_k)\varepsilon_k e^{\beta\omega_k}.
\label{I1234}
\end{aligned}
\end{eqnarray}
where $f_{B}(x)=1/(e^{\beta x}-1)$, $W_k=(1/2+f_B(E_k))$, $\omega_k=\varepsilon_k-\mu+2U\rho$ and $E_{k}^{2}=\varepsilon_k(\varepsilon_k+2\mu_{eff})$.

Equation\,(\ref{CSM}) includes $(\partial\rho_{1}/\partial\mu_{eff})$,  $(\partial\sigma/\partial\mu_{eff})$ and $(\partial\mu_{eff}/\partial U)$. The derivatives with respect to $\mu_{eff}$ can be found by differentiating
Eqs.\,(\ref{MainEqs}) to give
\begin{eqnarray}
\left( \frac{\partial\rho_1}{\partial\mu_{eff}}\right) =\sum_k\frac{\varepsilon_k}{E_{k}^2}\left[  \frac{\mu_{eff}W_k}{E_k}+\frac{(\varepsilon_k+\mu_{eff})W_{k}^{\prime}}{4}\right]\, , \\
\left( \frac{\partial\sigma}{\partial\mu_{eff}}\right) =-\sum_k\frac{\varepsilon_k}{E_{k}^2}\left[ \frac{\mu_{eff}W_{k}^{\prime}}{4}+\frac{(\varepsilon_k+\mu_{eff})W_{k}}{E_k}\right]\, ,
\label{drho1dsig}
\end{eqnarray}
where $W_{k}^{\prime}=\beta(1-4W_{k}^{2})$.

As to $(\partial\mu_{\rm eff}/\partial U)$, it can be obtained by differentiation of both sides of the equation
$\mu_{\rm eff}(U)=\mu+2U\left[ \sigma(\mu_{eff}(U))-\rho_1(\mu_{eff}(U))\right]$ with respect to $U$
and then solving it for $(\partial\mu_{eff}/\partial U)$. This leads to
\begin{eqnarray}
\left( \frac{\partial\mu_{\rm eff}}{\partial U}\right) =\frac{\mu_{eff}-\mu}{U D}\, , \\
D=1+2U\left[\left( \frac{\partial\rho_1}{\partial\mu_{eff}}\right) -\left( \frac{\partial\sigma}{\partial\mu_{eff}}
\right)  \right]\, .
\label{dmueffdU}
\end{eqnarray}
Now, we proceed to evaluate $(\partial C/\partial H)=-2m^2U^2g\mu_B(\partial\rho/\partial U)$ and
$(\partial C/\partial T)=-2m^2U^2(\partial S/\partial U)$ near the critical temperature $T_c$.

From the explicit expression
\begin{equation}
S(T>T_c)=\sum_k \ln(1-e^{-\beta\omega_k})+\beta\sum_k\omega_k f_B(\omega_k),
\label{STbig}
\end{equation}\normalsize
where $\omega_k=\varepsilon_k-\mu+2U\rho(U)$, one obtains
\begin{equation}
\left( \frac{\partial S}{\partial U}\right)_{T=T_c+0}=-\beta^{2}\sum_k e^{\beta\omega_k}\omega_k\omega^{\prime}_{k}f_B^2(\omega_k)\, ,
\label{dSdUTbig}
\end{equation}\normalsize
where $\omega^{\prime}_k=(\partial\omega_k/\partial U)=2\rho+2U(\partial\rho/\partial U)$ and
$(\partial\rho/\partial U)$ is given by Eq.\,(\ref{drodSUBg}). It is straightforward to show that the integral $I_2$ in
Eq.\,(\ref{drodSUBg}) is divergent at $T=T_c$, where $\omega_k=\varepsilon _k\sim k^2/2m$, and hence
$\left.(\partial\rho/\partial U)\right|_{T=T_c}=-\rho/U$. This means that $\omega^{\prime}_k(T=T_c^+)=0$ and
$\left.(\partial S/\partial U)\right|_{T=T_c^+}=0$. Therefore $\left.(\partial C/\partial T)\right|_{T=T_c^+}=0$ and
$\left.(\partial\rho/\partial U)\right|_{T=T_c^+}=-\rho/U$.

Replacing $\omega_k$ by $E_k$ in Eq.\,(\ref{STbig}) and taking the derivative one finds
\begin{equation}
\left( \frac{\partial S}{\partial U}\right)_{T<T_c}=-\beta^2\left(\frac{\partial\mu_{\rm eff}}{\partial U}\right)\sum_k\varepsilon _k e^{\beta E_k}f^2_B(E_k).
\end{equation}\normalsize
Now, taking into account that $\mu_{\rm eff}(T_c)=0$,
$\rho(T_c)=\rho_c=\mu/2U$, $E_k(T_c)=\varepsilon _k$ and using Eq.\,\re{dmueffdU} we obtain
\begin{equation}
\left( \frac{\partial S}{\partial U}\right)_{T=T_c^-}=\frac{\mu\beta^2}{UD}\sum_k\varepsilon _k e^{\beta\varepsilon _k}f^2_B(\varepsilon _k)\,
\end{equation}\normalsize
which proves the Eq.\,(\ref{RHURT}). The derivative $\left.(\partial\rho/\partial U)\right|_{T=T_c^-}$ may be obtained from
$ {\rho(U)=(\mu_{\rm eff}+U)/2U} $ as
\begin{equation}
\left( \frac{\partial\rho}{\partial U}\right)_{T=T_c^-}=-\frac{\mu(1+D)}{2U^2D}\, ,
\end{equation}\normalsize
where $D$ is given by Eq.\,(\ref{dmueffdU}) and $(\partial \mu_{eff})/(\partial U)$ should be evaluated at $T=T_c$.

The expansion of entropy and the effective chemical potential near QCP may be found in the appendix
of Ref.\,[\onlinecite{ourMCE}].

\section*{Appendix B}
\defB.\arabic{equation}{B.\arabic{equation}}
\setcounter{equation}{0}

{Here we outline a sketch of a possible method of  measuring Tan's contact
for spin gapped magnets. We presume that this can be achieved in a similar way to that for atomic gases.}

{Wild et al. \cite{wild2012} measured $C$ for $^{85}$Rb atoms in a gas with a
$60 \%$ condensate fraction. The atoms are in the $|F=2,m_F=-2>$ state, where $F$ is the total atomic
spin and $m_F$ is the spin projection. Then by applying an rf pulse the atoms are driven from
$|2,-2>$ to $|2,-1>$ state. The rate for transferred atoms to the final spin state is given by \cite{Br9,wild2012}
\begin{equation}
\lim_{\omega\rightarrow\infty}\Gamma(\omega)= \frac{\Omega^2}{4\pi}F(a,\omega,m)C\, ,
\lab{Gamma}
\end{equation}
where  $\Omega$ is the Rabi frequency, $F(a,\omega,m)$ is a
given function of scattering length $a$, the frequency $\omega$, and  the  atomic  mass $m$. Tan's contact has been
extracted from this equation after a direct measurement of $\Gamma$ as a function of $\omega$.}

\begin{figure}[!t]
\begin{center}
\includegraphics[scale=0.4]{fig4.pdf}
\end{center}
\caption{Cartoon of the spin levels of a typical spin gapped quantum magnet showing that the
there are $S_z=\pm 1$ excited doublet and $S_z=0$ ground states.
They are separated by a gap $\Delta=g\mu_B H_c$, which can be collapsed
due to Zeeman effect for $H>H_c$.
}
\end{figure}

{
Now we consider the spin gapped magnets. Here for $H>H_c$ we have two spin levels with
$S_z=1$ (ground state) and $S_z=-1$ (nearest excited state), as illustrated in Fig.\,4. The direct transition
between the spin singlet and triplet states is, in principle, forbidden in magnetic dipole transitions.
Nonetheless, such transitions have been observed in many spin gapped systems by means of high frequency electron
spin resonance (ESR) measurements \cite{kimura2004,kimura2018,kimura2020,amr}.
It was shown that, these transitions are, in fact,  driven
by ac electric fields, and observed electric dipole transitions can be explained by spin dependent polarization.
Very recently, Matsumoto et al. \cite{amr} developed a theoretical description of
singlet-triplet transitions in spin gapped magnet KCuCl$_3$ in the framework of spin-wave theory,
and proposed practical formulas for the transition probabilty $W$, using Fermi's Golden rule (see e.g.
Eq.\,(18) of Ref.\,[\onlinecite{amr}]). We surmise that if
one uses an effective Hamiltonian similar to Eq.\,\re{Ham} the transition probability from Ref.\,[\onlinecite{amr}]
can be related to the Tan's contact as in Eq.\,\re{Gamma}, and hence one will be able to obtain $C$ from measurements of $W$, using proper normalization of observed quantities.
}


\section{The CBM experiment and its muon detector}
\label{sec:cbm_much}

CBM is a fixed-target experiment for the study of energetic nuclear collisions at beam momenta from 3.5 -- 12 GeV per nucleon for heavy nuclei~\cite{cbmpaper}. It will be operated with external beams from the FAIR accelerator complex from 2025 on. CBM comprises two detector configurations: a setup for measurements of hadrons and electrons, with a RICH and a TRD detector for electron identification and a TOF detector for hadron identification, and a setup for measurements of muons, in which the RICH detector is replaced by the muon system MuCh. Both configurations share the central tracking system STS inside the field of a superconducting dipole magnet, and the forward calorimeter PSD. The muon configuration of CBM, with RICH in parking position, is shown in Fig~\ref{fig:muon_setup}.

\begin{figure}[htb]
\begin{center}
\includegraphics*[width=0.7\linewidth]{figures/CBMCatia.jpeg}
\caption{Layout of the CBM experiment in its muon configuration. RICH will be in parking position during MuCh dating data. The beam enters from the left. The magnet houses the target, MVD and the central tracking system STS. The muon system is placed directly after the magnet, outside of the magnetic field. Further downstream are the TRD, TOF and PSD detectors. }
\label{fig:muon_setup}
\end{center}
\end{figure}

The MuCh detector~\cite{muchwiki, muchtdr} serves for the identification of muons by filtering out other particle species in massive absorbers. The absorbers are arranged as five segments with three detector layers in between each pair of segments. The segmented absorber system allows to trace particle trajectories through the entire setup, which would not be possible by a monolithic absorber because of small-angle scattering of the muons in the absorber material. A thick absorber also likely to absorb muons without generating signal. Figure~\ref{fig:much_geometry} (left) shows the MuCh geometry as implemented in cbmroot.

The hit rates in the detector layers differ significantly because of the successive absorption of particles by the absorber segments. Thus, different detector technologies are employed: the two upstream stations (detector triplets) will use Gas Electron Multiplier (GEM detector), while the two downstream stations will be constructed from high-rate Resistive Plate Chambers (RPC). In both cases, the detector readout planes are segmented into pads in a $r-\phi$ geometry, the granularity depending on the layer position and the radial distance from the beam. The readout segmentation was optimised based on efficiency and signal-to-background ratio for the detection of muon pairs in heavy-ion collisions at FAIR energies~\cite{muchtdr, segmentation}. The pad layout is illustrated in Fig.~\ref{fig:much_geometry} (right) on the example of the first layer of the first station.

\begin{figure}[htb]
\begin{center}
\includegraphics*[height=0.45\linewidth]{figures/MUCHSetup.png}
\hfill
\includegraphics*[height=0.45\linewidth]{figures/pad_layout}
\caption{Left: Schematic layout of the MuCh detector, consisting of a set of absorber segments inter-laid with detector stations. Right: Pad segmentation of the first layer of the first station. $1^{\circ}$ in  $\phi$, \& $\delta$r $=$ $r\delta\phi$.}
\label{fig:much_geometry}
\end{center}
\end{figure}

\section{Detector response simulation}
\label{sec:digitizer}

\subsection{Simulations in cbmroot}

In cbmroot, simulations are performed in two separate steps with intermediate file output. In the first step called ``transport simulation‘‘, the primary particles are traced through the detector, taking into account their trajectory in the magnetic field, their interaction with the materials and - for unstable particles - their decay. Secondary particles created by decay or interaction with materials are created and traced as well. This step employs external engines; the framework makes use of the ROOT TVirtualMC features which allow to choose between GEANT3 and GEANT4 on the discretion of the user~\cite{vmc}, however, GEANT3 is used for all the results mentioned in this article. At transport level, a realistic description of the detector geometries, the proper material properties and a map of the magnetic field are required~\cite{MuchGeometry}. As output of the transport simulation, the geometric intersections of particle trajectories with active detector elements (``MCPoints'') are recorded, along with the time-of-flight from the event start and the energy deposit in the active material. This information is stored to file and provides the input for the second simulation step.

The task of the detector response simulation is to model the physics processes in the active detector material, the readout and the digitization in the front-end electronics. The output objects are called ``digi'' represent the simulated measurements of a single read-out channel - corresponding to one read-out pad in the case of the MuCh. This atomic data unit contains the information on the respective channel address, the time of the measurement termed as time-stamp, and the digitized charge. The digi objects are then forwarded to the software emulation of the CBM DAQ system as implemented in cbmroot. At this step, the association of digis to events is destroyed, resulting in a data stream similar to that expected from the actual experiment. In the time-based simulation, based on the dead-time discussed later, a single pad might be recorded as two digis with different time-stamp.

In the context of a free-streaming readout as described in section~\ref{sec:intro}, the correct description of the timing behaviour is of utmost importance. The effects contributing to the finally registered time stamp are:
\begin{itemize}
    \item Event time: the time corresponding to the actual collision. It is generated by the cbmroot framework from the time profile of the beam and the interaction probability of the beam particles in the target.
    \item Time-of-flight: from the event start to the intersection with the detector element. This time is provided by the transport engine and stored in the MCPoint object.
    \item Drift time: the time taken by the primary electrons in the GEM or RPC detectors to drift to the readout planes.
    \item Time response of the readout ASIC to the analog charge collected in the readout pads.
\end{itemize}

An important figure is the dead time of the readout electronics, i.e., the time after a hit in which the corresponding readout channel is blocked, such that a second hit arriving within this time is neglected. This leads to a detector inefficiency and track pile-up which are obviously dependent on the interaction rate, i.e., the time separation between two subsequent collisions.

\subsection{Analogue response simulation}

The working principle of a multi-layer GEM detector is illustrated in Fig.~\ref{fig:analog_response}. A traversing charged particle creates a number of primary electrons in the drift gap through ionisation of the gas. The primary electrons drift in the electrical field towards the GEM foils, where amplification through the creation of avalanches takes place. This process is repeated in various transfer gaps between successive GEM foils. Finally, the produced electrons are further amplified in the induction gap and collected on the pads of the readout PCB.

Our analogue simulation implements a simplified scheme as shown in the right panel of Fig.~\ref{fig:analog_response}. The trajectory of the charged particle in the drift gap is approximated by a straight line, the coordinates of the entry and exit points being provided from the transport simulation stage. For each particle, the number of primary electrons is sampled from a Landau distribution, the parameters of which depend on particle type, energy, track length in the active volume and specifications of the gas mixture~\cite{gemparameters}. The primary electrons are randomly generated along the trajectory in the active volume with a uniform distribution.

For each primary electron arriving at the amplification gap, the number of secondary electrons is calculated from the gas gain settings, and the charge depositions to the gaps are calculated assuming a spot radius of the avalanche. Depending on the track inclination, the gas gain and the avalanche profile, a particle can activate one or several pads. The drift time is the time taken by the primary electron to traverse the active volume (up to the first GEM foil). The avalanche travels with close to the speed of light, such that its propagation time is negligible in comparison to the drift time. Primary electrons generated in the passive volume do not contribute measurably to the total signal since no avalanche production takes place.

\begin{figure}[htb]
\begin{center}
\includegraphics*[height=0.3\linewidth]{figures/gem_response_real}
\hfill
\includegraphics*[height=0.3\linewidth]{figures/gem_response_simple}
\caption{Left: working principle of a multi-layer GEM detector; Right: simplified scheme as implemented in the analogue response simulation}
\label{fig:analog_response}
\end{center}
\end{figure}

The parameter input to this analogue simulation consists of the parameters of the Landau distribution for the number of primary electrons, the drift velocity of the primary electrons, the gas gain, and the avalanche spot radius. All of these parameters are tunable in the software. The size of the drift gap is taken as per geometry of GEM or RPC. The analogue simulation process is the same for GEM and RPC detectors, but with different parameters. The default settings for both detector types are summarised in Table~\ref{tab:analog_parameters}.

\begin{table}[h]
    \begin{center}
        \caption{\label{tab:analog_parameters}Parameters for the analogue simulation of GEM and RPC detectors}
        \begin{tabular}{l|l|l}
        \hline
        Parameter & GEM & RPC \\
        \hline
        MPV for primary Landau distribution & from HEED & 12 \\
        Drift gap (fixed as per detector geometry)  & 3 mm & 2 mm\\
         Drift Velocity of primary electrons ($\mu$m/ns)& 100 & 120\\
        Gas gain & 5000 & 30000\\
        Spot Radius & 500 $\mu$m & 2 mm\\
        \hline
        \end{tabular}
    \end{center}
\end{table}

The analogue simulation as described above is performed separately for each input MCPoint.
The resulting charge depositions to the readout pads, however, cannot directly be digitised, since they block the respective electronics channel for a certain amount of time (dead time), thus potentially influencing later measurements. Such incidences of different particles contributing charge to the same pad can happen within one event (conventional pile-up), but in our free-running scenario also for particles from different events, provided the difference in the event times is comparable or smaller than the dead time. To cope with this situation, the registered charges per pad are internally buffered and keeps a match object associating a link per corresponding incident track. Later this match information is used for finding the primary and secondary track contribution towards generation of digi or pile-up. The stored information also contains the pad address, the signal time, the charge, and the duration of the active signal.

\subsection{Digital response simulation}

The digital response to the accumulated charge in the readout pads is modelled following the properties of the readout ASIC, the STS/MUCH-XYTER (SMX)~\cite{iaea,nxyter}. This ASIC features two shaping channels, the fast shaper used to determine the time stamp of the measurement, and the slow shaper which measures the charge amplitude by time-over-threshold. Once a signal in the fast shaper crosses a pre-defined threshold, the signal in the slow shaper is evaluated and a message to the DAQ system is issued when the amplitude in the slow shaper falls below a second threshold. The signal shape in the slow shaper thus defines under which conditions two subsequent signals in the same readout channel can be disentangled.

This situation is modelled in our software by a delta function for the signal shape, the start time being the time stamp and the duration being the dead time of the ASIC. The amplitude is the total charge. Two overlapping signals (i.e., the second one arriving within the dead time of the first) are merged into one, having the start time of the first signal, the stop time of the second, and the amplitude the sum of both signals. The generated signals are transiently stored in a readout buffer. When the system time (the time of the currently processed event) exceeds the stop time of the signal, it is digitised, since then it is guaranteed that it will not be influenced by subsequently arriving signals.

All analogue signals exceeding the threshold are digitised, using the parameters of the ADC built into the SMX. The charge is linearly discretized, above a tunable threshold, into 32 bins corresponding to the 5 bit resolution. The time of the digital signal is smeared by a Gaussian distribution representing the time resolution of the ASIC. The parameter input to the digital response simulation are thus the number of ADC channels, the ADC threshold and dynamic range, and the time resolution. The currently used values of these parameters are summarised in Table~\ref{tab:digital_parameters}. It should be noted that these parameters are still subject to changes; in particular, the dead time, which is a critical parameter for the performance of the detector system, is expected to be smaller for SMX v2.2 currently under development~\cite{SMX}.

\begin{table}[h]
    \begin{center}
        \caption{\label{tab:digital_parameters}Parameters for the digital response simulation of GEM and RPC detectors}
        \begin{tabular}{l|l|l}
        \hline
        Parameter & GEM & RPC \\
        \hline
        Number of ADC channels& 32 (5 bit) & 32 (5 bit)\\
        ADC threshold    & 2 fC & 30 fC \\
        ADC dynamic range  & 80 fC & 130 fC \\
        Dead time   &  400 ns   &   400 ns  \\
        Time resolution  &  5 ns  &  5 ns \\
        \hline
        \end{tabular}
    \end{center}
\end{table}

\subsection{Noise simulation}

Our simulation also covers thermal noise generating random signals from the readout electronics. Thermal noise charge is always present in the readout channel; the distribution of its value is approximately Gaussian. Whenever the noise charge exceeds the threshold, a message is generated just as in case of charge originating from a traversing particle. The occurrence of this noise is random and not correlated in time to events.

In our implementation, the noise is described by a tunable parameter, the noise rate per pad. From this parameter, the number of noise signals is sampled for each module for the time period from the previous to the current event time from a Poissonian distribution with the expectation value noise rate $\times$ number of pads per module x time interval. The time of the noise signals is sampled from a uniform distribution in the respective time interval; their charge is sampled from the Gaussian noise distribution.

Noise signals are inserted into the readout buffer and thus into the simulated data stream. They can interfere with ''real`` signals from traversing particles in the same way as two ''real``signals (see the previous subsection).

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

The data acquisition systems of most experiments in high-energy particle or nuclear physics are based on a hardware trigger, where a signal generated by a suitable set of hardware indicates that a collision took place and triggers the timely readout of the front-end electronics. The hardware trigger thus defines ``event'' as data collections of moderate size representing a separate single collision, which may then be either written to permanent storage or subjected to further inspection and selection by higher-level triggers, e.g., on FPGA or in software. The software framework used for simulation and analysis of such experiments are thus designed on an event-by-event scheme, where each event is treated as a separate and independent entity.

Several next-generation experiments, however, face the difficulty that such a triggered readout scheme is not feasible, e.g., because of very high interaction rates and/or complex trigger topologies which are not well suited to be evaluated in hardware logic and time consuming. An example is the Compressed Baryonic Matter experiment (CBM)~\cite{cbm} at the upcoming Facility for Antiproton and Ion Research (FAIR)~\cite{fair} in Darmstadt, Germany. CBM intends to inspect up to $10^7$ nuclear collisions per second, each producing several hundreds of particles to be registered in the detectors. These conditions will use a trigger-less, free running data acquisition, where self-triggered front-end electronics elements register signals above threshold caused by particles traversing the respective detectors and autonomously push the data forward~\cite{TriggerLess}. Such a system is not limited by latency, i.e. the time needed to generate a hardware trigger, but by data throughput bandwidth. It results in a continuous data stream in contrast to a series of events defined by the hardware trigger in conventional readout schemes. The association of raw data to physical events is based on a precise time stamp provided to each measurement by a central timing system.

The software to analyse this data stream - both in real-time and offline - must cope with this situation~\cite{friese2012}. The same holds also for the software to simulate the experiment, which must generate such a data stream from appropriate physics models of high-energy nuclear collisions. The simulation thus has to convert a series of events into a stream of data incorporating response of detectors, front-end electronics and data acquisition to these events. The CBM software framework cbmroot~\cite{cbmroot}, based on the FairRoot simulation and analysis framework~\cite{FlorianFairRoot}, provides the appropriate structures for this task, in particular the software emulation of the data acquisition collecting the free-streaming data from the detector front-ends~\cite{friese2017}. These framework structures must of course be filled by the implementation of the actual detector response of the various CBM sub-systems.

In this article, we describe the time based simulation of the CBM muon detector system MuCh. In section~\ref{sec:cbm_much} the CBM experiment and its muon detector are introduced. Section~\ref{sec:digitizer} then describes the implementation of the detector response for the MuCh detector. Verification and performance of the simulation are discussed in section~\ref{sec:results}. Finally, section~\ref{sec:summary} gives a summary and conclusions.

\section{Results and detector performance}
\label{sec:results}

The result of the MuCh simulation is a stream of digi objects integrated by the DAQ software into the streams of the other CBM detector systems. Figure~\ref{fig:datastream} shows the digi rate obtained from the simulation. The MC-true event time generated by the cbmroot simulation framework is overlaid for comparison. We see that events can be identified by peaks in the distribution of digis~\cite{eventbuilding}.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.8\linewidth]{figures/digi_time_dist_minbias_10agev_au_au_noise_on_bin_size_1ns.pdf}
\caption{The time distribution of the MuCh digis for a certain time period ( within a time slice), with several events in and with moderate noise switched on, resembling self-triggered free running data stream.}
\label{fig:datastream}
\end{center}
\end{figure}

The time distribution of digis within a single event is shown in figure~\ref{fig:time_event} together with the corresponding Monte-Carlo true time of the track in the detector. Included is also the reconstructed hit time, where a hit results from a cluster-finding procedure which groups simultaneously active neighbouring pads. The MC true time shows a narrow peak and an approximately exponential tail. The peak originates from both primary and secondary tracks, which are registered within a time span of less than 10 ns. Due to interaction with absorber material of muon system, 94.5\% digis, registered in MuCh detector layers, are generated from secondary tracks. The offset of the peak (about 7 ns) reflects the time-of-flight from the main interaction vertex to the MuCh detectors. The exponential tail comprises secondary tracks only. In comparison, the digi time distribution shows a broadening of the peak due to the time resolution and an additional offset stemming from the drift time of the electrons in the drift gap of the GEM counters (see Fig.~\ref{fig:analog_response}). Note that the integrals of the distributions of MC points and digis are different since one track can activate more than one pad and more than one MC points can generates one digi. The hit time, on the other hand, is corrected for the electron drift and thus shows no offset with respect to the MC origin.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.8\linewidth]{figures/time_event}
\caption{Time distribution of Monte-Carlo points, digis and hits in a single event in the first station of the MuCh detector. Hits are reconstructed from clusters of digis in neighbouring active pads. The origin of the time axis corresponds to the Monte-Carlo event time.}
\label{fig:time_event}
\end{center}
\end{figure}

The CBM experiment intends to operate at high interaction rates of up to $10^7$ events per second. At such rates, the average time between two subsequent events (100 ns for $10^7$ events/s) becomes comparable to the dead time of the readout electronics (400 ns). Consequently, losses due to pile-up from particles originating from different events are expected. The developed simulation software as described in the previous section allows to study such timing effects and to quantify rate-dependent losses, which is essential to determine the performance of the detector.

For the following studies, we simulated both minimum-bias and central Au+Au collisions at 10 GeV/$c$. In the minimum-bias event sample, the collision impact parameter $b$ is realistically distributed according to the geometric cross section ($\propto b\mathrm{d}b$). The time sequence of such events corresponds to the actual experimental situation. Central events ($b=0$) create the highest track multiplicity and thus the highest data load on the detectors. A time series of central events only does not correspond to the physical situation, but is used here to verify the simulation software since it magnifies effects arising from the finite dead time.

A common detector parameter figure is the occupancy, conventionally defined as the probability for a single channel to be activated in one event. In our free-running scenario, we calculate this quantity as the number of digis divided by the number of pads and the number of simulated events. Figure~\ref{fig:occupancy} shows this occupancy (averaged over one detector layer) for the first layer of the first station at an event rate of $10^7$/s as function of the single-channel dead time. The decrease of occupancy with dead time, moderate for minimum-bias events and better visible for central events, is an expected consequence of the data loss due to signals arriving in the dead time of a previous signal.

\begin{figure}[htb]
\begin{center}
\includegraphics*[width=0.6\linewidth]{figures/occupancy}
\caption{Average single-channel occupancy in the first layer of the first station as function of dead time for Au+Au collisions at $10^7$ events/s}
\label{fig:occupancy}
\end{center}
\end{figure}

Rate-dependent data losses can also be qualified by studying the pile-up occurrences. In our simulation as earlier mentioned, the reference to the Monte-Carlo origin is kept for each digi object. Pile-up thus happened if a digi object has references to more than one Monte-Carlo particles. We define the pile-up fraction as the number of digis with multiple MC reference divided by the total number of digis. This fraction is shown in Fig.~\ref{fig:pileup_deadtime} again as a function of dead time for an interaction rate of $10^7$/s.  For minimum-bias events, the fraction slowly varies with the dead time; for a value of 400 ns, corresponding to the current SMX design, it amounts to about 13\%. It should be mentioned that as discussed later the maximum contribution to this pile-up fraction is from the secondary tracks generated in the material. This does not represent the pile-up for all tracks (primary \& secondary) associated with the incident primaries, which will govern the capabilities of tracking of primary tracks in the detector.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.6\linewidth]{figures/pileup_deadtime.pdf}
\caption{Average pile-up fraction in the first layer of the first station as function of dead time for Au+Au collisions at $10^7$ events/s.}
\label{fig:pileup_deadtime}
\end{center}
\end{figure}

The pile-up fraction is shown in Fig.~\ref{fig:pileup_sector} differentially for each sector in the first layer of the first station for different interaction rates, varying from $10^4$/s to $10^7$/s. Here one sector denotes one ring containing 360 pads which is according to progressive geometry (see right fig.~\ref{fig:much_geometry}). Again, the rate effect is better visible for central events (right-hand panel) than for minimum-bias events (left-hand panel). Pile-up at low rates happens between tracks within the same event and is thus irreducible by varying the dead time; it corresponds to conventional ''double hits‘‘ and is fixed by the granularity (pad layout) of the detector. From an event rate of about $10^6$/s on, the additional pile-up between tracks from different events becomes visible and increases with increasing rate.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.49\linewidth]{figures/pileup_sector_mbias.pdf}
\hfill
\includegraphics[width=0.49\linewidth]{figures/pileup_sector_central.pdf}
\caption{Pile-up fraction in each sector of the first layer of the first MuCh station for a dead time of 400 ns and for different interaction rates. The left panel shows the (realistic) case of minimum-bias events, the right panel the (artificial) case of central events only.}
\label{fig:pileup_sector}
\end{center}
\end{figure}

The situation in the other three MuCh stations is comparable to the one in the first station, as shown in Fig.~\ref{fig:pileup_station}. In all stations, an increase of inter-event pile-up with rate is observed on top of the irreducible in-event pile-up. Note that for constructional reasons, station 2 has the same angular pad segmentation as station 1, although the hit density is significantly reduced in the absorber segment in between which can be seen in left panel of fig.~\ref{fig:occupancy_sector}. This results in a lower occupancy and a lower amount of pile-up losses compared to station 1. The pad layout of stations 3 and 4 was chosen such as to result in an occupancy similar to station 1. The right panel of fig.~\ref{fig:occupancy_sector} shows the digi occupancy for first layer of $3^{rd}$ and $4^{th}$ MuCh station as function of sector.

\begin{figure}[htb]
\begin{center}
    \includegraphics[width=0.6\linewidth]{figures/pileup_station.pdf}
\caption{Average pile-up fraction as a function of interaction rate for minimum-bias Au+Au collisions at 10$A$GeV/c. The single-channel dead time was taken to be 400 ns.}
\label{fig:pileup_station}
\end{center}
\end{figure}

The consequence of the pile-up for the sensitivity of the detector to physics observables must be assessed by proper simulations including the full reconstruction of the simulated data stream in terms of hits and tracks. Such studies are yet to come. However, identifying the Monte-Carlo origin of the pile-up cases, we find that 92.47\% of all pile-ups happen between two secondary tracks in a shower created close to the downstream edge of an absorber layer. Such shower tracks are concentrated in a small emission cone are thus more likely to deliver charge into the same read-out pad. About 7.4\% of pile-ups occur between a primary and a secondary track; the fraction of pile-ups between two primary tracks (from the main event vertex) is negligible. This suggests the impact of pile-up on the physics performance of the detector to be minute.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.49\linewidth]{figures/occupancy_minbias_auau_1st_2nd_station.pdf}
\hfill
\includegraphics[width=0.47\linewidth]{figures/occupancy_minbias_auau_3rd_4th_station.pdf}
\caption{Digi occupancy as function of sector of the first layer of $1^{st}$ and $2^{nd}$ MuCh station (left panel) and $3^{rd}$ and $4^{th}$ MuCh station (right panel).}
\label{fig:occupancy_sector}
\end{center}
\end{figure}

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

We have presented a time-based simulation scheme for the MuCh detector of the CBM experiment which provides a realistic detector response of detector and read-out electronics in the environment of a trigger-less, free running data acquisition system. The key to this scheme is to enable interference of tracks from different events with small time separation, using an intelligent buffering procedure and a proper prescription for the treatment of signals arriving close in time in the same read-out channel (pad). The described MuCh simulation software is integrated into the common cbmroot simulation framework, which produces a data stream similar to that expected from the real experiment. The treatment of thermal, uncorrelated noise is included.

The simulation software allows to study the performance of the CBM-MuCh detector in a realistic scenario of a time sequence of minimum-bias nuclear collisions. In particular, rate-dependent effects like data losses arising from pile-up between subsequent events can be addressed, which is essential to assess the reconstruction efficiency of the detector under various operating conditions. Our findings of pile-up as functions of dead time and of the interaction rate are in qualitative agreement with our expectations, which to a certain extent verifies the software implementation. Quantitatively, we see significant losses from inter-event pile-up, above the irreducible in-event double hits, for interaction rates from $10^6$/s on.

The developed software also allows to optimize the read-out of the detector in terms of parameters like single-channel dead time, time resolution or the applied threshold, thus providing feedback to the designers of the read-out ASIC SMX. The software will be validated against the real detector behaviour using data from laboratory and in-beam tests of detector prototypes. In particular, we plan to improve the currently simplified treatment of double hits in an ASIC channel, replacing the now used Theta function with the real signal shape in the slow channel of the SMX.

The quantitative effects of the rate-dependent data losses on the sensitivity of the detector to physics observables will be the subject of further studies. This requires the application of the full, time-based reconstruction (hit and track finding) on the generated data stream. The simulation software described in the article provides the required tools for such studies.