\section{Introduction}
\label{sec:intro}
The investigation of highly redshifted absorption systems in the line-of-sight of quasars, by means of high-resolution spectroscopic observations, has become an established and powerful method to constrain a possible variation of the laws of physics over cosmological time. The variation of dimensionless fundamental constants, appearing as key building blocks in the laws of physics, is targeted in a comparison between astrophysical observation with laboratory measurement. For example, the fine-structure constant $\alpha=e^{2}/(4\pi\epsilon_{0}\hbar c)$ is probed via metal absorption searching for temporal \mbox{\citep{Webb1999}} and spatial variations \mbox{\citep{Webb2011}}.

The other dimensionless constant of nature determining the structure of molecular matter, the proton-to-electron mass ratio \mbox{$\mu \equiv m_{p}/m_{e}$}, may be targeted by observing absorption lines of a wide variety of molecules \mbox{\citep{Jansen2014}}. While radio astronomical observations of ammonia~\citep{Murphy2008, Kanekar2011} and of methanol molecules~\citep{Bagdonaite2013a,Bagdonaite2013b,Kanekar2015} have proven to be very sensitive, constraining $|\Delta\mu/\mu|$ at the level of \mbox{$\sim 10^{-7}$}, there are are only two radio sources at extragalactic distances where these molecules are observed, and these lie at redshifts $z<1$.

Molecular hydrogen (H$_{2}$) is a target to investigate a possible variation of $\mu$ at higher redshifts, where it is  observed in larger numbers of absorption systems. H$_{2}$ has many spectral lines that are sensitive to a variation in $\mu$~\citep{Thompson1975,Varshalovich1993,Ubachs2007}. Molecular hydrogen absorbing galaxies often display up to 100 H$_{2}$ spectral lines, which helps in improving the statistical basis of measurements made with this technique. Most importantly, H$_{2}$ absorbing systems can be observed, using ground-based telescopes, at redshifts $z>2$, even extending to $z=4.22$~\citep{Bagdonaite2015}, where the Lyman and Werner bands are redshifted into the optical band.

So far eight H$_{2}$ absorbing systems at $z>2$ have been analyzed for $\mu$-variation, most notably and most accurately the \mbox{J2123--005} system at \mbox{\emph{z}=2.05}~\citep{Malec2010,Weerdenburg2011}, the \mbox{HE0027--1836}~\citep{Rahmani2013} and the \mbox{Q2348--011} system~\citep{Bagdonaite2012} both at \mbox{\emph{z}=2.42}, the \mbox{Q0405--443} system at \mbox{\emph{z}=2.59}~\citep{Reinhold2006,King2008,Thompson2009}, the \mbox{Q0642--504} system at \mbox{\emph{z}=2.66}~\citep{Bagdonaite2014,Vasquez2014}, the \mbox{Q0528--250} system at \mbox{\emph{z}=2.81}~\citep{Reinhold2006,King2008,King2011}, the \mbox{Q0347--383} system at \mbox{\emph{z}=3.02}~\citep{Reinhold2006,King2008,Thompson2009,Wendt2011,Wendt2012}, and the \mbox{Q1443+272} system at \mbox{\emph{z}=4.22}~\citep{Bagdonaite2015}. A general conclusion from these studies is that the proton-to-electron mass ratio is constrained to \mbox{$|\Delta\mu/\mu| < 10^{-5}$} for redshifts in the range \mbox{\emph{z}=2-4.2}, corresponding to look-back times of up to 90 per cent of the age of the Universe. A review of the results of the analyses of H$_{2}$ absorbing systems has been given by \cite{Ubachs2011} and by \cite{Wendt2014}.

In the present study H$_{2}$ absorption at \emph{z}=2.69 is investigated in the line of sight toward quasar \mbox{SDSS J123714.60+064759.5}, hereafter \mbox{J1237+0647}, for further constraining a variation of the proton-to-electron mass ratio. The importance of this system is that carbon monoxide (CO) absorption in the optical domain is observed alongside with H$_2$~\citep{Noterdaeme2010}. This may be the basis for a simultaneous $\mu$-variation analysis using H$_2$ and CO~\citep{Salumbides2012}. As a first step the H$_2$ absorption in this system is analyzed for constraining $\Delta\mu/\mu$, results of which are presented in this paper, and this will serve as the basis for an analysis of the CO lines in a later paper.

\section{Data}
\label{sec:data}
The analysis presented in this work is based on observations made under four different observing programs, all carried out using the Ultraviolet and Visual Echelle Spectrograph (UVES) mounted on the 8.2m Very Large Telescope (VLT) at Paranal, Chile \citep{Dekker2000}. UVES is a cross-dispersed echelle spectrograph with two arms that are functionally identical: one covers the wavelengths in the range \mbox{300--500 nm} (Blue) and the other covers the range \mbox{420--1100 nm} (Red). Two dichroic beam splitters can be used to work in parallel with the two arms, centering them at different wavelengths. The details of the observational campaigns are presented in Table \ref{dataset}.

The raw 2D data were reduced following the procedure used by \cite{Bagdonaite2014}. The Common Pipeline Language version of the UVES pipeline was used to bias correct and flat field the exposures and then to extract the flux.  The ThAr lamp exposures were used for wavelength calibration; the ThAr flux was extracted  by using the same object profile weights as in a corresponding quasar exposure. After the standard reduction, the custom software \textsc{UVES\_popler}\footnote{\url{http://astronomy.swin.edu.au/~mmurphy/UVES_popler/}} was used to combine the extracted echelle orders into a single 1D spectrum which was then manually inspected and cleaned from bad pixels and other spectral artifacts, and the continuum was fitted with low-order polynomials.

The `final spectrum', obtained by combining all the exposures together, covers the wavelengths from 329.0 to \mbox{960.0 nm} with a signal-to-noise ratio (S/N) of \mbox{26 per 2.5\,km\,s$^{-1}$ per pixel} at \mbox{$\sim$400 nm} in the continuum.

\begin{table*}
\centering
\caption{Details of the J1237+0647 quasar observations with UVES/VLT included in the present analysis. The total integration time, summed over 17 exposures, is 19.9 hrs. Most of the listed exposures were followed by attached ThAr calibrations and additional `supercalibrations' as indicated (see Section \ref{sec:supercali}). The rest of the data were calibrated using the regular ThAr exposures taken at the end of the night. The CCD binning was \mbox{2$\times$2} and the slit width was {1.0 arcsec} for all frames. Data collected under programs 082.A-0544(A) and 083.A-0454(A) were retrieved from the ESO archive.}
\label{dataset}
\begin{tabular}{*{6}{c}}
\hline
Program ID & Date & Integration time [s] & Central $\lambda$ [nm] & ThAr & Supercalibration\\
\hline
\hline
082.A-0544(A) & 27-03-2009 & 5400 & 390 + 564 & Regular & No\\
 & 27-03-2009 & 5400 &  & Regular & No\\
	& 29-03-2009 & 5400 & & Regular & No\\
 & 29-03-2009 & 5400 & & Regular & No\\
083.A-0454(A) & 27-04-2009 & 4500 & 390 + 775 & Regular & No\\
 & 27-04-2009 & 4500 & & Regular & No\\
091.A-0124(A) & 15-05-2013 & 4800 & 390 + 580 & Attached & Yes\\
 & 15-05-2013 & 1727 & & Attached & No\\
 & 15-05-2013 & 4800 & & Attached & Yes\\
 & 16-05-2013 & 4800 & & Attached & Yes\\
 & 16-05-2013 & 2025 & & Attached & No\\
093.A-0373(A) & 23-03-2014 & 4800 & 390 + 580 & Attached & Yes\\
 & 03-04-2014 & 4800 & & Attached & Yes\\
 & 28-05-2014 & 4800 & & Attached & Yes\\
 & 29-05-2014 &  826 & & Attached & No\\
 & 31-05-2014 & 4800 & & Attached & Yes\\
 & 02-06-2014 & 2800 & & Attached & No\\
\hline
\end{tabular}
\end{table*}

\subsection{Data from 2013 and 2014}
\label{subsec:2013/4}
The quasar J1237+0647 was observed in visitor mode in May 2013, program 091.A-0124(A), and in service mode in the period March-June 2014, program 093.A-0373(A). The total integration time amounted to \mbox{11.5 hrs}. Each exposure was nominally 4800 s in duration and was expected to deliver a S/N of \mbox{13 per 2.5\,km\,s$^{-1}$ per pixel} at \mbox{400 nm} and a resolving power \emph{R} of \mbox{$\sim$40~000}, with a seeing of \mbox{0.8 arcsec}, a slit of \mbox{1.0 arcsec}, an airmass of 1.2 and \mbox{2$\times$2} binning. Due to bad weather conditions, four exposures were shorter than scheduled and, therefore, the predicted S/N level could not be reached in those cases.

Each of the quasar exposures had an `attached' standard ThAr calibration, which is a lamp exposure taken immediately after the science exposure leaving all the instrument parameters unchanged, and was immediately followed by a `supercalibration' exposure without allowing any grating reset. The supercalibration is a method to quantify any long-range wavelength distortion and it involves observing objects with a solar-like spectrum, like asteroids or `solar twin' stars (see Section \ref{sec:supercali}).

\subsection{Data from 2009}
\label{subsec:2009}
Part of the total dataset used in this work was retrieved from the ESO data archive\footnote{\url{http://archive.eso.org/eso/eso_archive_main.html}}. This was the case for observations from programs 082.A-0544(A) and 083.A-0454(A), which were carried out in March-April 2009 and were reported by \cite{Noterdaeme2010}. The main difference with the most recent observations of J1237+0647 is that these exposures only had the regular `non-attached' ThAr calibrations taken at the end of each night and dedicated observations of supercalibration targets were not carried out. These observations of \mbox{J1237+0647} from 2009 added \mbox{8.5 hrs} of integration on the target.

\section{Method}
\label{sec:method}
\subsection{Theory}
\label{subsec:theory}
\cite{Thompson1975} originally proposed that a possible temporal variation of the proton-to-electron mass ratio \mbox{$\mu = m_{p}/m_{e}$} can be detected using absorption spectra of molecular hydrogen observed at high redshift. The observed wavelength $\lambda^{z}_{i}$ of the $\emph{i}$-th transition will show a shift given by:
\begin{equation}
\lambda^{z}_{i}=\lambda^{0}_{i}(1+\ensuremath{z_\textrm{\scriptsize abs}})(1+K_{i}\frac{\Delta\mu}{\mu}),
\label{wl_shift}
\end{equation}
where $\lambda^{0}_{i}$ is the rest wavelength, $\ensuremath{z_\textrm{\scriptsize abs}}$ is the redshift where the absorption occurs, $\Delta\mu \equiv \mu_{z} - \mu_{0}$ is the difference between the proton-to-electron mass ratio in the system observed and the value measured in the laboratory, and $K_{i}$ is the sensitivity coefficient which, for a given transition, determines the shifting power and sign due to varying $\mu$. These coefficients are defined as
\begin{equation}
K_{i} = \frac{\mathrm d \ln\lambda_{i}}{\mathrm d \ln\mu}.
\label{kcoeff}
\end{equation}
The sensitivity coefficients for the H$_{2}$ molecule used in this analysis were calculated within a semi-empirical framework by \cite{Ubachs2007}, including effects beyond the Born-Oppenheimer approximation.

The H$_{2}$ laboratory wavelengths, $\lambda^{0}_{i}$, used as a reference in this analysis were measured by \cite{Salumbides2008} with a fractional wavelength accuracy $\Delta\lambda/\lambda\sim5\times10^{-8}$, and they can be considered exact in comparison with the uncertainties of the lines in the quasar spectrum. The complete list of laboratory data required for fitting the H$_{2}$ absorption lines observed in this work, including rest wavelengths $\lambda^{0}_{i}$, oscillator strengths \emph{f}$_{i}$, damping coefficients $\Gamma_{i}$ and the sensitivity coefficients \emph{K}$_{i}$, was compiled by \cite{Malec2010}.

\subsection{Fitting Method}
\label{subsec:fitmethod}
A comprehensive fitting method \mbox{\citep{Malec2010,King2011,Bagdonaite2014}}, which involves a simultaneous treatment of all the lines, was used to model the H$_{2}$ spectrum. Furthermore this technique allows to tie part of the fitting parameters together, resulting in a smaller number of free parameters. Simultaneous fitting and parameter tying allows to include also those H$_{2}$ transitions that are overlapped by metal or \mbox{H \textrm{\sc{i}} }lines from the Lyman $\alpha$ forest which can be modelled at the same time.

The program used is the non-linear least-squares Voigt profile fitting program \textsc{vpfit}\footnote{\url{http://www.ast.cam.ac.uk/~rfc/vpfit.html}}, developed specifically for quasar spectra analysis. A Voigt profile is obtained  from the convolution of a Lorentzian profile, which describes the natural line broadening of the transitions and is specific for each molecular transition considered, and a Gaussian profile, which reflects the physical conditions within the absorbing clouds by describing the velocities, both thermal and turbulent, of the molecules. These two profiles are eventually convolved with a model for the instrumental profile, which is assumed to be Gaussian.

In \textsc{vpfit}, the profile of each velocity component, of each transition, is described by three free parameters: the column density \emph{N}, the redshift \emph{z} at which the absorption occurs, and the Doppler width \emph{b}. These parameters are used in addition to the molecular physics input for every \emph{i}-th transition, like $\lambda^{0}_{i}$, \emph{f}$_{i}$, whose product with \emph{N} gives the optical depth of the absorption line, and $\Gamma_{i}$, which defines the lifetime broadening. Transitions probing the same rotational ground state share the same population, and as a consequence they are all described using the column density \emph{N}$_{J}$. An important underlying assumption is that absorbing features at the same redshift originate from the same cloud, hence they share certain parameters related to the physical conditions of the cloud. In particular, they are assumed to share the redshift \emph{z} and the width \emph{b}, which reflects the turbulent motion and the kinetic temperature of the gas in the cloud. This is achieved by tying in \textsc{vpfit} the \emph{z} and \emph{b} parameters of all the absorption lines associated with the gas in a single absorbing cloud.

The fitting program finds the best match between the data and the model iteratively by minimizing the $\chi^{2}$ parameter. The initial values of the free parameters are inserted by the user. Then, in each iteration, \textsc{vpfit}  changes the values of the free parameters of the model and checks the changes in relative $\chi^{2}$, reporting convergence once a stopping criterion, which is user-defined, is met.

Since introducing too many free parameters may lead to an overfit of the spectrum, models with a different number of velocity components per H$_{2}$ transition were compared using the $\chi^{2}_{\nu}$ and the Akaike Information Criterion~\citep[\emph{AICC},][]{Akaike1974} parameter. In particular, the latter parameter is defined as
\begin{equation}
AICC = \chi^{2} + 2p + \frac{2p(p+1)}{n-p-1}
\end{equation}
where \emph{p} is the number of free parameters and \emph{n} is the number of spectral points included in the fit. $\Delta$\emph{AICC}$>5$ and $\Delta$\emph{AICC}$>10$ are considered to be a strong and a very strong evidence respectively that the model with the lower \emph{AICC} is statistically preferred.

Multiple fits, with a different number of velocity components each, were performed in order to develop a robust absorption model. Only after the results of the fitting process are stable, the possible variation of the proton-to-electron mass ratio expressed by the variable $\Delta\mu/\mu$ is invoked as an extra free parameter in the fit besides \emph{N$_{J}$}, \emph{z} and \emph{b}. This fourth parameter is not introduced earlier during the fitting process to avoid that any over/underfitting can be compensated by an artificial variation of $\mu$. The program derives the value of $\Delta\mu/\mu$ according to Eq. (\ref{wl_shift}), therewith also using the sensitivity coefficients of Eq. (\ref{kcoeff}). The consistency of the derived $\Delta\mu/\mu$ constraint is then tested for sensitivity to a number of assumptions made in the fitting process. The results will be discussed in detail in Section \ref{subsec:constraint}.

\section{Analysis}
\label{sec:analysis}
Quasar J1237+0647 is located at redshift \emph{z}=2.78, thereby defining the extent of the \mbox{Lyman $\alpha$} forest, which comprises multiple series of neutral hydrogen transitions arising from the intergalactic medium at all redshifts below the quasar redshift. Toward the quasar there is one major absorption system at redshift $z=2.69$, which is a sub-damped \mbox{Lyman $\alpha$} system (DLA) with a neutral hydrogen column density of \mbox{$\log N=20.0\pm0.15$ cm$^{-2}$}, that features many absorption lines from atomic, ionized and molecular species \mbox{\citep{Noterdaeme2010}}.

Molecular hydrogen absorption lines associated with the DLA are spread in the range from 352.1 to \mbox{409.5 nm}. H$_{2}$ is found in three absorption features at redshift \mbox{\emph{z}=2.688001}, \mbox{\emph{z}=2.68868} and \mbox{\emph{z}=2.68955}; their velocities relative to the redshift of the strongest feature at \mbox{\emph{z}=2.68955} are respectively \mbox{\emph{v}=$-125$\,km\,s$^{-1}$}, hereafter Z1, \mbox{\emph{v}=$-72$\,km\,s$^{-1}$}, hereafter Z2, and \mbox{\emph{v}=0\,km\,s$^{-1}$}, hereafter Z3. The presence of these three features close to each other results in the line profiles shown in Fig. \ref{lprof}. The presence of clearly separated velocity features is not uncommon and has been observed in various systems in the past. The number of such components can vary from one up to the seven distinct features observed in the absorbing system toward quasar Q2348--011 \mbox{\citep{Petitjean2006,Bagdonaite2012}}.

Z3 has the highest column density for every \emph{J}-level, resulting in heavily saturated absorption features for low \emph{J}-levels with \mbox{\emph{J}$\le1$}. Furthermore, Z3 shows also deuterated molecular hydrogen (HD) absorption, where only the \mbox{\emph{J}=0 level} was detected \mbox{\citep{Noterdaeme2010}}. The HD molecule is sensitive to a variation of $\mu$ \mbox{\citep{Ivanov2008,Ivanov2010}}, and its detected transitions were included in the present $\mu$-variation analysis.
\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{velprofileJ3.pdf}
\caption{Typical line profile for the three absorption features toward J1237+0647. The transition plotted is L1P(3) at \mbox{$\lambda_{0}=109.978$ nm}. The vertical dashed lines (red) indicate the positions of the three H$_{2}$ spectral features fitted in this study. The velocity frame is centered on the strongest feature Z3.}
\label{lprof}
\end{figure}

\subsection{Creating an absorption model}
\label{subsec:model}

\subsubsection{Selecting and fitting spectral regions}
\label{ssubsec:regions}
Given the broad absorption profile, spanning \mbox{$\sim150$\,km\,s$^{-1}$}, overlaps between profiles of different molecular hydrogen transitions are common as well as overlaps with the H$\,$\textsc{i} transitions from the Lyman $\alpha$ forest. Within the comprehensive fitting method, all kinds of overlaps can be handled. However, no relevant information was gained in the case of a complete overlap with a saturated H$\,$\textsc{i} line, hence such absorption features were not included in this study. Less frequently, some H$_{2}$ features partially overlap with narrow metal transitions. However, these lines were not excluded from the dataset considered in this work. A large number of metal absorption features were identified by \mbox{\cite{Noterdaeme2010}}.

All the 137 potentially useful H$_{2}$ and HD transitions were contained in 66 spectral regions, in the range \mbox{351.0--413.5 nm}, which were selected trying to avoid, whenever possible, \mbox{H$\,$\textsc{i}} transitions and too-heavily saturated low \mbox{\emph{J}-level} H$_{2}$ absorption features in Z3. The useful H$_{2}$ transitions, probing six different rotational quantum states with \emph{J}=0--5, are listed in Table \ref{transitions}. Absorption features originating at the same redshift had both the \emph{z} and the \emph{b} parameters tied together, while the column densities \emph{N}$_{J}$ were allowed to vary independently from each other. Z3 shows saturated H$_{2}$ lines, particularly in the low \emph{J}-levels, but there is an apparent evidence of a partial coverage with a non-zero residual flux at the base of its absorption features. To take into account this effect, a zero-level correction was included as a free parameter in the fit for each H$_{2}$ region.
\begin{table*}
\centering
\caption{List of the molecular transitions used in this analysis.}
\label{transitions}
\begin{tabular}{cp{7.5cm}p{5cm}c}
\hline
\emph{J}-level & Lyman transitions & Werner transitions & n$_{\textrm{transitions}}$ \\
\hline
\hline
H$_{2}$ \emph{J} = 0 & L0R(0), L2R(0), L3R(0), L4R(0), L6R(0), L7R(0), L8R(0), L10R(0) & W0R(0), W1R(0) & 10 \\
H$_{2}$ \emph{J} = 1 & L0P(1), L0R(1), L1P(1), L1R(1), L2P(1), L2R(1), L3P(1), L7P(1), L9R(1), L10R(1), L12R(1) & W0R(1), W1R(1) & 13 \\
H$_{2}$ \emph{J} = 2 & L0P(2), L0R(2), L1P(2), L1R(2), L2P(2), L2R(2), L3P(2), L4P(2), L4R(2), L5P(2), L6P(2), L7P(2), L7R(2), L8P(2), L9P(2), L10P(2), L10R(2), L11P(2), L11R(2), L13P(2), L13R(2) & W0P(2), W0Q(2), W1Q(2), W2P(2), W2Q(2), W2R(2) & 27 \\
H$_{2}$ \emph{J} = 3 & L0P(3), L0R(3), L1P(3), L1R(3), L2R(3), L3P(3), L3R(3), L5P(3), L5R(3), L6P(3), L6R(3), L7P(3), L7R(3), L8P(3), L8R(3), L9R(3), L10P(3), L12R(3), L13P(3), L13R(3), L14P(3) & W0P(3), W0Q(3), W1R(3), W2P(3), W2Q(3), W2R(3), W3P(3) & 28\\
H$_{2}$ \emph{J} = 4 & L0P(4), L0R(4), L1P(4), L1R(4), L2R(4), L3P(4), L3R(4), L4P(4), L4R(4), L5R(4), L6P(4), L7P(4), L7R(4), L8P(4), L9P(4), L9R(4), L11P(4), L11R(4), L12R(4), L13P(4), L14P(4) & W0P(4), W0Q(4), W0R(4), W1P(4), W2Q(4), W2R(4), W3Q(4) &  28\\
H$_{2}$ \emph{J} = 5 & L0R(5), L1P(5), L2P(5), L2R(5), L3P(5), L3R(5), L4P(5), L5P(5), L5R(5), L7P(5), L8P(5), L9P(5), L10P(5), L11R(5), L12R(5), L13P(5), L14P(5), L15P(5) & W0P(5), W0R(5), W1Q(5), W2Q(5), W2R(5), W3P(5) & 24 \\
\hline
HD \emph{J} = 0 & L3R0, L5R0, L8R0, L12R0, L15R0 & W1R0, W3R0 & 7 \\
\hline
\multicolumn{3}{l}{Total number of lines} & 137 \\
\hline
\end{tabular}
\end{table*}

In order to build a robust model, the neutral hydrogen absorption lines occurring in the Lyman $\alpha$ forest and included in the selected regions must be accounted for. The \mbox{H$\,$\textsc{i}} transitions were modelled by assigning to each of them a set of free parameters (\emph{N}, \emph{z}, \emph{b}) in \textsc{vpfit}. Initial guesses of the free parameters were user-provided for each H$_{2}$, HD and \mbox{H$\,$\textsc{i}} transition.

\subsubsection{Extra velocity components}
\label{ssubsec:extravcs}
The H$_{2}$ absorption spectrum toward J1237+0647 displays three clearly distinct features in velocity space. In order to fit the spectrum as accurately as possible, it is imperative to consider possible underlying substructures  in velocity space. Physically this may relate to the existence of closely separated distinct clouds, or inhomogeneities in the clouds in the line-of-sight. To investigate such underlying velocity structure, the fitting residuals from the 66 selected regions were normalized, shifted to a common velocity scale and averaged together, creating a composite residual spectrum \mbox{\citep[CRS,][]{Malec2010}}. Only the H$_{2}$ absorption features that were not overlapped by \mbox{H$\,$\textsc{i}} features were included in the CRS, in order to avoid that the residuals of the \mbox{H$\,$\textsc{i}} model could introduce any spurious effect. The CRS, by combining the residuals of many transitions, highlighted any over/underfitted structure in the absorption features.

Z1 could be described using only one velocity component (VC), while indications of underfitting were observed in the CRS for absorption features Z2 and Z3, hence they were modelled by adding extra VCs. They were added to the model until \textsc{vpfit} started to reject them. Then the models with a different number of VCs  were compared using the $\chi^{2}$ and the \emph{AICC} parameters, in order to decide which one better fitted the data. The CRS for the three absorption features Z1, Z2 and Z3 are shown in Figs. \ref{res1}--\ref{res3}. Z2 could be modelled using two VCs, labelled as `a' and `b' in Fig. \ref{res2}, and Z3 was modelled using three VCs, labelled as `a', `b' and `c' in Fig. \ref{res3}. Parameters for the model are summarized in Table \ref{param} and the VCs are shown in Fig. \ref{abs1}. Previous studies found no strong evidences for a systematic shift in $\Delta\mu/\mu$ due to a segregation in different rotational levels \mbox{\citep{Malec2010,King2011}}, hence it was assumed that the velocity structure is the same among all the \emph{J}-levels in each absorption feature. The complete absorption model is shown in Section \ref{app:model}.

\begin{figure}
\includegraphics[width=\columnwidth]{resex1.pdf}
\caption{Top panel: normalized composite residual spectrum for Z1 formed from 30 H$_{2}$ absorption features. The dashed lines represent the $\pm 1\sigma$ boundaries. Bottom panel: an example transition shown on the same velocity scale. The velocity scale is centred at \mbox{$z = 2.688009$}. The vertical dashed line shows the position of the velocity component.}
\label{res1}
\end{figure}•
\begin{figure}
\includegraphics[width=\columnwidth]{resex2.pdf}
\caption{Top panel: normalized composite residual spectrum for Z2 formed from 33 H$_{2}$ absorption features obtained from the 2 VCs model. The dashed lines represent the $\pm 1\sigma$ boundaries. Bottom panel: an example transition shown on the same velocity scale. The velocity scale is centred at \mbox{$z = 2.688664$}, on the strongest velocity component. The vertical dashed lines show the positions of the two velocity components at $v = 0$\,km\,s$^{-1}$ (a) and +7\,km\,s$^{-1}$ (b).}
\label{res2}
\end{figure}•
\begin{figure}
\includegraphics[width=\columnwidth]{resex3.pdf}
\caption{Top panel: normalized composite residual spectrum for Z3 formed from 33 H$_{2}$ absorption features obtained from the 3 VCs model. The dashed lines represent the $\pm 1\sigma$ boundaries. Bottom panel: an example transition shown on the same velocity scale. The velocity scale is centred at \mbox{$z=2.689555$}, on the strongest velocity component. The vertical dashed lines show the positions of the three velocity components at $v = -4$\,km\,s$^{-1}$ (a), 0\,km\,s$^{-1}$ (b) and +13\,km\,s$^{-1}$ (c).}
\label{res3}
\end{figure}•

\begin{table*}
\centering
\caption{Column densities \emph{N$_{J}$}, redshift \emph{z} and Doppler parameters \emph{b} of the H$_{2}$ transitions in the three clearly distinguishable absorption features and in their underlying velocity components a, b, c. Levels with \emph{J} = 4 and 5 are too weak to be detected in Z1 and Z2. Levels with \emph{J} = 0 and 1 are heavily saturated in Z3 and they have been discarded from the model.}
\label{param}
\begin{tabular}{c|c|cc|ccc}
 & Z1 & \multicolumn{2}{c}{Z2} & \multicolumn{3}{c}{Z3} \\
 \hline
\emph{J}-levels & $\log{N_{J}}$ [cm$^{-2}$] & \multicolumn{2}{c}{$\log{N_{J}}$ [cm$^{-2}$]} & \multicolumn{3}{c}{$\log{N}_{J}$ [cm$^{-2}$]} \\
 & & a & b & a & b & c \\
\hline
\emph{J} = 0 & $15.74 \pm 0.11$ & $15.44 \pm 0.11$ & $13.78 \pm 0.17$ & \multicolumn{3}{c}{Discarded} \\
\emph{J} = 1 & $16.29 \pm 0.08$ & $16.23 \pm 0.11$ & $14.56 \pm 0.10$ & \multicolumn{3}{c}{Discarded} \\
\emph{J} = 2 & $15.35 \pm 0.04$ & $15.27 \pm 0.07$ & $13.90 \pm 0.15$ & $15.67 \pm 0.10$ & $17.64 \pm 0.10$ & $15.26 \pm 0.16$ \\
\emph{J} = 3 & $15.21 \pm 0.02$ & $14.88 \pm 0.05$ & $14.31 \pm 0.08$ & $15.66 \pm 0.11$ & $18.02 \pm0.06$ & $14.64 \pm 0.28$ \\
\emph{J} = 4 & \multicolumn{3}{c}{N/D} & $14.79 \pm 0.07$ & $15.34 \pm 0.04$ & $13.93 \pm 0.07$ \\
\emph{J} = 5 & \multicolumn{3}{c}{N/D} & $14.50 \pm 0.07$ & $14.79 \pm 0.04$ & $13.33 \pm 0.17$ \\
HD \emph{J} = 0 & \multicolumn{3}{c}{N/D} & N/D & $14.36 \pm 0.03$ & N/D \\
\hline
\emph{z} & 2.68801(06) & 2.68866(09) & 2.68872(82) & 2.68950(96) & 2.68955(14) & 2.68972(52) \\
\hline
\emph{b} [kms$^{-1}$] & $2.94 \pm 0.08$ & $2.32 \pm 0.13$ & $6.99 \pm 0.65$ & $11.36 \pm 0.48$ & $5.78 \pm 0.11$ & $2.54 \pm 0.48$ \\
\hline
\end{tabular}
\end{table*}

\begin{figure*}
\centering
\includegraphics[width=2\columnwidth]{vcs.pdf}
\caption{Example transitions for every \emph{J}-level. Each velocity component is plotted against the fitted continuum. \mbox{H$\,$\textsc{i}} contributions to the line shapes are not highlighted. Top panel: \emph{J}-levels 0-4 for Z1 (blue), the velocity scale is centred at \emph{z}=2.688009. Middle panel: \emph{J}-levels 0-4 for Z2. The velocity component at \emph{z}=2.688664, Z2a, is depicted with the (blue) solid line and the one at \emph{z}=2.688717, Z2b, with the (red) dashed line. The velocity scale is centred on Z2a. Bottom panel: \emph{J}-levels 2-5 for Z3. The three velocity components at redshifts \emph{z}=2.689498, Z3a, 2.689555, Z3b, and 2.689723, Z3c, are shown with the dashed (red), solid (blue) and dotted (green) line respectively. The velocity scale is centred on Z3b.}
\label{abs1}
\end{figure*}

\subsection{Constraining $\Delta\mu/\mu$}
\label{subsec:constraint}
Once an optimal model was found, an extra free parameter, corresponding to $\Delta\mu/\mu$, was added to the H$_{2}$ transitions in \textsc{vpfit}. All the transitions were included in the fit, without any distinction between their bands, their \emph{J}-levels or absorbing clouds where they arose. The model returns \mbox{$\Delta\mu/\mu=(-1.1 \pm 6.3_{\textrm{\small{stat}}})\times10^{-6}$}. The statistical error reported here was derived from the appropriate diagonal term of the final covariance matrix for the fit. Given the absorption model, it represents only the uncertainty in $\Delta\mu/\mu$ stemming from the photon statistics of the quasar spectrum, i.e. its signal-to-noise ratio. However, to evaluate the robustness of the result, hereafter the statistical value of $\Delta\mu/\mu$, a variety of consistency tests and sources of systematic error which contribute to the overall uncertainty budget were explored, results of which are shown in Fig. \ref{tests} and discussed below.

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{tests1.pdf}
\caption{$\Delta\mu/\mu$ constraints obtained from various tests (triangles) compared with the statistical value from the entire dataset (circle). The solid (green) horizontal line indicates the value of the statistical constraint and the two dashed lines represent its $1 \sigma$ statistical uncertainties. The square represents the $\Delta\mu/\mu$ value obtained only from the exposures taken in 2014 after applying to each of them a long-range distortion correction, as discussed in Section \ref{sec:supercali}. The diamond (red) represents the fiducial value obtained for the entire dataset (2009+2013+2014) after applying the distortion corrections.}
\label{tests}
\end{figure}

\subsubsection{Isolating single absorbing features}
\label{ssubsec:singleabs}
Each of the absorbing features Z1, Z2 and Z3 can be investigated separately. The constraints on $\Delta\mu/\mu$ returned are \mbox{$\Delta\mu/\mu\vert_{\textrm{\small{Z1}}}=(11.4 \pm 10.2_{\textrm{\small{stat}}})\times10^{-6}$} for Z1, \mbox{$\Delta\mu/\mu\vert_{\textrm{\small{Z2}}}=(-17.7 \pm 14.0_{\textrm{\small{stat}}})\times10^{-6}$} for Z2 and \mbox{$\Delta\mu/\mu\vert_{\textrm{\small{Z3}}}=(-6.1 \pm 14.8_{\textrm{\small{stat}}})\times10^{-6}$} for Z3. The resulting constraints are in agreement with each other and with the statistical value for the entire dataset within $1 \sigma$ boundaries.

\subsubsection{Isolating Lyman and Werner transitions}
\label{ssubsec:lymanwerner}
All the H$_{2}$ transitions considered in this analysis belong either to the Lyman- or to the Werner band systems. Any long-range distortion of the wavelength scale will be nearly degenerate with a shift of the molecular hydrogen lines due to a non-zero value of $\Delta\mu/\mu$. This degeneracy is partly broken by fitting together H$_{2}$ transitions from both Lyman and Werner bands. The effect of such a distortion can be investigated by fitting separately transitions from the two band systems.

The values returned in this way are \mbox{$\Delta\mu/\mu\vert_{\textrm{\small{L}}} = (-0.6 \pm 6.4_{\textrm{\small{\small{stat}}}}) \times 10^{-6}$} for Lyman transitions and \mbox{$\Delta\mu/\mu\vert_{\textrm{\small{W}}} = (-10.3 \pm 33.9_{\textrm{\small{stat}}}) \times 10^{-6}$} for Werner transitions. The larger error of the value returned by Werner transitions is due to three effects. Firstly, there are fewer Werner transitions than Lyman transitions in the spectrum, with a ratio between the former and the latter of $\sim$1:4. Then the spread of the sensitivity coefficients is larger for Lyman transitions, \mbox{$\Delta K_{\textrm{\small{L}}} = 0.092$}, than the one of Werner transitions, \mbox{$\Delta K_{\textrm{\small{W}}} = 0.054$}. Finally, the Werner transitions fall in the bluest part of the spectrum, where the S/N is much lower than in the region where Lyman transitions are. Such differences lead the statistical value of $\Delta\mu/\mu$ to be dominated by the Lyman transitions.

\subsubsection{Isolating cold and warm states}
\label{ssubsec:coldwarm}
Using \textsc{vpfit} it is possible to calculate a different value of $\Delta\mu/\mu$ for transitions originating in each \emph{J}-level, in order to obtain the relative contributions of different levels to the statistical value. This is particularly useful to investigate the impact of possible temperature inhomogeneities in the absorbing clouds. Due to the para-ortho distribution, the \emph{J}$=1$ rotational state is significantly populated even at the lowest temperatures \mbox{\citep{Ubachs2007}}. As a consequence, the H$_{2}$ transitions can be divided into two sets: rotational states with \emph{J}=0,1 (cold) and with \emph{J}$\ge2$ (warm).

The constraints returned by the model are \mbox{$\Delta\mu/\mu\vert_{\textrm{\small{cold}}} = (-15.1 \pm 11.5_{\textrm{\small{stat}}}) \times 10^{-6}$} for the cold states and \mbox{$\Delta\mu/\mu\vert_{\textrm{\small{warm}}} = (4.4 \pm 7.1_{\textrm{\small{stat}}}) \times 10^{-6}$} for the warm states. The two values match within their uncertainties, which means that any temperature inhomogeneity in the absorbing clouds does not have a significant impact on the $\mu$-variation.

Furthermore, the slightly larger error derived from the cold states reflects the fact that levels with \emph{J}$\le3$ are saturated in Z3 and, in particular, cold states are heavily saturated. As a consequence of the saturation, there are fewer constraining pixels contributing to the fitting process and there are fewer transitions in the cold states dataset.

\subsubsection{Separating exposures from 2009}
\label{ssubsec:20092014}
The dataset used in this work can be divided between the exposures taken in 2013 and 2014, which had an attached ThAr calibration and a supercalibration, and the exposures taken in 2009, retrieved from the ESO archive, which only had regular non-attached ThAr calibrations. In order to see how the absence of the individual calibrations affects the statistical value, two values for $\Delta\mu/\mu$ were derived from the two subsets.

The constraints obtained in this way are \mbox{$\Delta\mu/\mu\vert_{\textrm{\small{2009}}} = (-9.1 \pm 7.6_{\textrm{\small{stat}}}) \times 10^{-6}$} for the exposures taken in 2009 and \mbox{$\Delta\mu/\mu\vert_{\textrm{\small{2014}}} = (14.3 \pm 10.3_{\textrm{\small{stat}}}) \times 10^{-6}$} for the exposures taken in 2013 and 2014. The two values do not agree within their combined 1$\sigma$ statistical uncertainties, although they do within \mbox{$1.5 \sigma$}.
However, this comparison does not take into account the long-range distortion effect, which will be discussed in the next section.

\section{supercalibrations}
\label{sec:supercali}
In order to constrain a variation of $\mu$, it is crucial to accurately wavelength calibrate the quasar observations. Attempts to perform advanced tests of the effective accuracy of UVES were performed by \mbox{\cite{Molaro2008}} by comparing a UVES-recorded asteroid spectrum with a highly accurate solar reference spectrum. These measurements revealed zero offsets up to \mbox{$\approx$ 50 ms$^{-1}$}, likely due to a non-uniform slit illumination, but found no evidence for a wavelength-dependent velocity shift in the spectrum. Similar studies, that are now referred to as `supercalibrations', showed that spectra taken with UVES and HIRES spectrographs suffer from intra-order calibration shifts on scales of single echelle orders but they did not uncover any evidence for a long-range distortion in the spectra on the scale of the spectrographs arms \mbox{\citep{Whitmore2010,Griest2010}}. \mbox{\cite{Rahmani2013}} detected for the first time a long-range distortion in some UVES spectra up \mbox{to $\sim$ 400 ms$^{-1}$} on a scale \mbox{of $\sim$ 600 \AA}. More recently, \mbox{\cite{Whitmore2015}} found that significant wavelength calibration distortions, both intra-order and long-range, are ubiquitous across UVES and HIRES history, most likely due to a different light path of the astronomical object and the ThAr lamp within the instrument.

The presence of long-range distortions between the attached ThAr calibration and the main target spectrum velocity scales introduces a wavelength-dependent shift that mimics a non-zero $\Delta\mu/\mu$. This happens because the sensitivity coefficients \emph{K$_{i}$} decrease with increasing wavelength (see Fig. 1 of \mbox{\cite{Bagdonaite2014}}) and a long-range distortion of the wavelength scale would produce a strong systematic effect which will be nearly degenerate with the effect produced by a variation of $\mu$. In principle, this degeneration can be broken by fitting Lyman and Werner transitions together; however, the paucity of Werner transitions and the low S/N in the spectral region where they fall prevent the degeneracy to be broken in this way, as illustrated by \mbox{\cite{Malec2010}} and discussed in Section \ref{ssubsec:lymanwerner}.
Over most of UVES's history, the sign of measured long-range distortions was such that they spuriously pushed $\Delta\mu/\mu$ to more positive values \mbox{\citep{Whitmore2015}}.

\subsection{Method}
\label{subsec:supercalimethod}
The first supercalibration performed with UVES consisted of a comparison between asteroid spectra and a Fourier Transform Spectrometer (FTS) solar atlas \mbox{\citep{Molaro2008}}. \mbox{\cite{Rahmani2013}} cross-correlated asteroid spectra with a laboratory solar spectrum and, more recently, \mbox{\cite{Whitmore2015}} improved the technique by forward modelling the FTS spectra to match observed asteroid spectra, thereby allowing information to be derived from much shorter spectral ranges, and by showing that the same technique can be applied on solar twin stars as well. Their method is used in this study to supercalibrate the quasar spectrum.

The supercalibration process consists in the comparison of a ThAr calibrated spectrum with a reference spectrum from a FTS, whose frequency scale is expected to be much more accurate than the UVES one, hence it is considered `absolutely calibrated' for the purpose of this study. The reference used is the FTS vacuum solar spectrum from \cite{Chance2010}\footnote{Available at \url{http://kurucz.harvard.edu/sun/irradiance2005/irradthu.dat}}. The targets for the supercalibrations were asteroids, which reflect the solar light and hence have the same spectrum as the Sun, and `solar twin' stars, i.e. objects whose spectra are almost identical to the solar one \mbox{\citep{Melendez2009,Datson2014}}. Since the supercalibration targets are astronomical objects, no changes are required in the focus or in the slit alignment.

The two spectra were divided into small regions \mbox{of $\approx$ 500\,km\,s$^{-1}$} each, which corresponds to \mbox{$\approx$ 8 \AA}. Their small size allows for sampling each echelle order with \mbox{$\approx$ 10} regions, and for identifying any relative velocity distortion down to \mbox{$\sim$ 30 ms$^{-1}$} \mbox{\citep{Whitmore2015}}. The spectral regions of the two spectra were compared using a $\chi^{2}$-minimization technique looking for any wavelength-dependent calibration distortions across each echelle order and across the wavelength range covered by all the echelle orders collectively. The overall velocity distortion was obtained by a linear fit of all the spectral regions and it was used to distortion-correct the corresponding quasar exposure.

\subsection{Supercalibration data}
\label{subsec:supercalidata}

\subsubsection{Data from 2013 and 2014}
\label{subsubsec:supercali2014}
The exposures taken in 2013 and in 2014 were recorded with an attached ThAr wavelength calibration. Since the problems of long-range distortions and their effect on $\mu$-variation analysis from H$_2$ spectra had been well identified and reported \mbox{\citep{Rahmani2013,Bagdonaite2014}}, supercalibration exposures targeting asteroids and solar twins were purposely recorded immediately after the science exposures with their attached ThAr calibrations. In view of their accessibility during observations of J1237+0647, supercalibration spectra of the Eunomia asteroid, as well a solar twin stars HD~76440, HD~147513, HD~097356 and HD~117860, were used. Only some exposures, which had to be aborted due to bad weather conditions, were not supported by a supercalibration.

For each supercalibration spectrum, its long-range wavelength distortion was determined in terms of a slope parameter, which was found to range between 300 and \mbox{725 ms$^{-1}$} per \mbox{1000 \AA}. In particular, during each observing run the values of the distortion slopes showed an apparent stability, with a spread in the slope values of, at most, \mbox{$\sim$200 ms$^{-1}$} per \mbox{1000 \AA} in May 2014. These values, and in particular the positive sign, are commensurate with values found previously in such measurements \mbox{\citep{Rahmani2013,Whitmore2015}}. Results of the present analysis are plotted in \mbox{Fig. \ref{supercals}} and numerical values listed in \mbox{Table \ref{supercaldataset}}. Thereupon each individual J1237+0647 science exposure was corrected for its long-range distortion by introducing a counter-distortion based on these slopes, after which the exposures were rebinned and combined using \textsc{UVES\_popler}. Exposures without a directly attached supercalibration were corrected using the average value of all the distortion slopes from the supercalibrations taken during the same night, in view of the stability discussed above.

\subsubsection{Data from 2009}
\label{subsubsec:supercali2009}
The data for J1237$+$0647 taken in 2009, retrieved from the ESO archive, were not supported by an attached ThAr calibration, and they were therefore wavelength calibrated against a regular ThAr calibration recorded at the end of the night. No supercalibration spectra were recorded with the J1237+0647 science exposures. In the ESO archive some suitable asteroid observations were found as three exposures of the Ceres asteroid, recorded in March and April 2009 within one week of the quasar exposures, under program 080.C-0881(B) (PI Dumas). These exposures did not share the same telescope settings as the J1237+0647 observations. In particular they had different grating settings, covering only half of the H$_{2}$ window, as shown in Fig. \ref{supercals}, and they did not have an attached ThAr calibration. Unlike the quasar observations, these exposures were not taken with the slit perpendicular to the horizon and hence an atmospheric dispersion corrector was used not to lose the blue part of the flux due to the atmospheric dispersion. Such differences can affect the light path of Ceres within the telescope, resulting in a different one from the quasar light path. Since the causes of the long-range distortions are unclear \mbox{\citep{Whitmore2015}}, it is not possible to estimate the impact of those differences on the supercalibration process.

Analysis of the Ceres exposures yields two distortion slopes of \mbox{of $\sim 150$ ms$^{-1}$ per 1000 \AA}, while a third one shows a negative distortion \mbox{of $\sim -500$ ms$^{-1}$ per 1000 \AA}. The presence of two  values showing distortions with opposing signs for the same night is remarkable, in particular since these two Ceres spectra were calibrated using the same non-attached ThAr spectrum, and it prevents one to calculate a reliable average value of the distortion correction for these exposures. For these reasons, it was decided to not perform a counter-distortion analysis of the wavelength scale.

Nevertheless, a magnitude of the wavelength distortion can be estimated using the limiting values of the distortion slopes, given in Table \ref{supercaldataset}. Both limiting values of opposing signs were considered as the largest distortions that could have affected the UVES exposures of J1237+0647 in 2009. Therefore the full sub-spectrum formed from 2009 exposures was first counter-distorted using the limit values of the distortion slope and then was combined with the distortion-corrected spectrum formed from exposures taken in 2013 and 2014 to derive the constraints on $\Delta\mu/\mu$. The two values derived were \mbox{$\Delta\mu/\mu\vert_{\textrm{\small{low}}} = (-4.6 \pm 6.1_{\textrm{\small{stat}}}) \times 10^{-6}$} and \mbox{$\Delta\mu/\mu\vert_{\textrm{\small{up}}} = (0.1 \pm 6.1_{\textrm{\small{stat}}}) \times 10^{-6}$} for the lower and for the upper limit respectively. The spread in resulting values for \mbox{$\Delta\mu/\mu$} was interpreted in terms of a systematic uncertainty amounting to $2.4 \times 10^{-6}$.

\subsubsection{Other sources of systematics}
\label{subsec:systs}
The supercalibrations allow to address the issues of both the intra-order and the long-range distortions that are known to affect UVES\citep{Molaro2008,Griest2010,Whitmore2010,Rahmani2013,Whitmore2015}. All the supercalibration targets used in this analysis are presented in Fig. \ref{supercals}. In each panel, the correction that needs to be applied to the ThAr wavelength is shown; note that each echelle order is characterized by several (typically 4-9) measurements separated \mbox{by $\sim$8 \AA}. The vertical spread in each echelle order reflects the magnitude of the intra-order distortions, while the slopes, obtained by least-squares fits to the mean velocity shift correction from each echelle order as a function of wavelength \mbox{\citep[as in][]{Whitmore2015}}, represents the long-range distortions.

The impact of the intra-order distortions is, typically, not expected to be dominant in the systematic error budget on the fiducial value of $\Delta\mu/\mu$, because the H$_{2}$ and HD transitions are spread across multiple orders \mbox{\citep{Malec2010,King2011,Bagdonaite2014}}. Nevertheless, the supercalibrations taken in 2014 show intra-order distortions that are \mbox{$\sim$ three} times larger than the distortions in the exposures taken in 2009 and 2013. An estimation of the impact of the intra-order distortions on the systematic error on $\Delta\mu/\mu$ is made using \mbox{$\delta(\Delta\mu/\mu) = [(\Delta v/c)/\sqrt{N}]/\Delta K_{i}$}, where \mbox{$\Delta v/c$} is the mean amplitude of the intra-order distortions, \mbox{$N=137$} is the number of H$_{2}$ and HD transitions detected and \mbox{$\Delta K_{i}=0.065$} is the spread in the sensitivity coefficients.

The supercalibration exposures return a mean velocity shift of \mbox{$\sim$ 25 ms$^{-1}$} for 2009 and 2013 and \mbox{$\sim$ 60 ms$^{-1}$} for 2014. The shifts translate into systematic uncertainties on $\Delta\mu/\mu$ of \mbox{$1.2\times10^{-6}$} and \mbox{$3.2\times10^{-6}$} respectively. The average value of $\delta(\Delta\mu/\mu)=2.2\times10^{-6}$ is taken as the systematic uncertainty on the $\mu$-variation due to the intra-order distortions. However, this estimation does not take into account various factors like the different sensitivity of the transitions, where they fall with respect to the order centre and the variation of the S/N across the H$_{2}$ window. As a consequence, this value should be considered as an upper limit \mbox{\citep{Malec2010}}.

Long-range and intra-order distortions are not the only contributors to the systematic error budget. In addition, other effects, like the spectral redispersion, the presence of exposures without any ThAr attached calibration and the calibration residuals, contribute to the total systematic error budget. Previous studies addressed these possible error sources, finding that the absence of attached ThAr calibrations introduces an error of  $\sim 0.7 \times 10^{-6}$ \mbox{\citep{Bagdonaite2014}}, the calibration residuals have an effect of, at most, \mbox{$\sim2.0\times10^{-6}$} \mbox{\citep{Murphy2007,Bagdonaite2014}}, and the spectral redispersion can contribute up to \mbox{$\sim 1.4\times10^{-6}$} for a pixel size of \mbox{$\sim 2.5$\,km\,s$^{-1}$} \mbox{\citep{King2011}}. Adding in quadrature all these contributions to the systematic error budget yields a systematic error on the fiducial value of $\Delta\mu/\mu$ of $\sim 4.0 \times 10^{-6}$.
\begin{figure}
\centering
\includegraphics[width=\columnwidth]{supercali.pdf}
\caption{Distortion maps of the UVES wavelength scale derived from asteroid and solar twin targets used for the supercalibration of science exposures of J1237+0647. For each supercalibration exposure, the plot shows the velocity shift required to align the ThAr-calibrated wavelength scale with that derived from the FTS solar spectrum. One observes both intra-order and long-range distortions of the wavelength scale. Only the slopes of the latter, whose values are listed in Table \ref{supercaldataset} with a one-to-one correspondence, are the physically relevant parameters used in the distortion analysis. They were determined by least-squares fits to the mean velocity shift correction from each echelle order as a function of wavelength. The two dashed vertical lines show the window in which the H$_{2}$ absorption lines toward J1237+0647 fall. The radial velocities are unimportant for the analysis of the slopes and have been removed here; the distortion maps have all been shifted to zero velocity shift at 3800 \AA. Top panel: three supercalibrations extracted from Ceres exposures taken in 2009 and retrieved from the ESO archive. Lower panel: seven supercalibrations extracted from asteroid and solar twin observations performed in 2013 and 2014. Eunomia, HD76440 and HD147513 were observed in visitor mode in May 2013 while the solar twins HD097356 and HD117860 were observed in service mode in 2014.
\label{supercals}
\end{figure}

\subsection{Constraining $\Delta\mu/\mu$ from the distortion-corrected spectrum}
\label{subsec:velshiftdm}
After having distortion-corrected the exposures taken in 2013 and 2014, an updated value of \mbox{$\Delta\mu/\mu\vert_{\textrm{\small{2014}}} = (1.2 \pm 11.8_{\textrm{\small{stat}}}) \times 10^{-6}$} was derived from the sub-spectrum formed from those exposures only. This value, which agrees within 1$\sigma$ with the constraint derived only from exposures from 2009, is represented as a square in Fig. \ref{tests}.

The `final spectrum', obtained combining the distortion corrected exposures together with the uncorrected exposures from 2009, delivers a fiducial constraint on the proton-to-electron mass ratio of \mbox{$\Delta\mu/\mu = (-5.4 \pm 6.3_{\textrm{\small{stat}}} \pm 4.0_{\textrm{\small{syst}}}) \times 10^{-6}$}, which is the final result of the present study.
\begin{table*}
\centering
\caption{Observations with VLT/UVES of the supercalibration targets used in this work. Data from program 080.C-0881(B)  were retrieved from the ESO archive and have only the regular ThAr calibration taken at the end of the night. The slit width was \mbox{1.0 arcsec} for all exposures. The uncertainty on the slopes is $\sim$ 30 ms$^{-1}$ per 1000 \AA.}
\label{supercaldataset}
\begin{tabular}{ll*{3}{c}}
\hline
Target & Program ID & Date & Grating [nm] & Distortion slope \\
 & & & & [ms$^{-1}$ per 1000 \AA] \\
\hline
\hline
Ceres & 080.C-0881(B) & 23-03-2009 & 346 & 120 \\
 & & 01-04-2009 & 346 & 200 \\
 & & 01-04-2009 & 346 & $-510$ \\
Eunomia & 091.A-0124(A) & 15-05-2013 & 390+580 & 300 \\
HD76440 & & 15-05-2013 & 390+580 & 400 \\
HD147513 & & 16-05-2013 & 390+580 & 370 \\
HD097356 & 093.A-0373(A) & 23-03-2014 & 390+580 & 650 \\
 & & 03-04-2014 & 390+580 & 720 \\
HD117860 & & 28-05-2014 & 390+580 & 460 \\
 & & 31-05-2014 & 390+580 & 260 \\
\hline
\end{tabular}
\end{table*}

\section{Conclusion}
\label{sec:conclusions}
In this work, the analysis of molecular hydrogen absorption in the system at redshift $z=2.69$, in the line of sight toward quasar J1237+0647, is presented in order to constrain a possible variation of the proton-to-electron mass ratio $\mu$. 137 H$_{2}$ and HD transitions, found in three distinct velocity features associated with the DLA, some even exhibiting further velocity fine-structure, were analyzed. The large number of absorption features considered in the dataset includes partially overlapped features, due to the complex velocity structure of the system, as well as strongly saturated lines. This was possible because of the comprehensive fitting method used in this analysis. Intra-order and long-range distortions, which are known to affect the UVES spectra, were taken into account by applying the supercalibration technique presented by \mbox{\cite{Whitmore2015}} to the exposures taken in 2013 and 2014. For exposures taken in 2009, and retrieved from the ESO archive, the impact of such distortions was estimated using asteroid observations taken within one week of the quasar ones, and it was translated into a contribution to the systematic error budget. The resulting value of a constraint on a varying proton-to-electron mass ratio is \mbox{$\Delta\mu/\mu = (-5.4 \pm 6.3_{\textrm{\small{stat}}} \pm 4.0_{\textrm{\small{syst}}}) \times 10^{-6}$}.

This constraint can be improved using the upcoming new generation of high resolution spectrographs, like the Echelle SPectrograph for Rocky Exoplanet and Stable Spectroscopic Observations \mbox{\citep[ESPRESSO,][]{Pepe2010}}. The spectral resolution achieved in this work allows to resolve components with Doppler broadening parameters \mbox{\emph{b} $\sim2.5$\,km\,s$^{-1}$},  while the high resolution of the new instruments, roughly three times higher than used here, will allow to resolve components with Doppler broadening parameters of the order \mbox{of $\sim$ 1.0\,km\,s$^{-1}$}. This may help further resolve finer structure within individual spectral features and hence result in a more precise $\Delta\mu/\mu$ constraint from \mbox{J1237$+$0647}. Moreover, the new instrument with fiber-fed and frequency-comb support should not suffer from the long-range wavelength distortions  that are affecting UVES, and should lead to smaller systematic uncertainties in $\mu$-variation analyses.

In Fig.~\ref{summary}  the result of the present work is compared with other constraints, taken from literature, on a varying $\mu$ obtained from previous investigations on other H$_2$ absorbers at medium-to-high redshifts. The weighted mean, obtained considering, where possible, both the statistical and the systematic errors, resulting from this larger set is \mbox{$\Delta\mu/\mu=(2.9 \pm 1.7)\times10^{-6}$}. This is consistent with no variation of $\mu$ over a cosmological time-scale at a level of $\sim$10$^{-5}$ for a look-back time of \mbox{$\sim$10.5--12.5 Gyrs}.

\begin{figure}
\centering
\includegraphics[width=1\columnwidth]{summary.pdf}
\caption{Overview of results from investigations on a varying proton-to-electron mass ratio using molecular hydrogen absorbing systems. The result on \mbox{J1237+0647} presented in this work is indicated with a (red) square. The circles show results from previous analyses on systems: a) \mbox{J2123--0050}~\citep{Malec2010,Weerdenburg2011}, b) \mbox{HE0027--1836}~\citep{Rahmani2013}, c) \mbox{Q2348--011}~\citep{Bagdonaite2012}, d) \mbox{Q0405--443}~\citep{King2008}, e) \mbox{B0642--5038}~\citep{Bagdonaite2014}, f) \mbox{Q0528--250}~\citep{King2008,King2011}, g) \mbox{Q0347--383}~\citep{King2008,Wendt2012} and h) \mbox{Q1443+272}~\citep{Bagdonaite2015}. Note that multiple values for the constraint on $\Delta\mu/\mu$ were derived from systems J2123--0050 (labelled as `a'), Q0528--250 (labelled as `f') and Q0347--383 (labelled as `g'). To avoid overlaps between observations of the same system, their points are presented with an offset of +0.05 on the \emph{z}-scale. The dotted line represents the zero level, while the dashed line shows the weighted mean of all the $\Delta\mu/\mu$ values and the shaded area shows its $1\sigma$ boundaries.
}
\label{summary}
\end{figure}•

\section*{Acknowledgments}
The authors thank the Netherlands Foundation for Fundamental Research of Matter (FOM) for financial support. MTM thanks the Australian Research Council for \textsl{Discovery Project} grant DP110100866 which supported this work. The work is based on observations with the ESO Very Large Telescope at Paranal (Chile).

\bibliographystyle{mn2e}


\section{Introduction}
The dynamical relationship between predators and prey has been studied by several authors for a long time. In those researches, various forms of functional responses have been used. Here a functional response means  the average number of prey killed per individual predator per unit of time. Some biologists have argued that in many situation, especially when predators have to search for food, functional responses should depend on both prey's and predator's densities, see \cite{APS,Do,JE,SG} and references therein.

Let us consider a population of three species, say a Lotka-Volterra model,  with the following properties:
\begin{itemize}
\item [\rm (i)] one species is a predator of two competitive other species.
\item [\rm (ii)]  the predator consumes prey with the functional response given by Beddington \cite{Be} and DeAngelis et al. \cite{DGO}
\end{itemize}
 There are many models having the property (i) or (ii) with diffusion in a constant environment \cite{CC1,CC2,CD}, \cite{Linh}, \cite{TVT1,TVT2,TVT3}.
 However, natural environments are usually periodic in time due to the periodicity of seasons. Therefore,  the parameters in these models should be periodic in time. This paper devotes to studying such a
 Lotka-Volterra model which is performed by a nonlinear system of differential equations:
\begin{eqnarray}
\begin{array}{rcl}
x'_1&=&x_1[a_1(t)-b_{11}(t) x_1-b_{12}(t)x_2]\\
&&- \frac{c_1(t)x_1x_3}{\alpha(t)+\beta(t) x_1+\gamma(t) x_3},\\
x'_2&=&x_2[a_2(t)-b_{21}(t) x_1-b_{22}(t)x_2]\\
&&- \frac{c_2(t)x_2x_3}{\alpha(t)+\beta(t) x_2+\gamma(t) x_3},\\
x'_3&=&x_3\Big[-a_3(t)+ \frac{d_1(t)x_1}{\alpha(t)+\beta(t) x_1+\gamma(t) x_3}\\
&&+\frac{d_2(t)x_2}{\alpha(t)+\beta(t) x_2+\gamma(t) x_3}\Big].
\end{array}
\label{E1}
\end{eqnarray}
Here $x_i(t)$ represents the population density of species $X_i$ at time $t \, (i\geqslant 1),$ $X_3$ is a predator species and $ X_1, X_2$ are  competitive prey  species.
At time $t$, $a_i(t)$ is the intrinsic growth rate of  $X_i \, (i=1,2)$ and $a_3(t)$ is the death rate of $X_3$; $b_{ij}(t)$  measures the amount of competition between $X_1$ and $X_2 \, (i\ne j, i,j\leqslant 2),$ and $b_{ii}(t)\, (i\leqslant 2)$ measures the inhibiting effect of environment on $X_i.$  The predator consumes prey with functional responses:
$$ \frac{c_1(t)x_1x_3}{\alpha(t)+\beta(t) x_1+\gamma(t) x_3} \text{ and } \frac{c_2(t)x_2x_3}{\alpha(t)+\beta(t) x_2+\gamma(t) x_3};$$
 and contributes to its growth with amounts:
 $$\frac{d_1(t)x_1}{\alpha(t)+\beta(t) x_1+\gamma(t) x_3} \text{ and } \frac{d_2(t)x_2}{\alpha(t)+\beta(t) x_2+\gamma(t) x_3}.$$
Furthermore, we assume that the parameters $a_i(t), b_{ij}(t), $ $c_i(t),$ $ d_i(t),  \alpha(t), \beta(t), \gamma(t) (1\leqslant i, j\leqslant 3) $ are $\omega$-periodic and continuous in $t$ and bounded below by some positive constants.

 In the next section, we present our main results. First, we use the continuation theorem in coincidence degree theory to show existence of  positive periodic solutions of \eqref{E1}. Second, by using   Lyapunov functions we verify global asymptotic stability  of  boundary periodic solutions. Finally, we give numerical examples.
\section{Main results}
For biological reasons we only consider \eqref{E1} with nonnegative initial values, i.e. $x_1(0), x_2(0),x_3(0)\geqslant 0.$
Let $g(t)$ be a function, for a brevity, instead of writing $g(t)$ we write $g$. If $g$ is a bounded continuous function on $\mathbb R$, we denote\\
\centerline{$g^u= \sup_{t \in \mathbb R}\ g(t), \ g^l=\inf_{t \in \mathbb R} g(t),$}
and
$\hat g=\frac{1}{\omega} \int_0^\omega g(t)dt,$
if $g$ is a periodic function with period $\omega$.
\begin{definition}
A nonnegative solution $x^*(t)$ of \eqref{E1} is called a global asymptotic stable solution  if it attracts any other solution $x(t)$ of \eqref{E1} in the sense that\\
\centerline{$\lim_{t \to \infty}\sum_{i=1}^3|x_i(t)-x_i^*(t)|=0.$}
\end{definition}
\subsection{Existence of  positive periodic solutions}
In this subsection, we shall study existence of  periodic solutions of  \eqref{E1}.  It is not difficult to verify  global existence and uniqueness of nonnegative solutions of  \eqref{E1}.  To show the existence of a positive periodic solution, we shall use the continuation theorem in coincidence degree theory which has been used for some mathematical models of Lotka-Volterra type \cite{L3,TVT2} and references therein. The following are some concepts and results taking from \cite{GM}.

Let $\mathbb X$ and $\mathbb Y$ be two Banach spaces. A linear mapping $L\colon{\mathcal D}(L) \subset \mathbb X \to \mathbb Y$ is called {\it Fredholm} if  it satisfies two conditions:
\begin{itemize}
\item [\rm (i)] $\mathop{\rm Im} L$ is closed and has finite codimension;
\item [\rm (ii)] Ker $L$ has finite dimension.
\end{itemize}
The {\it index} of $L$ is the integer $\mathop{\rm dim} \mathop{\rm Ker} L -\mathop{\rm codim}\mathop{\rm Im} L$.
If $L$ is Fredholm  of index zero, there exist continuous projections $ P\colon  \mathbb X \to \mathbb X$ and $ Q\colon \mathbb Y \to \mathbb Y$ such that $\mathop{\rm Im} P=\mathop{\rm Ker} L, \mathop{\rm Im} L=\mathop{\rm Ker} Q=\mathop{\rm Im} (I-Q),$ and an isomorphism $J\colon  \mathop{\rm Im} Q \to \mathop{\rm Ker} L$.
It follows that\\
\centerline{$L_p=L|_{{\mathcal D}(L) \cap \mathop{\rm Ker} P}\colon  (I-P) {\mathbb X }\to \mathop{\rm Im} L$}
is invertible. We denote the inverse of that map
by $K_p$. Let $\Omega$ be an open bounded subset of $\mathbb X$. A continuous mapping $N\colon\mathbb X \to\mathbb Y$ is said to be $L$-{\it compact} on $\bar\Omega$ if the following two conditions take place:
\begin{itemize}
\item [\rm (i)] the mapping $QN\colon  \bar\Omega \to \mathbb Y$ is continuous and bounded;
\item [\rm (ii)] $K_p(I-Q)N\colon\bar \Omega\to \mathbb X$ is compact, i.e. it is continuous and   $K_p(I-Q)N(\bar \Omega)$ is relatively compact.
\end{itemize}
 To introduce the definition of the degree of  $N$ in $\Omega$, for simplicity we assume that ${\mathbb X}={\mathbb R}^N.$  Suppose furthermore that $N$ is smooth on $\bar \Omega$. Let $p\notin \partial \Omega$ be a regular value of $N$, i.e. the equation $N(x)=p$ on $\bar \Omega$ has only a finite number of solutions $x_1,\dots,x_n\in \Omega$ with  nonsingular $DN(x_i)$ for each $i=1,\dots,n$ where $DN(x_i)$ is the Jacobi matrix of $N$ at $x_i$. Then the degree $\deg(N,\Omega,p) $ of $N$ in $\Omega$ at $p$ is defined by the formula
$$
\deg(N,\Omega,p)=\sum_{i=1}^n \mathop{\rm sgn}\{\det DN(x_i)\}.
$$
\begin{lemma} [Continuation theorem \cite{GM}] \label{lem3}
Let L be a Fredholm mapping of index $0$. Assume that $N\colon  \bar \Omega \to \mathbb Y$ is $L$-compact
on $\bar \Omega$ and satisfies conditions:
\begin{itemize}
\item [\rm (a)] for each $\lambda \in (0, 1)$ every solution of $Lx = \lambda Nx$ is such that $x \notin \partial \Omega$,
\item [\rm (b)] $QN x 	\ne {\bf 0}$ for each $x \in \partial\Omega \cap \mathop{\rm Ker} L$, and
$\deg \{JQN, \Omega \cap \mathop{\rm Ker} L, {\bf 0}\} \ne 0.$
\end{itemize}
Then the operator equation $Lx = Nx$ has at least one solution in ${\mathcal D}(L) \cap $ $\bar\Omega.$
\end{lemma}

We now put
\begin{equation*}
\begin{aligned}
&L_{i1}=\ln\frac{\hat a_i}{\hat b_{ii}}, H_{i1}=L_{i1}+ 2 \hat a_i \omega \quad (i=1,2),\\
&L_{12}=\ln \Big\{\frac{\hat a_1-\hat b_{12} e^{H_{21}}-\widehat{(\frac{c_1}{\gamma })}}{\hat b_{11}}\Big\}, \\
&H_{12}=L_{12}-2 \hat a_1 \omega,\\
&L_{22}=\ln \Big\{\frac{\hat a_2-\hat b_{21} e^{H_{11}}-\widehat{(\frac{c_2}{\gamma })}}{\hat b_{22}}\Big\}, \\
&H_{22}=L_{22}-2 \hat a_2 \omega,\\
&L_{31}=\ln\Big\{\frac{\hat d_1e^{H_{11}}+\hat d_2e^{H_{21}}-\hat a_3 \alpha^l}{\hat a_3 \gamma ^l}\Big\}, \\
&H_{31}=L_{31}+2 \hat a_3 \omega,\\
&L_{32}=\ln[(\hat d_1-\hat a_3 \beta^u)e^{H_{12}}\\
&\hspace{1.2cm}+(\hat d_2-\hat a_3 \beta^u)e^{H_{22}}-2\hat a_3 \alpha^u]-\ln(2\hat a_3 \gamma ^u),\\
& H_{32}=L_{32}-2 \hat a_3 \omega.
\end{aligned}
\end{equation*}
The convention here is that $\ln x = -\infty$ if $x\leqslant  0$. Under the  conditions
\begin{equation} \label{H2}
\begin{cases}
\hat b_{11}\hat b_{22}\ne \hat b_{12}\hat b_{21},\\
\hat a_1-\hat b_{12} e^{H_{21}}-\widehat{(\frac{c_1}{\gamma })}>0,\\
\hat a_2-\hat b_{21} e^{H_{11}}-\widehat{(\frac{c_2}{\gamma })}>0,\\
\hat d_1e^{H_{11}}+\hat d_2e^{H_{21}}-\hat a_3 \alpha^l>0,\\
(\hat d_1-\hat a_3 \beta^u)e^{H_{12}}+(\hat d_2-\hat a_3 \beta^u)e^{H_{22}}>2\hat a_3 \alpha^u,
\end{cases}
\end{equation}
we shall verify existence of an $\omega$-periodic solution of \eqref{E1}.
\begin{theorem} \label{thm:5}
Let   \eqref{H2} be satisfied. Then \eqref{E1} has at least one positive $\omega$-periodic solution.
\end{theorem}
\begin{proof}
  By putting
$x_i(t)=e^{u_i(t)} (i\geqslant 1)$,  \eqref{E1} becomes
\begin{equation}\label{E14}
\begin{cases}
u'_1=a_1-b_{11}e^{u_1} -b_{12}e^{u_2}- \frac{c_1e^{u_3}}{\alpha+\beta e^{u_1}+\gamma e^{u_3}},\\
u'_2=a_2-b_{21}e^{u_1} -b_{22}e^{u_2}- \frac{c_2e^{u_3}}{\alpha+\beta e^{u_2}+\gamma e^{u_3}},\\
u'_3=-a_3+ \frac{d_1e^{u_1}}{\alpha+\beta e^{u_1}+\gamma e^{u_3}}+\frac{d_2e^{u_2}}{\alpha+\beta e^{u_2}+\gamma e^{u_3}}.
\end{cases}
\end{equation}
Let
\begin{equation*}
\begin{aligned}
\mathbb X&=\mathbb Y\\
&=\{u=(u_1, u_2, u_3)^T\in  C^1(\mathbb R,\mathbb R^{\text{$3$}})\,{\text{such that}}\\
& \hspace{1cm}u_i(s)=u_i(s+\omega) \text{ for } s\in \mathbb R \text{ and   $i$} \geqslant \text{$1$}\},
\end{aligned}
\end{equation*}
with norm
$$ ||u||=\sum_{i=1}^3 \text{\max}_{s\in [0,\omega]} |u_i(s)|, \hspace{1cm} u\in \mathbb X.$$
Then both $\mathbb X$ and $ \mathbb Y$ are  Banach spaces.
Let
\begin{equation*}
\begin{aligned}
&N
\begin{bmatrix}
u_1\\
u_2\\
u_3
\end{bmatrix} (s)
=
\begin{bmatrix}
N_1(s)\\
N_2(s)\\
N_3(s)
\end{bmatrix}\\
&
=
\begin{bmatrix}
a_1-b_{11}e^{u_1} -b_{12}e^{u_2}- \frac{c_1e^{u_3}}{\alpha+\beta e^{u_1}+\gamma e^{u_3}}\\
a_2-b_{21}e^{u_1} -b_{22}e^{u_2}- \frac{c_2e^{u_3}}{\alpha+\beta e^{u_2}+\gamma e^{u_3}}\\
-a_3+ \frac{d_1e^{u_1}}{\alpha+\beta e^{u_1}+\gamma e^{u_3}}+\frac{d_2e^{u_2}}{\alpha+\beta e^{u_2}+\gamma e^{u_3}}
\end{bmatrix},\\
&L
\begin{bmatrix}
u_1\\
u_2\\
u_3
\end{bmatrix}
=
\begin{bmatrix}
u_1'\\
u_2'\\
u_3'
\end{bmatrix},\\
&P
\begin{bmatrix}
u_1\\
u_2\\
u_3
\end{bmatrix}=
Q
\begin{bmatrix}
u_1\\
u_2\\
u_3
\end{bmatrix}=
\begin{bmatrix}
\frac{1}{\omega} \int_0^\omega u_1(s) ds\\
\frac{1}{\omega} \int_0^\omega u_2(s) ds\\
\frac{1}{\omega} \int_0^\omega u_3(s) ds
\end{bmatrix}.
\end{aligned}
\end{equation*}
Hence,
$$\mathop{\rm Ker} L= {\mathbb R^{\text 3}}, \mathop{\rm Im} L=\{u\in \mathbb Y\,{|}\, \text{$\int_0^\omega u_i(s)ds=0, \, i\geq1\}$},$$
and $\mathop{\rm dim} \mathop{\rm Ker} L=3=\mathop{\rm codim} \mathop{\rm Im} L.$
Then, it is easy to obtain the following conclusions.
\begin{enumerate}
  \item  $L$ is a Fredholm mapping of index zero, since $\mathop{\rm Im} L$ is closed in $\mathbb Y$.
  \item $P$ and $ Q$ are continuous projections such that $\mathop{\rm Im} P=\mathop{\rm Ker} L, \mathop{\rm Im} L=\mathop{\rm Ker} Q=\mathop{\rm Im} (I-Q)$.
  \item The generalized inverse (to $L$) $K_P\colon\mathop{\rm Im} L \to {\mathcal D}(L)\cap \mathop{\rm Ker}  P$ exists and is given by
\begin{equation*}
K_P
\begin{bmatrix}
u_1\\
u_2\\
u_3
\end{bmatrix} (\nu)
=
\begin{bmatrix}
\int_0^\nu u_1(s)ds-\frac{1}{\omega}\int_0^\omega \int_0^\nu u_1(s)dsd\nu\\
\int_0^\nu u_2(s)ds-\frac{1}{\omega}\int_0^\omega \int_0^\nu u_2(s)dsd\nu\\
\int_0^\nu u_3(s)ds-\frac{1}{\omega}\int_0^\omega \int_0^\nu u_3(s)dsd\nu\\
\end{bmatrix}.
\end{equation*}
\item
 $QN$ and $K_P(I-Q)N$ are  continuous.
\item  $N$ is $L$-compact on $\bar \Omega$ with any open bounded set $\Omega\subset \mathbb X.$
\end{enumerate}

Let us now find an appropriate open, bounded subset $\Omega$ for application of the continuation theorem.  Obviously, from   \eqref{H2},  $-\infty< L_{i2}\leqslant L_{i1}<\infty \, (i\geqslant 1).$ Corresponding to the  equation $Lu=\lambda N u, \lambda\in (0,1),$ we have
\begin{align}
&u_1'=\lambda\Big[a_1-b_{11}e^{u_1} -b_{12}e^{u_2}- \frac{c_1e^{u_3}}{\alpha+\beta e^{u_1}+\gamma e^{u_3}}\Big], \notag\\
&u_2'=\lambda\Big[a_2-b_{21}e^{u_1} -b_{22}e^{u_2}- \frac{c_2e^{u_3}}{\alpha+\beta e^{u_2}+\gamma e^{u_3}}\Big],   \label{E15}\\
&u_3'=\lambda\Big[-a_3+ \frac{d_1e^{u_1}}{\alpha+\beta e^{u_1}+\gamma e^{u_3}}\notag\\
&\hspace{1.2cm}+\frac{d_2e^{u_2}}{\alpha+\beta e^{u_2}+\gamma e^{u_3}}\Big].\notag
\end{align}
Suppose that $(u_1, u_2, u_3)^T\in \mathbb X$ is an arbitrary solution of \eqref{E15}. Integrating both the hand sides of \eqref{E15} over the interval $[0,\omega]$, we obtain
\begin{align}
&\hat a_1 \omega=\int_0^\omega\Big[b_{11}e^{u_1}+b_{12}e^{u_2}+ \frac{c_1e^{u_3}}{\alpha+\beta e^{u_1}+\gamma e^{u_3}}\Big]dt,  \notag\\
&\hat a_2 \omega=\int_0^\omega\Big[b_{21}e^{u_1}+b_{22}e^{u_2}+ \frac{c_2e^{u_3}}{\alpha+\beta e^{u_2}+\gamma e^{u_3}}\Big]dt,  \label{E16}\\
&\hat a_3 \omega=\int_0^\omega\Big[\frac{d_1e^{u_1}}{\alpha+\beta e^{u_1}+\gamma e^{u_3}}+\frac{d_2e^{u_2}}{\alpha+\beta e^{u_2}+\gamma e^{u_3}}\Big]dt. \notag
\end{align}
Combining the first equations of \eqref{E15} and \eqref{E16}, we observe that
\begin{equation*} \label{E17}
\begin{aligned}
&\int_0^\omega|u_1'|dt \\
&\leqslant \lambda\Big[\int_0^\omega a_1 dt+\int_0^\omega b_{11}e^{u_1}dt\\
 &\hspace{1cm} +\int_0^\omega b_{12}e^{u_2}dt+ \int_0^\omega \frac{c_1e^{u_3}}{\alpha+\beta e^{u_1}+\gamma e^{u_3}}\Big]dt\\
 &<  2 \hat a_1 \omega.
 \end{aligned}
\end{equation*}
Similarly, we have
$\int_0^\omega |u_2'(t)|dt < 2 \hat a_2 \omega,$
and
\begin{align*}
&\int_0^\omega|u_3'(t)|dt \leqslant   \lambda\Big[\int_0^\omega a_3 dt+ \int_0^\omega \frac{d_1e^{u_1}}{\alpha+\beta e^{u_1}+\gamma e^{u_3}}dt\\
 &+ \int_0^\omega\frac{d_2e^{u_2}}{\alpha+\beta e^{u_2}+\gamma e^{u_3}}\Big]dt< 2 \hat a_3 \omega.
\end{align*}
Since $u \in \mathbb X$, there exist $\xi_i, \eta_i\in [0, \omega] \, (i\geqslant 1)$ such that
\begin{equation}\label{E18}
u_i(\xi_i)=\text{\min}_{t\in [0,\omega]} u_i(t), \quad u_i(\eta_i)=\text{\max}_{t\in [0,\omega]} u_i(t).
\end{equation}
From the first equation of \eqref{E16} and \eqref{E18}, we obtain
\begin{equation*}
\begin{aligned}
\hat a_1 \omega &\geqslant \int_0^\omega b_{11}e^{u_1(\xi_1)}dt+\int_0^\omega b_{12}e^{u_2(\xi_2)}dt\\
&=\hat b_{11} \omega e^{u_1(\xi_1)}+\hat b_{12} \omega e^{u_2(\xi_2)},
\end{aligned}
\end{equation*}
which implies that $u_1(\xi_1)< L_{11}.$ Hence,  for all $ t\geqslant 0$
$$u_1(t)\leqslant u_1(\xi_1)+\int_0^\omega|u_1'(t)|dt < L_{11}+ 2 \hat a_1 \omega=H_{11}.$$
Similarly, we have $u_2(t) < H_{21}$ for all $t\geqslant 0$.\\
 On the other hand, from the first equation of \eqref{E16} and \eqref{E18},
\begin{equation*}
\begin{aligned}
&\hat a_1 \omega \\
&\leqslant \int_0^\omega b_{11}e^{u_1(\eta_1)}dt+\int_0^\omega b_{12}e^{u_2(\eta_2)}dt+\int_0^\omega \frac{c_1(t)}{\gamma (t)}dt\\
&=\hat b_{11} \omega e^{u_1(\eta_1)}+\hat b_{12} \omega e^{u_2(\eta_2)}+\widehat{(\frac{c_1}{\gamma })}\omega\\
&\leqslant \hat b_{11} \omega e^{u_1(\eta_1)}+\hat b_{12} \omega e^{H_{21}}+\widehat{(\frac{c_1}{\gamma })}\omega.
\end{aligned}
\end{equation*}
Hence,
$$u_1(t)\geqslant u_1(\eta_1)-\int_0^\omega |u_1'(t)|dt \geqslant H_{12},  \hspace{1cm} \forall  t\geqslant 0.$$
Similarly,  $u_2(t) \geqslant H_{22}$ for all $t\geqslant 0.$ Therefore, by putting
$B_i=\text{\max} \{|H_{i1}|, |H_{i2}|\}$, we conclude that
$$\text{\max}_{t\in [0, \omega]} |u_i(t)| \leqslant B_i, \hspace{1cm} i=1,2.$$

Let us give estimates for $u_3(t)$. It follows from the third equation of \eqref{E16} and \eqref{E18} that
\begin{equation*}
\begin{aligned}
\hat a_3 \omega &\leqslant \int_0^\omega\Big[\frac{d_1(t)e^{H_{11}}}{\alpha^l+\gamma ^l e^{u_3(\xi_3)}}+\frac{d_2(t)e^{H_{21}}}{\alpha^l+\gamma ^l e^{u_3(\xi_3)}}\Big]dt\\
&=\frac{[\hat d_1e^{H_{11}}+\hat d_2e^{H_{21}}]\omega}{\alpha^l+\gamma ^l e^{u_3(\xi_3)}}
\end{aligned}
\end{equation*}
and
\begin{equation*}
\begin{aligned}
\hat a_3 \omega \geqslant &\int_0^\omega\Big[\frac{d_1(t)e^{H_{12}}}{\alpha^u+\beta^u e^{H_{12}}+ \gamma ^u e^{u_3(\eta_3)}}\\
&+\frac{d_2(t)e^{H_{22}}}{\alpha^u+\beta^u e^{H_{22}}+\gamma ^u e^{u_3(\eta_3)}}\Big]dt\\
= &\frac{\hat d_1e^{H_{12}}\omega}{\alpha^u+\beta^u e^{H_{12}}+ \gamma ^u e^{u_3(\eta_3)}}\\
&+\frac{\hat d_2e^{H_{22}}\omega}{\alpha^u+\beta^u e^{H_{22}}+\gamma ^u e^{u_3(\eta_3)}}\\
\geqslant &\frac{[\hat d_1e^{H_{12}}+\hat d_2e^{H_{22}}]\omega}{2\alpha^u+\beta^u [e^{H_{12}}+e^{H_{22}}]+2\gamma ^u e^{u_3(\eta_3)}}.
\end{aligned}
\end{equation*}
Hence, $u_3(\xi_3) \leqslant L_{31}$ and $u_3(\eta_3) \geqslant L_{32}.$ We then observe that \\
\centerline{$
u_3(t) \leqslant u_3(\xi_3)+ \int_0^\omega|u_3'(t)| dt \leqslant H_{31}
$}
and
$$
u_3(t) \geqslant u_3(\eta_3)- \int_0^\omega|u_3'(t)| dt \geqslant H_{32}.$$
Therefore, by putting $B_3=\text{\max} \{|H_{31}|, |H_{32}|\},$ we get
$$\text{\max}_{t\in [0, \omega]}|u_3| \leqslant B_3.$$

By the above estimates, for any solution $u\in \mathbb X$ of \eqref{E16} we have
$||u|| \leqslant \sum_{i=1}^3 B_i.$
Clearly, $B_i \, (i\geqslant 1)$ are independent of $\lambda.$ Take $B=\sum_{i=1}^4 B_i$ where $B_4$ is taken sufficiently large such that
\text{$B_4\geqslant \sum_{i=1}^3\sum_{j=1}^2 |L_{ij}|.$}
Let $\Omega=\{u\in \mathbb X \,{|}\, \text{$ ||u||<B$}\},$
then $\Omega$ satisfies the condition (a) of Lemma \ref{lem3}.

Let us verify that the condition (b) of Lemma \ref{lem3} is also satisfied. Consider the homotopy
$$H_{\mu}(u)=\mu QN(u)+(1-\mu)G(u), \hspace{1cm} \mu\in[0,1]$$
where $G\colon \mathbb R^{\text 3} \to \mathbb R^{\text 3}$,
$$G(u)=
\begin{bmatrix}
\hat a_1-\hat b_{11} e^{u_1}-\hat b_{12} e^{u_2}\\
\hat a_2-\hat b_{21} e^{u_1}-\hat b_{22} e^{u_2}\\
-\hat a_3+f(u_1,u_3)+g(u_2,u_3)
\end{bmatrix}$$
with
$f(u_1, u_3)=\frac{1}{\omega} \int_0^\omega \frac{d_1e^{u_3}dt}{\alpha+\beta e^{u_1}+\gamma e^{u_3}}$ and
 $g(u_2,u_3)= \frac{1}{\omega}\int_0^\omega \frac{d_2e^{u_3}dt}{\alpha+\beta e^{u_2}+\gamma e^{u_3}}.$
It is easy to see that
\begin{align*}
&H_{\mu}(u)\\
&=
\begin{bmatrix}
\hat a_1-\hat b_{11} e^{u_1}-\hat b_{12} e^{u_2}-\frac{1}{\omega}\int_0^\omega \frac{\mu c_1e^{u_3} dt}{\alpha+\beta e^{u_1}+\gamma e^{u_3}}\\
\hat a_2-\hat b_{21} e^{u_1}-\hat b_{22} e^{u_2} -\frac{1}{\omega}\int_0^\omega \frac{\mu c_2e^{u_3}dt}{\alpha+\beta e^{u_2}+\gamma e^{u_3}}\\
-\hat a_3+f(u_1,u_3)+g(u_2,u_3)
\end{bmatrix}.
\end{align*}
By carrying out similar arguments as above, we observe  that any solution $u^*$ of the equation $H_{\mu}(u)={\bf 0} \in \mathbb R^{\text 3}$ with $\mu\in [0,1]$ satisfies the estimate
\begin{equation} \label{E22}
L_{i2}\leqslant u_i^*\leqslant L_{i1}, \hspace{1cm} i\geqslant 1.
\end{equation}
Thus, ${\bf 0}\notin H_{\mu}( \partial \Omega \cap \mathop{\rm Ker} L)$ for $\mu\in [0,1]$. Consequently, by taking $\mu=1$,  we conclude that ${\bf 0}\notin QN(\partial\Omega \cap \mathop{\rm Ker} L)$. Note that the isomorphism $J$ can be the identity mapping $I$, since $\mathop{\rm Im} P=\mathop{\rm Ker} L.$ By the invariance property of homotopy, we obtain that
\begin{equation*} \label{E23}
\begin{aligned}
&\deg (JQN,  \Omega  \cap  \mathop{\rm Ker} L, {\bf 0})=\deg(QN,  \Omega \cap \mathop{\rm Ker} L, {\bf 0})\\
&=\deg(QN,  \Omega \cap {\mathbb R^{\text 3}}, {\bf 0})=\deg(G,  \text{$\Omega$} \cap \mathbb R^{\text 3}, {\bf 0})\\
&=\mathop{\rm sgn}\det \Lambda\\
&=\mathop{\rm sgn}\left\{ (\hat b_{11}\hat b_{22}-\hat b_{12}\hat b_{21})  \frac{\partial [f(u_1,u_3)+g(u_2,u_3)]}{\partial u_3}\right\}
\end{aligned}
\end{equation*}
where
\begin{align*}
&\Lambda=
\begin{bmatrix}
-\hat b_{11} e^{u_1}& -\hat b_{12} e^{u_2} &0\\
-\hat b_{21} e^{u_1}& -\hat b_{22} e^{u_2} &0\\
\frac{\partial f(u_1,u_3)}{\partial u_3} & \frac{\partial g(u_2,u_3)}{\partial u_3} & \frac{\partial f(u_1,u_3)}{\partial u_3}+ \frac{\partial g(u_2,u_3)}{\partial u_3}
\end{bmatrix}.
\end{align*}
Since  both functions $f(u_1, u_3)$ and $ g(u_2,u_3)$  increase in $u_3$,  $
\frac{\partial f(u_1,u_3)}{\partial u_3}+ \frac{\partial g(u_2,u_3)}{\partial u_3}>0.
$ Hence, by using the first condition in \eqref{H2}, we conclude that \\
\centerline{$\deg (JQN,  \Omega  \cap  \mathop{\rm Ker} L, {\bf 0})\ne 0.$}

By now we have proved that $\Omega$ verifies all requirements of Lemma \ref{lem3}. Therefore, the equation  $Lu=Nu$ has at least one solution in ${\mathcal D}(L)\cap \bar \Omega,$ i.e. \eqref{E14} has at least one $\omega$-periodic solution $u^*$ in ${\mathcal D}(L)\cap \bar \Omega.$ Set $x_i^*=e^{u_i^*} (i\geqslant 1),$ then $x^*$ is an $\omega$-periodic solution of \eqref{E1} with strictly positive components. It completes the proof.
\end{proof}
\subsection{Global asymptotic  stability  of  boundary periodic solutions}
In this subsection, we shall establish a sufficient criteria for   global asymptotic  stability of boundary $\omega$-periodic solutions  of \eqref{E1}.  Consider  the boundary dynamics of \eqref{E1} where $X_3$ is absent, i.e. $x_3(t)=0$ for every $t\geqslant 0$. We then consider the  periodic competitive model of two prey $X_1, X_2$:
\begin{equation}\label{E25}
\begin{cases}
x'_1=x_1\left[a_1(t)-b_{11}(t) x_1-b_{12}(t)x_2\right],\\
x'_2=x_2\left[a_2(t)-b_{21}(t) x_1-b_{22}(t)x_2\right].
\end{cases}
\end{equation}
Denote by $\bar X_i(t)$  the unique positive $\omega$-periodic solution  of the logistic equation:\\
\centerline{$X'=X\left[a_i(t)-b_{ii}(t) X\right].$}
Then
$\bar X_i(t)=\frac{e^{\int_0^\omega a_i(s)ds -1}}{\int_t^{t+\omega} b_{ii}(s) e^{-\int_s^t a(\tau)d\tau}ds}.$
Due to \cite{Lb}, if
\begin{equation} \label{E26}
\hat {a_i} > \widehat {b_{ij} \bar X_j}\hspace{1cm}  (i\ne j, i, j=1,2),
\end{equation}
then  \eqref{E25}  has a positive $\omega$-periodic solution $(\bar x_1, \bar x_2)$. Furthermore,
if
\begin{equation} \label{E27}
\hat A_{12}<0
\end{equation}
then $(\bar x_1, \bar x_2)$ is globally asymptotically stable, where
$a_{ij}(t)=b_{ij}(t) \bar x_j(t) \,\, ( i\ne j, i,j=1,2)$ and
$$A_{12}(t)=\text{\max}\left\{ \frac{(a_{ij}+a_{ji})^2}{4 a_{ii}}-a_{jj}, i\ne j\right\}.$$
 Our result is as follows.
\begin{theorem}
If \eqref{E26} and \eqref{E27} hold then $\bar x=(\bar x_1, \bar x_2, 0)$ is a  $\omega$-periodic boundary solution of \eqref{E1}. Furthermore,
\begin{itemize}
\item [{\rm (i)}]
If $b_{ij}<b_{jj} \, (1\leqslant i\ne j\leqslant 2)$ and  $ c_1+c_2+d_1+d_2< \beta a_3$ then $\bar x$ is globally asymptotically stable.
\item[{\rm (ii)}]  If $ c_1+c_2+d_1+d_2< \beta a_3$ then $\bar x$ attracts any solution $x$ of \eqref{E1} which satisfies the condition
 $$[x_1(t)-\bar x_1(t)] [x_2(t)-\bar x_2(t)] \geqslant 0, \hspace{1cm} \forall t\geqslant 0.$$
\item[{\rm (iii)}] If $ d_1+d_2< \beta a_3$ then $\bar x$ attracts any solution $x$ of \eqref{E1} which satisfies the condition
 $$x_i(t) \geqslant \bar x_i(t), \hspace{1cm} \forall t\geqslant 0, i=1, 2.$$
\end{itemize}
\end{theorem}
\begin{proof}
The first statement is obvious. To prove (i), let  $x$ be any other solution of \eqref{E1}.
Consider a Lyapunov function defined by
 $V(t)=\sum_{i=1}^2 |\ln x_i-\ln \bar x_i|+ x_3,  t \geqslant 0.$
Calculating of the right derivative $D^{+}V(t)$ of  $V(t)$ along the solutions of \eqref{E1} gives
\begin{align} \label{E28}
&D^+V(t) \notag\\
= &\sum_{i=1}^2 \mathop{\rm sgn}(x_i-\bar x_i) (\frac{x_i'}{x_i}-\frac{\bar x_i'}{\bar x_i}) +x_3' \notag\\
= &\sum_{i\ne  j}^2 \Big\{\mathop{\rm sgn}(x_i-\bar x_i) \Big[(a_i-b_{ii}x_i-b_{ij}x_j\notag\\
&-\frac{c_i x_i x_3}{\alpha+\beta x_i+\gamma x_3})-(a_i-b_{ii} \bar x_i-b_{ij} \bar x_j)\Big] \Big\}\notag\\
&+(-a_3+\sum_{i=1}^2 \frac{d_i x_i}{\alpha+\beta x_i+\gamma x_3}) x_3\notag\\
=&\sum_{i\ne  j}^2 \left\{[ -b_{ii}|x_i-\bar x_i|-b_{ij} (x_j-\bar x_j) \mathop{\rm sgn}(x_i-\bar x_i)  ]\right\}\notag\\
&+ x_3 \Big [-a_3+\sum_{i=1}^2 \Big\{\frac{d_i x_i}{\alpha+\beta x_i+\gamma x_3}\notag\\
&\hspace{3cm}-\frac{c_i x_i \,\mathop{\rm sgn}(x_i-\bar x_i)}{\alpha+\beta x_i+\gamma x_3}\Big\}\Big].
\end{align}
Then
\begin{align}
D^+V(t) \leqslant & \sum_{i\ne  j}^2 (b_{ij}-b_{jj}) |x_j-\bar x_j| \notag\\
&+\frac{(c_1+c_2+d_1+d_2-\beta a_3) x_3}{\beta}. \label{E13}
\end{align}
By assumptions in (i) and the periodicity of parameters, there exist $\mu_1>0$ such that
\begin{align*}
\text{\max}_{t\in [0,\omega], 1\leqslant i\ne j\leqslant 2} &\left\{ \frac{c_1+c_2+d_1+d_2-\beta a_3}{\beta}, b_{ij}-b_{jj}\right\}\\
&<-\mu_1.
\end{align*}
Thus, by integrating both the hand sides of \eqref{E13} from  $0$ to $t$, we observe that\\
\centerline{$V(t)+\mu_1 \int_0^{t}\sum_{i=1}^3 |x_i-\bar x_i|ds \leqslant V(0) <\infty$}
for every $t \geqslant 0$.
Hence,
$ \sum_{i=1}^3 |x_i-\bar x_i| \in L^1([0, \infty)).$

On the other hand, by the periodicity,  $x_i$ and  $\bar x_i \, (i\geqslant 1)$ have bounded derivatives on $[0,\infty)$. As a consequence,
$ \sum_{i=1}^3 |x_i-\bar x_i|$ is uniformly continuous on $[0, \infty)$.
Therefore, by using the Barbalat lemma \cite{TVT2}, we conclude that\\
 \centerline{$\underset{t \to \infty}{\lim}\sum_{i=1}^3 |x_i-\bar x_i|=0,$}
i.e.  $\bar x$ is globally asymptotically stable.

Similarly,  we obtain the conclusions in (ii) and (iii) by using the following inequalities, respectively.
\begin{align*}
&D^+V(t)\\
 =&\sum_{i\ne  j}^2 [ -(b_{ii}+b_{ji}) |x_i-\bar x_i|]+ x_3 \Big [-a_3\notag\\
&+\sum_{i=1}^2 \Big\{\frac{d_i x_i}{\alpha+\beta x_i+\gamma x_3}-\frac{c_i x_i \,\mathop{\rm sgn}(x_i-\bar x_i)}{\alpha+\beta x_i+\gamma x_3}\Big\}\Big]\\
\leqslant & \sum_{i\ne  j}^2 [ -(b_{ii}+b_{ji}) |x_i-\bar x_i|]\notag\\
&+\frac{(c_1+c_2+d_1+d_2-\beta a_3) x_3}{\beta},
\end{align*}
and
\begin{align*}
&D^+V(t) \\
=&\sum_{i\ne  j}^2 [ -(b_{ii}+b_{ji}) |x_i-\bar x_i|]+ x_3 \Big [-a_3\notag\\
&+\sum_{i=1}^2 \Big\{\frac{d_i x_i}{\alpha+\beta x_i+\gamma x_3}-\frac{c_i x_i }{\alpha+\beta x_i+\gamma x_3}\Big\}\Big]\\
\leqslant & \sum_{i\ne  j}^2 \left[ -(b_{ii}+b_{ji}) |x_i-\bar x_i|\right]+\frac{(d_1+d_2-\beta a_3) x_3}{\beta}.
\end{align*}
We complete the proof.
\end{proof}
\subsection{Numerical examples}
In this subsection, we exhibit some numerical examples which show the convergence of  positive solutions of  \eqref{E1} to  periodic solutions of  \eqref{E1}.
Set $a_1=3+\sin(8t); a_2=5.5-0.2\cos(8t); a_3=0.4-0.3\cos(8t); b_{11}=2+\cos(8t); b_{22}=5+0.4\sin(8t); b_{12}=0.04-0.02\sin(8t); b_{21}=0.15-0.1\cos(8t); c_1=0.5-0.4\sin(8t);$ $ c_2=0.4-0.3\sin(8t); \alpha=0.03-0.02\cos(8t);
\beta=0.3+0.2\cos(8t);
\gamma=2-\sin (8t);
d_1=3+2 \sin(8t);
 d_2=3-2 \sin(8t);
$ and an initial value $(x_1(0),x_2(0),x_3(0))=(0.5,0.7,1).$ Figure \ref{Fig1} shows the behavior of the solution of \eqref{E1}. It is seen that the solution converges to a positive periodic solution of \eqref{E1}.
 \begin{figure}[h]
\begin{center}
\includegraphics[scale=0.7]{PositivePeriodicSol.eps}
\end{center}
\vspace*{9pt}
\caption{ A solution of  \eqref{E1} which converges to a positive periodic solution of \eqref{E1}}
\label{Fig1}
\end{figure}

We now set $a_3=4-0.3\cos(8t)$ and $\beta=3+0.2\cos(8t)$ and retain other parameters as above. Figure \ref{Fig2} gives the behavior of the positive solution of \eqref{E1}. It converges to the boundary periodic solution of \eqref{E1}.

 \begin{figure}[h]
\begin{center}
\includegraphics[scale=0.7]{BoundaryPeriodicSol.eps}
\end{center}
\vspace*{9pt}
\caption{A  positive solution of  \eqref{E1} which converges to the boundary periodic solution of \eqref{E1}}
\label{Fig2}
\end{figure}
\section*{Acknowledgement}
 The work of the last author is supported by JSPS KAKENHI Grant Number 20140047. The authors would like to thank the anonymous referees for their helpful suggestions which improved the
paper.