\section{Introduction}
Since their initial demonstration by Maiman~\cite{maiman1960stimulated}, lasers have gradually become perhaps the most important tool for neutral gas, plasma and nanoparticle diagnostics, both stagnant and in a flow. Exhibiting the significant advantages of performing fast, non-intrusive and non-perturbative measurements, the ultimate goal of laser diagnostics is to replace the use of mechanical probes that are still the state of the art for relevant measurements, such as Pitot tubes and hot-wires for neutral gas flows~\cite{shinder2015nist,4198449,RUSSO201167,gatski2013compressibility,morris2012measurement,JENSEN2004}, Langmuir and emissive probes in the case of plasma flows~\cite{10.1063/1.5016229,cherrington1982use} and mechanical extractors in the case of nanoparticles~\cite{FANG2016273}.

Mechanical probes are inherently intrusive to the system they measure, and they exhibit generally slow response times while sometimes they do not survive the environment which they are called to measure. To overcome these limitations, a plethora of laser based diagnostic techniques have been developed and demonstrated over the years to help measure thermodynamic quantities of interest in all states of matter. Those designed to perform measurements on neutral gases include, to name but a few, Spontaneous Rayleigh Scattering (SRayS)~\cite{Miles_2001,Kampmann:93}, Filtered Rayleigh Scattering (FRS)~\cite{doi:10.1366/000370208784909526,gustavsson2005filtered,Doll:17,Hoffman:96,forkey1993time,KRISHNA2021329}, Coherent Rayleigh-Brillouin Scattering (CRBS)~\cite{Pan2,Gerakis2, doi:10.1063/1.4959778}, Spontaneous Raman Scattering (SRS)~\cite{Karpetis:96,doi:10.1002/jrs.1965}, Coherent Anti-Stokes Raman Scattering (CARS)~\cite{Tolles:s,zheltikov,Lempert_2014,doi:10.1002/jrs.1489,Lempert:11,doi:10.1002/jrs.5339,doi:10.1063/1.3483871}, Laser Induced Fluorescence (LIF) and Planar LIF (PLIF)~\cite{hanson1990planar,Lee:93,crimaldi2008planar,jiang2012no,doi:10.2514/2.1939}.
In addition to diagnostics, advanced neutral and plasma manipulation concepts such as laser acceleration of neutral species~\cite{maher2012laser}, optical deceleration of atomic~\cite{https://doi.org/10.48550/arxiv.2206.07034} or molecular species~\cite{PhysRevLett.93.243004,Yang_2018} to even antihydrogen transport~\cite{Madsen_2021}, as well as dual laser pulse ignition schemes~\cite{doi:10.1063/1.5043295,f96eb5b913034b98b912310b23394b33}, can also benefit from the use of advanced laser sources.

Although the techniques mentioned above have already been demonstrated using mostly standard, commercial laser sources, their efficiency has often been subject to limitations posed by the performance specifications and capabilities of the laser sources used in the respective studies. Specifications such as pulse duration, energy per pulse, temporal shape of the laser pulses and, in some cases, achievable chirp of the pulses are mainly fixed parameters, essentially limiting the exploration and determination of the optimum operating conditions per use case.

To enable better efficiency and enhancement of the operational parameter space in each of these areas, a novel laser system has been developed, applicable (non exclusively) to all of the aforementioned regimes and techniques, and is presented here. This system is capable of achieving (i) variable pulse durations ranging from $\sim1$ ns to $\sim1$ $\upmu$s with (ii) arbitrary temporal profiles at (iii) a dual wavelength configuration, while (iv) exhibiting faster chirping rates to what has been reported in the relevant literature, up to a full chirping range of $\sim10$ GHz. Additionally, the system presented here (v) can be easily scalable, able to deliver pulse energies ranging from a few nanoJoules per pulse up to a few Joules per pulse, with no physical limitation on the achievable upper limit. A laser system like the one presented here can bring substantial advances to the fields discussed previously, by not only performing more optimally versus current laser designs but can also, and more importantly, lead to novel modi operandi of the techniques where it is being used.

\section{A dual-color, frequency and pulse duration agile laser system}
The overall optical layout of the laser system which features dual color, frequency and pulse duration agile capabilities is presented in Fig.~\ref{fig:lasersys}. A short overall description of the laser's operation will be provided first, before breaking down into the respective laser sub-systems, as these are noted in Fig.~\ref{fig:lasersys}.

\begin{figure}[h!]
\centering
\includegraphics[width=\textwidth]{1_Laser_system.pdf}
\caption{\label{fig:lasersys} Conceptual schematic of the dual color, frequency, pulse duration and shape agile laser system presented in this manuscript. The four main sub-systems of the laser are noted in Arabic numbers while the thickness of the line indicating the laser beam path is proportional to the laser energy/power of the beam at that particular point.}
\end{figure}

The laser chain starts with a custom designed and constructed, continuous wave (CW) $1064$ nm beam emitted from an Nd:YVO$4$ microchip laser, whose frequency can be rapidly chirped by an intracavity LiTaO$_3$ electro-optic crystal. This laser, which is referred to as the \textit{master}, is injection-locked onto a diode laser, termed the \textit{slave}, which produces a stable intensity CW beam while retaining the frequency information from the master laser. After its power is amplified to $\sim1$ W by a commercial fiber amplifier, this CW beam is split in two beams. Each beam is independently fed into a respective intensity electro-optic pulse shaper where the beam can be chopped up into pulses ranging from a temporal duration of $\sim1$ ns to $\sim1$ $\upmu$s with any arbitrary temporal shape, driven by an arbitrary waveform generator (AWG). These output optical pulses are each amplified by separate diode-pumped Nd:YAG amplifier chains, eventually producing two main $1064$ nm pulsed beams, each of $\sim450$ mJ energy per pulse. Part of one of the main beams is picked up to generate the third beam, which can be frequency doubled to generate a $532$ nm beam for the dual color configuration, maintaining the exact temporal characteristics as the first harmonic.

The described overall functions of the system are realized mainly by four sub-systems, as noted on Fig.~\ref{fig:lasersys}: 1) Master-slave injection locking, 2) Pulse shaping, 3) Nd:YAG amplifier chains and 4) laser beam output. Each sub-system of the laser system performs specific roles in the laser parameter configuration. The details of the laser sub-systems are discussed in detail in the following sections.

\subsection{Master-slave injection locking}
\subsubsection{A custom-built Nd:YVO$_4$ microchip master laser}

A neodymium based dioded-pumped solid state (DPSS) laser serves as a master laser in the developed system (Section 1 in Fig.~\ref{fig:lasersys}). The main components of the master laser cavity are an Nd:YVO$_4$ crystal,
acting as the lasing medium\footnote{$2.95$x$3.05$x$0.5$ mm $3$\% Nd doped, AR coated for $808$ nm on the incident side and HR coated for $1064$ nm on the cavity side; CASIX} and a plano-concave partially reflecting (PR) output coupler\footnote{$95\%$ PR@$1064$ nm on the cavity side; LAYERTEC} (OC). The Nd:YVO$_4$ crystal is optically pumped by an $808$ nm pumping diode\footnote{L$808$P$1000$MM; Thorlabs}. Most importantly, an intracavity LiTaO$_3$ electro-optic modulator\footnote{$1$x$3$x$5$ mm, AR coated for $1064$ nm on both laser propagation sides, gold plated on the two orthogonal sides; CASIX} (EOM), comprises almost the whole length of the laser cavity: by applying voltage to the crystal, its refractive index is changed and thus the cavity length, allowing for a fast change of the output laser frequency.

\begin{figure}[h!]
\centering
\includegraphics[width=0.7\textwidth]{2_cavity.pdf}
\caption{\label{fig:cavity}Schematic of the microchip laser cavity.}
\end{figure}

Figure~\ref{fig:cavity} presents a schematic of the microchip laser cavity of the master laser. A shorter laser cavity for the Nd:YVO$_4$ microchip DPSS laser is chosen, since it allows for a larger free spectral range (FSR) of the laser, thus ensuring single mode operation of the laser over the chirped duration of the final laser pulse. Apart from it being readily commercially available, the Nd:YVO$_4$ crystal is preferred as the active lasing medium over an Nd:YAG crystal since it has a higher stimulated emission cross-section and experiences strong, broadband absorption at $808$ nm\cite{Koechner1988PropertiesMaterials}. Additionally, due to its large absorption bandwidth, it is more tolerant to temperature variations of the pumping diode, enabling overall a more stable laser operation. It is noteworthy that a plano-concave (PC) OC is chosen for the cavity design: compared to a more common plano-parallel configuration which fully relies on thermal lensing for lasing, the PC OC allows for the lasing condition to fall well within the stable region of the stability criterion for the laser resonator. The cavity stability criterion is given as $0 \le g_1 g_2 \le 1$, where $g = 1-L/R$ with $L$ representing the cavity length and $R$ the radii of curvature (ROC) of the cavity surfaces, the subscripts representing each cavity forming surface~\cite{Koechner1988PropertiesMaterials}. With the flat surface on the Nd:YVO$_4$ side ($R_1 = \infty$), the PC OC equaling 100 mm ($R_2$ = 100 mm) and the cavity length equaling the length of the intracavity crystal ($L\approx5$ mm): $g_1 = 1$ and $g_2 < 1$, thus satisfying the stability criterion, ensuring stable laser operation.

To optimize the performance of the DPSS laser, two OCs having different partial reflectivities (PR) of 95\% and 97.5\%, respectively, were tested. An experimental test like this is preferred over an accurate simulation method for such DPSS lasers as it is much faster and more economical, requiring a direct characterization of the laser performance with actual parts in place. Figure~\ref{fig:IP} shows characteristic I-P curves of the output optical power (P) of the master laser versus the $808$ nm pump diode current (I). Initially, without the EOM in place (i.e. without an extra loss mechanism), the I-P efficiencies for $95$\% and $97.5$\% PR OCs are obtained, leading to the selection of the 95\% PR OC. Additionally, the final master laser efficiency with the $95$\% PR OC and the EOM in place was also measured. All measurements were performed when the laser cavity was aligned to produce a TEM$_{00}$ beam.

\begin{figure}[h!]
\centering
\includegraphics[width=0.7\textwidth]{3_IP_Curves.pdf}
\caption{\label{fig:IP}Characteristic master laser output power versus pump diode current IP curves for two different output coupler partial reflections: $95$\% (black) and $97.5$\% (blue). EOM ON/OFF indicates whether the EOM is in place or not. The IP curve for when the EOM was placed inside the cavity is shown in red. All curves are obtained for the pump diode at 14 \textcelsius~ which is a nominal temperature.}
\end{figure}

One of the important roles of the master laser in this laser system is to enable frequency chirping, i.e. for relative frequency excursions between the two eventually produced laser beams of the system. The frequency chirping is achieved by the intracavity LiTaO$_3$ EOM crystal. An electrical voltage via a waveform generator\footnote{3314A; Agilent} is amplified by a voltage amplifier\footnote{A-304; A.A. Lab system}, and this amplified voltage is applied onto the gold plated top and bottom side of the crystal to induce a change in the refractive index $n_1$ of the crystal with respect to time, and hence in its optical pathlength $n_1l_1$, where $l_1$ is its physical length. The applied voltage waveform can be of any waveform shape achieved by the waveform generator; usually a linear ramp is preferred, since this is translated to a linear chirp across the final laser pulse duration. The applied voltage on the intracavity crystal effectively equates to a change in the cavity optical path length $l$ and thus a change to the output laser frequency $f$, eventually leading to the laser frequency chirping. The achieved frequency change $\delta f$ in the master laser output with respect to the applied voltage $\delta V$ is given by\cite{Repasky2002}:
\begin{equation}
    \frac{\delta f}{\delta V} = \frac{f n_1^2 r_{33}}{2 d} \frac{n_1 l_1}{l_\mathrm{cav}},
    \label{eqn:chirprate}
\end{equation}
where $f$ is the optical frequency in the cavity, $n_1$ and $l_1$ are the refractive index and length of the LiTaO$_3$ crystal, and $l_\mathrm{cav}$ is the total optical path length given as $l_\mathrm{cav} = n_1 l_1 + n_2 l_2 + n_3 l_3$ where the subscripts 2 and 3 are for the Nd:YVO$_4$ and the air gap, respectively. Thus, as the total optical pathlength of the intracavity LiTaO$_3$ EOM crystal within the cavity increases, the frequency chirping rate (measured in MHz/V) increases.

\color{black}
\begin{figure}[h!]
\centering\includegraphics[width=0.6\textwidth]{4_Chirp_Rate.pdf}
\caption{\label{fig:chirpmod} Frequency chirp rate (in MHz/V) as a function of the modulation frequency of the applied voltage. The measurement was performed by applying a continuous, sinusoidal voltage of varying voltage amplitude and frequency on the crystal, and measuring the laser output through a Fabry-Perot interferometer.

\end{figure}

The frequency chirp rate of the present system as a function of the modulation frequency of the electric voltage applied to the intracavity LiTaO$_3$ EOM crystal was measured and is shown in Fig.~\ref{fig:chirpmod}. It should be noted that the modulation frequency referred to here is distinct from the laser (optical) frequency $f$ that is used when estimating the chirp rate. The modulation frequency is the frequency of the sinusoidal voltage signal sent across the EOM crystal, or the inverse of the period for a single chirp. The chirp rate, as defined in Eq.~(\ref{eqn:chirprate}), is the change in the laser frequency over the applied peak-to-peak voltage on the EOM crystal. A high chirp rate is desirable in order to achieve the widest frequency difference between the two eventual laser beams, typically on the order of a few gigahertz; it is favorable if this is achieved using the lowest voltage possible. To quantify the achieved chirp rate of this laser system, a Fabry-Perot interferometer is employed\footnote{SA200-8B, $1.5$ GHz bandwidth; Thorlabs}. The $\delta f$ measured for each individual voltage is proportional to the temporal width of the Fabry-Perot signal over its scanning window, measured as the full width half maximum (FWHM). The FWHM was measured for various peak to peak applied voltages and a linear fit performed, resulting in the chirp rate per modulation frequency. This process was repeated over various voltage modulation frequencies to find the optimal operating modulation frequency.

At low modulation frequencies, the chirp rate was found to steadily increase with higher modulation frequency, , increasing the rate at which the applied voltage changes the laser frequency. At approximately $700$ kHz modulation frequency the chirp rate starts to reach a plateau, reducing the effect of the modulation frequency on the chirp rate. Higher modulation frequency values lead to larger uncertainties in the chirp rate.
At the optimal operating conditions, determined through the process described above, a chirp rate of approximately $27$~MHz/V is obtained. For a $\pm 50$~V voltage sweep, this corresponds to an achievable frequency scanning range of $\pm 1.3$~GHz. This is approximately double the previously reported frequency chirping rate of $\sim15$~MHz/V for a similar system\cite{Gerakis2013}. Given the fact that we can supply up to $\pm 400$~V voltage sweep within the same pulse duration, we estimate that a full chirp range exceeding $10$ GHz can be achieved with this system. For this to be verified, a Fabry-Perot interferometer with a larger free-spectral-range (e.g. at $10$ GHz) is needed

We note that such rapid frequency chirping in a Nd based solid-state laser can induce relaxation oscillations\cite{Siegman1986Lasers}, which can lead to undesirable, strong laser output intensity oscillations. The relaxation oscillations are less strong with an Nd:YVO$_4$ crystal, since it has a shorter upper-state lifetime versus e.g. an Nd:YAG crystal\cite{Koechner1988PropertiesMaterials}. Nevertheless, the relaxation oscillations would cause severe intensity fluctuations at the final output of the laser system. To resolve this undesirable effect, a master-slave injection locking scheme was implemented and is described in the next section

\subsubsection{Slave laser}\label{sec:slave}
As mentioned in the previous section, with a master-slave injection locking scheme, undesirable laser intensity fluctuations due to the relaxation oscillations are avoided. This is because only the frequency information of the master laser is transmitted to the slave laser, while the intensity fluctuations are disregarded if the master laser's input power exceeds the slave laser's saturation power. When the frequencies of both the master and slave lasers are close enough, the chirped master laser forces the slave laser to emit at the exact same frequency through a strong, saturated coupling while maintaining a stable output intensity.

A $1064$ nm diode laser\footnote{M9-A64-0200; Thorlabs} is used as a slave laser with benchtop controllers\footnote{LDC205C and TED200C; Thorlabs} controlling its operating current and temperature. Temperature tuning of the laser diode is a crucial step in the process, since it enables the laser diode to reach an emitting wavelength close to that of the master laser: this is a key to a successful injection locking scheme. To prevent any back-reflection of the slave laser output into itself, an optical isolator\footnote{110-21022-0001-C1; Electro-Optics Technology} is placed in front of the slave laser. To maximize the coupling between the master and slave lasers, the output of the master laser is fed into the slave laser through a rejecting port of the optical isolator (see Fig.~\ref{fig:lasersys}). This improves the total transmission of the master laser into the slave laser. The status of the injection locking is monitored by a Fabry-Perot interferometer, which monitors the characteristic shape of the power spectrum of the injection-locked laser\cite{Troger1999Frequency-sweepLocking}.

\subsection{Temporal Pulse Shape Manipulation}
The frequency chirped CW laser output from the master-slave section is sent to the pulse shaping section of the laser system. Prior to being shaped into laser pulses, the CW beam coming from the master-slave section is initially further amplified by a commercial, Yb-based fiber amplifier\footnote{YAR-10K-1064-LP-SF; IPG Photonics} that is capable of emitting up to $10$ W power at the same frequency. The amplified output is split in two through a $50/50$ beamsplitter cube, with each resulting beam being fed into a respective commercial Mach-Zender type intensity EOM\footnote{AM1064b; Jenoptik} (iEOM). Each iEOM is used to shape a pre-amplified seed pump pulse. It should be noted that the adaptation of this type of iEOM enables the pulse shaping with commonly available arbitrary waveform generators since the required voltage for shaping is relatively low (range of -$5$V to $5$V). This advantage allows for avoiding the use of high voltage switching equipment such as Pockel's cells, which typically operate in the kilovolt range and induce a strong electromagnetic interference noise.  Pulse shaping begins by shaping the first pulse (Beam 1) with an iEOM, where initially the unchirped, constant frequency part of the laser beam is sent through the iEOM. At a later time, an electrical voltage is applied to the intracavity EOM in the master laser when shaping the second pulse (Beam 2) with the second iEOM: this enables inducing a relative frequency difference between the two pulses. This relative frequency difference is solely determined by the voltage waveform applied to the intracavity EOM. Additionally, since both pulses originate from the same master laser, as long as the master-slave laser system is stable within the hundreds of ns of the pulse shaping operation, the overall operation of the laser system will also be stable. Finally, to temporally align the two produced pulses, the first pulse goes  through an $85$ m optical fiber delay line.

An arbitrary waveform generator\footnote{AFG3102C; Tektronix} is used to drive the iEOMs to shape and manipulate the pulses from the CW beam. This provides the capability for additional delay control. In addition, and importantly, the arbitrary temporal pulse shaping allows for the generation of a flat top output pulse by accommodating the time dependent gain in the Nd:YAG amplification chains (discussed in Section~\ref{sec:amplifiers}).

\begin{figure}[h!]
\centering\includegraphics[width=\textwidth]{5_Pulse_shaping.pdf}
\caption{\label{fig:pul_mani} Temporal property manipulation of the two pulsed beams: (a) Manipulation of relative time delay between pulses and (b) manipulation of temporal pulse duration.}
\end{figure}

Temporal pulse delay manipulation and pulse duration manipulation examples are shown in Fig.~\ref{fig:pul_mani}. In Fig.~\ref{fig:pul_mani}~(a), the temporal position of Beam 2 is fixed in time, while the relative timing of Beam 1 is changed. In Fig.~\ref{fig:pul_mani}(b), the temporal duration of both pulses is manipulated by controlling the period of the voltage signal applied to the iEOMs. It is noteworthy that each pulse can be individually controlled, both in temporal shape/duration as well as in their relative delay. The obtained curves are taken by a fast photodiode located before the Nd:YAG amplifier chains (heterodyne detection $1$ in  Fig.~\ref{fig:lasersys}). Examples of arbitrary temporal pulse shape manipulation are shown in Section~\ref{sec:laser_output}, where the temporal shapes before and after the second harmonic generation (SHG) are compared (Fig.~\ref{fig:shg_shape}).

The freedom in the pulse delay control from the iEOMs enables the accurate control of the produced beams simultaneously in the test section. This is especially beneficial for applications involving multiple-wave mixing techniques (such as e.g. single shot CRBS) where all beams are required to propagate the scattering volume at the same time. Alternatively, two separate pulsed beams can also serve as individual laser sources for applications such as dual-pulse laser ignition~\cite{Mahamud2018Dual-pulseModel}. The controllable pulse duration is also highly advantageous in applications such as the ones involving strong background optical radiation. In these cases, the short probing pulse can be used to minimize the background collection for better signal-to-noise ratios. On the other hand, a longer duration pulse may be desirable for applications where the high peak laser intensity can be problematic. For instance, a longer duration pulse for the same laser pulse energy can avoid undesirable laser perturbation such as laser heating or photo-ionization, which could be problematic in e.g. Thomson scattering measurements on low density plasmas\cite{Carbone2015ThomsonChallenges}.

\subsection{Nd:YAG amplifier chains}
\label{sec:amplifiers}
Each of the two generated pulsed beams undergoes pulsed amplification by propagating through two arms, each comprising of a total of four diode pumped Nd:YAG amplifier stages in tandem. In each of these arms, the laser pulse at an initial energy per pulse of $\sim$100~nJ is amplified to a final energy per pulse of $\sim$450 mJ after all amplification stages. We note that the diode pumped amplifiers used here\footnote{In tandem in each stage: one RBAT20, two RBAT34, and one REA7006 from Northrop Grumman Cutting Edge Optronics} have several advantages compared to flash lamp pumped amplifiers since they are more compact in size, silent, and have good electrical-to-optical efficiency. Each amplifier stage in the chain is set in a double-pass single-stage amplifier configuration, with each stage consisting of a Faraday optical isolator, a half-wave plate, a thin-film polarizer, a diode-pumped Nd:YAG rod, a quarter-wave plate, and a high reflective mirror for $0\degree$ angle of incidence (refer to section 3 of Fig.~\ref{fig:lasersys}). The Nd:YAG amplifiers in each stage increase the energy from $\sim100$~nJ to $\sim12$~$\upmu$J in the first stage, to $\sim800$~$\upmu$J in the second, $\sim20$~mJ in the third, and finally to $\sim450$~mJ in the fourth, final stage. Additionally, we note that along the amplifier chain, the laser beam diameter is properly adjusted with the use of beam expanders to a beam diameter corresponding to $\sim80$\% of the Nd:YAG rod diameter as recommended by the manufacturer. By selectively operating the respective amplification stages, the final output laser pulse energy is easily controllable while, with the addition of more amplifier stages, the energy per pulse achievable with this system is highly scalable.

\subsection{Laser beam output section}
\label{sec:laser_output}
\subsubsection{First harmonic generation}
Throughout the previously described sections of the laser system, the initial first harmonic CW beam eventually provides two fundamental first harmonic $1064$ nm amplified beams (Beam 1 and Beam 2 in Fig.~\ref{fig:lasersys}) having a Gaussian spatial beam shape with a beam diameter of $\sim5.6$ mm with variable pulse energy and duration, which can be utilized for various applications. For instance, the typical characteristics of the beams used for CRBS application (discussed in Sec.~\ref{sec:apps}) are as follows: a pulse energy of $\sim450$ mJ, a pulse duration of $\sim 150$ ns with a flat-top temporal profile.

\begin{figure}[h!]
\centering\includegraphics[width=1\textwidth]{6_Hetero.pdf}
\caption{\label{fig:hetero} Heterodyne signal (solid black line) and analyzed chirped frequency (red dashed line) are shown for (a) before and (b) after amplification. This result demonstrates that the pulse amplification is not altering the chirp as set by the master laser.}
\end{figure}

Each beam is sampled and heterodyned on a fast InGaAs photodiode\footnote{G6854-01; Hamamatsu Photonics} at the heterodyne detection location (see Section 4 of Fig.~\ref{fig:lasersys}). The heterodyned signal of the two beams provides information on the final, relative chirped frequency between the two beams. Sample heterodyne signals before and after the amplification stages are shown in Fig.~\ref{fig:hetero}. The obtained beat patterns contain the chirped frequency information and are analyzed as described in Ref.~\cite{Fee1992OpticalSpectroscopy}. In Fig.~\ref{fig:hetero}, the extracted chirped frequency based on this method, for the heterodyned signal before and after amplification, is shown.
In the sample heterodyne signals shown in Fig.~\ref{fig:hetero}, the relative frequency difference between the two beams is initially large, then becomes zero at $t=0$, and then increases again: this practically translates to an optical lattice starting from a high velocity, decelerating to zero relative frequency difference (velocity of $0$ ms$^{-1}$) and accelerating again to a higher velocity

The repeatability and stability of the chirp rate before and after amplification is checked and shown in Fig.~\ref{fig:chirp}. The relative standard deviation over the total 19 shots was found to be $<4\%$ before the amplification and $<0.1\%$ after the amplification. The smaller deviation after the amplification is attributed to a better signal-to-noise ratio. The good agreement in both the heterodyne signals and the chirp rates before and after the amplification indicates that the amplification chains do not introduce any disturbance in the achievable chirp rate of the system.
\begin{figure}[h!]
\centering\includegraphics[width=0.6\textwidth]{7_Chirp_vs_shots.pdf}
\caption{\label{fig:chirp} Temporal chirp rate repeatability for 19 total sample shots are shown. Circles are from the pulses before amplification, and triangles are from those after amplification. Temporal chirp rates during both the frequency deceleration and acceleration phases are shown.}
\end{figure}

\subsubsection{Second harmonic generation}
For a dual wavelength configuration, part of Beam $2$ is sampled to produce the second harmonic beam (Beam 3 in Fig.~\ref{fig:lasersys}). The second harmonic beam is generated by a LBO crystal\footnote{$6$x$6$x$50$ mm; AR@$1064$+$532$ nm; Eksma Optics} operating with a commercial crystal oven\footnote{CO1-50-6/6; Eksma Optics} and a temperature controller\footnote{TC2; Eksma Optics}. The second harmonic $532$ nm probe energy per pulse as a function of the $1064$ nm beam input energy per pulse on the crystal is shown in Fig.~\ref{fig:shg}: approximately $27$\% conversion efficiency is obtained.

\begin{figure}[h!]
\vspace{-5mm}
\centering\includegraphics[width=0.7\textwidth]{8_Crystal_efficiency.pdf}
\vspace{-2mm}
\caption{\label{fig:shg}Second harmonic probe energy as a function of $1064$ nm probe beam. The data points are represented as square points along with the fitted curve for a visual guide.}
\end{figure}

It should be noted that, even though there is a natural loss in the achievable energy per pulse, a dual color option can effectively provide a stronger final signal for several applications such as linear laser light scattering measurements. For example, Rayleigh scattered signals scale as $\lambda^{-4}$, where $\lambda$ is the wavelength of the laser used. As a result, the signal for the $532$ nm probe beam will be sixteen times stronger than the $1064$ nm signal at the same energy. Additionally, detectors with a much higher quantum efficiency are available for $532$ nm than those for $1064$ nm. Another practical advantage is that, by adopting a dual color configuration, i.e. $1064$ nm pump beams and a $532$ nm probe beam in any phase matching required applications, the signal beam no longer shares the same optical path with any of the pump beams, which is often the case in a single color co-planar phase matching configuration; an example is single color coherent Rayleigh-Brillouin scattering (CRBS) versus dual color CRBS. All these advantages lead to a significant improvement of the achievable  signal-to-noise ratio.

\begin{figure}[h!]
\centering\includegraphics[width=1\textwidth]{9_Harmonic_four_mix.pdf}
\caption{\label{fig:shg_shape}Arbitrary temporal shapes of the third beam generated by the laser system presented here. Both a $1064$ nm beam, before the second harmonic generation stage (red line), and a $532$ nm beam, after the second harmonic generation stage (green line), are shown. Note that the slight differences on the temporal shapes is attributed to the different sensitivities of the photodiodes used for the respective wavelengths.}
\end{figure}

Importantly, we demonstrate here that the generated Beam $3$ at $532$ nm can maintain the temporal characteristics of the first harmonic. Figure~\ref{fig:shg_shape} shows several sample arbitrary temporal shapes for the second harmonic probe beam compared with the $1064$ nm first harmonic, normalized by the peak intensity. It is notable that the temporal shape is well preserved, thus enabling great flexibility on pulse temporal manipulation for both harmonics. The slight difference on the observed temporal shapes is attributed to the different sensitivity of the photodiodes used for the respective wavelengths. Finally, we note that both spatial beam profiles have good Gaussian shapes as shown in Fig.~\ref{fig:shg_profile}.

\begin{figure}[h!]
\centering\includegraphics[width=\textwidth]{10_before_after_shg.pdf}
\caption{\label{fig:shg_profile}Spatial beam profiles of (a) the 1064 nm probe beam and (b) the second harmonic 532 nm. Horizontal and vertical distributions in mm (along the white dashed lines) are given on the bottom and right sides.}
\end{figure}

\section{Applications of the laser system}\label{sec:apps}
The flexibility and capabilities of the laser system presented here are expected to be useful and potentially enabling for a variety of laser-aided diagnostics and optical manipulation experiments. Particularly, applications with optical lattices~\cite{Greiner2008OpticalLattices, maher2012laser, PhysRevLett.93.243004, Madsen_2021, PhysRevA.100.023401} that are created by the interference of two laser beams can benefit from the rapid frequency chirping capability provided with this all-optical setup, in comparison to the mechanical setups that are usually the current state-of-the-art. Additional applications can also include gas diagnostics such as single shot CRBS\cite{Gerakis2013,Gerakis2018,Gerakis2021}, micropropulsion by molecular manipulation\cite{Shneider2006}, optical microlinear deceleration and acceleration\cite{Barker2001,Barker2002SlowingDeceleration}, and particle spectroscopy\cite{Gerakis2022RamanLattices}. With a controllable pulse duration, the present system coupled with gated detection may benefit applications with accompanying strong background radiation such as spark or arc discharges\cite{Yatom2017DetectionIncandescence,Gerakis2018}. In addition, the system could be of benefit to applications where a longer pulse is desirable but the peak laser intensity is problematic, such as laser heating\cite{Carbone2012RevisionLAPD} and photo-ionization\cite{Carbone2015ThomsonChallenges}

With regard to CRBS in particular,, the flexibility of the pulse duration and chirping rate can accommodate a wider range of spectral lineshape changes resulting from temperature and pressure variations~\cite{Liang2017} or from the existence of a bulk flow velocity~\cite{Gerakis2021,Bak2021AchievingScattering}. Additionally, the dual color single shot CRBS configuration using a $532$ nm probe beam will enable access to gas information at much lower pressure applications than ever reported in the past~\cite{Graul2014,Pan2002b,Vieitez2010CoherentMixtures,Wang2020}. Furthermore, each beam produced by the system presented here could be utilized for a plethora of different applications such as laser induced plasma formation~\cite{LeDrogoff2004InfluenceProperties}, surface ablation~\cite{Amoruso1997AbsorptionPlasmas}, and dual-pulse laser ignition applications\cite{Shneider2011TailoringPulse,Mahamud2018Dual-pulseModel}, among others.

\begin{figure}[h!]
	\centering\includegraphics[width=\textwidth]{11_CRBS_sample.pdf}
	\caption{\label{fig:CRBScompa}Sample experimental CRBS lineshapes of CO$_2$ gas at room temperature and atmospheric pressure; (a) with a $1064$ nm probe beam and (b) with a $532$ nm probe beam. The dashed-dot line indicates the Rayleigh peak centered at $0$~m/s, and the black dashed line indicates the speed of sound for CO$_2$ ($267$~m/s at $293$~K and pressure of $1$~atm, estimated by the sound speed equation $v_\mathrm{sound}=\sqrt{\gamma R T}$, where $\gamma$ is the adiabatic constant, $R$ is the specific gas constant, and $T$ is the gas temperature).}
\end{figure}

As a demonstration of the capabilities of present laser system, we utilize all beams produced by the system to obtain single shot CRBS lineshapes, in single and dual color CRBS geometries. In this application, the two frequency chirped pump beams at $1064$ nm are crossed in air at an angle of $\sim178\degree$ to generate temporally decelerating and accelerating optical lattices within the $\sim150$ ns duration of the pulse: this enables the interaction of the optical lattice with the populations in all velocity domains of the velocity distribution function, within a single laser shot.
The third probe beam, either at $1064$ nm or at $532$ nm, is incident on the lattice at the Bragg angle ($89\degree$ and $\sim30\degree$ respectively for $1064$ nm and $532$ nm), generating a fourth signal beam containing the single shot CRBS lineshape that provides thermodynamic information of the gas properties, such as temperature, density and flow velocity. Figure~\ref{fig:CRBScompa} shows sample CRBS lineshapes of CO$_2$ gas at room temperature and atmospheric pressure, where (a) is obtained with a $1064$ nm probe beam, and (b) is obtained with a $532$ nm probe beam. Note that the lineshapes are averaged over five shots. The dashed-dot line and the black dashed line indicate the Rayleigh peak centered at $0$~m/s and the speed of sound in CO$_2$ ($267$~m/s at $293$~K and pressure of $1$~atm), respectively. Successful CRBS lineshape acquisition with the $532$ nm probe suggests the accessibility toward much lower pressure (gas density) regimes for gas property diagnostics, using the dual-color single shot CRBS configuration. This can in turn enable advances e.g. in low pressure plasma applications for determining departures from local thermodynamic equilibrium condition in a plasma, or for flow velocimetry in a hypersonic wind tunnel.

\section{Conclusion and future directions}
In this manuscript, a dual color, frequency and pulse duration agile laser system is presented. The system consists of four subsections where specific laser parameters are configured: master-slave injection locking (frequency manipulation), pulse shaping (pulse temporal duration and shape manipulation), Nd:YAG amplifier chains (pulse energy manipulation), and laser beam output (second harmonic generation). In summary, the system is capable of (i) achieving variable pulse durations ranging from $\sim1$ ns to $\sim1 \upmu$s with (ii) arbitrary temporal profiles at (iii) a dual wavelength output configuration, (iv) exhibiting faster chirping rates to what has been reported with chirp range up to $10$ GHz, and (v) with scalability of pulse energies from some nJ per pulse up to Joules per pulse. These capabilities of the presented laser system are demonstrated by obtaining single shot CRBS spectra in both single color and dual color configurations.
It is envisioned that the laser system could be of benefit to various applications for spectroscopy, gas and plasma diagnostics as well as particle manipulation, as one can tailor the specific laser parameters listed above with ease, as would be required by each application.

\section{Acknowledgements}
This material is based upon work supported by the National Science Foundation under Grant No. $1903481$.

\newpage


\section{Introduction}
The thermodynamical features of horizons have been extensively studied ever since Hawking predicted that black holes create particles in their vicinity \cite{haw75}. His seminal work inspired a wealth of research aiming at understanding the occurrence of quantum processes at causal boundaries. Although horizons act classically as semipermeable membranes, allowing only a one-way crossing, quantum effects might follow classically forbidden processes and exit the region trapped by the horizon. Hawking tentatively described particle production as a consequence of this through quantum tunneling, but it took twenty years until such description was made manifest \cite{par00,pad02}. In this picture, particles escape the causal enclosure of the horizon through gravitational tunneling on complex paths. The formalism proposed by \cite{par00,pad02} was substantiated by \cite{mas00}, and finally connected to the Hamilton--Jacobi formalism by \cite{di07,van11}, where thermodynamical properties were shown to be in agreement with previous results that use Bogolubov coefficients. In contrast to the latter, the Hamilton--Jacobi method links the pole structure of the semi-classical amplitude at the horizon to the imaginary part of the classical action, such that the tunneling rate can be compared with the Boltzmann distribution, from which thermal properties follow immediately.

The prevalence of horizon thermodynamics underlines how deeply this effect is enmeshed in the foundation of quantum field theory. Although seemingly simple, the tunneling picture offers various insights into the local quantum processes that govern the vicinity of the horizon. Especially after its application to dynamical horizons \cite{hay94,ash99, ash00,ash03}, this method offers a huge advantage over other approaches, since these concepts of horizons are quasi-local and independent from global properties such as asymptotic flatness. Due to its very explicit usage of the paths across the horizon, the tunneling picture allows to find generalizations in a very systematic way \cite{sen15,Kant:2009pm}. It also suggests a generalized notion of the Hawking effect \cite{gia20} that applies to all types of horizons, dynamical and static, and draws a connection to the consistency of quantum field theory.

In this article we extend the tunneling method to the case of universal horizons (UH) in Lorentz violating theories of Gravity. In particular, we study Einstein--Aether (EA) gravity \cite{Jacobson:2000xp}, where boost invariance is explicitly broken by the presence of a time-like unit four-vector $U^\mu$, the aether, which propagates a vector and a scalar degree of freedom on top of the usual transverse traceless graviton, when coupled to the Einstein--Hilbert action \cite{Jacobson:2004ts}. As a consequence of this, the symmetry of the theory is reduced -- once the time direction is identified with the integral lines of $U^\mu$ -- down to foliation preserving diffeomorphisms (FDiff), consisting of the direct product of time reparametrizations and time-dependent spatial diffeomorphisms
\begin{align}\label{eq:FDiff}
    \tau\rightarrow \tau'(\tau),\quad x^i\rightarrow x'^i(t,\tau)
\end{align}
where $\tau'(\tau)$ must be a monotonous function\footnote{From now on we use greek indices to denote space-time coordinates -- including time -- and latin indices to denote spatial directions only. We also use a mostly plus convention for the metric signature.}.

This allows for matter actions coupled to gravity to include higher derivative operators along \emph{spatial directions}, while keeping only two time derivatives, thus avoiding the presence of Ostrogradsky ghost instabilities. As a consequence, matter fields coupled to EA gravity propagate with modified dispersion relations of the generic form
\begin{align}
    \omega^2=k^2+\alpha_4\frac{k^4}{\Lambda^2}
    +\dots +\alpha_{2Z}\frac{k^{2Z}}{\Lambda^{2Z-2}},
\end{align}
which allow for superluminal propagation at high momentum $k> \Lambda$.

The former statement seems to imply that the usual notion of horizon looses its meaning, since rays of these fields can freely enter and exit the region enclosed by the event horizon of a black hole. However, this naive intuition is broken in the case of known stationary and spherically symmetric black hole solutions in EA Gravity \cite{berg12}. In these, the aether is hypersurface orthogonal, and thus it defines a preferred foliation in co-dimension one hypersurfaces, described by a scalar field $T$
\begin{align}\label{eq:aether_ortho}
    U_\mu = \frac{\partial_\mu T}{\sqrt{|\partial_\alpha T \partial^a T|}},
\end{align}
which has the interpretation of a preferred time direction that every motion has to follow. Incidentally, these are also solutions to the low energy action of Ho\v rava--Lifshitz Gravity \cite{Horava:2009uw,Blas:2009qj}, known as Khronometric Gravity \cite{Blas:2010hb}. The latter propagates only one extra scalar field, in contrast to the vector and scalar propagated by EA gravity. However, in this work we will not discuss the dynamics of gravitational perturbations, and thus all our results can be applied indistinctively to either EA or Ho\v rava--Lifshitz Gravity.

In certain space-times -- which include the black hole solutions of interest --, it might happen that a foliation leaf becomes a constant radius hypersurface. Since the foliation defines a preferred time direction that all world-lines must follow, this implies that nothing can escape the region enclosed by the leaf. Such a hypersurface is thus named \emph{universal horizon} and it is locally characterized by the following two properties \cite{Bhattacharyya:2015gwa}
\begin{align}
    (U\cdot\chi)=0,\quad (U\cdot a)\neq 0,
\end{align}
where $\chi^\mu$ is the time-like Killing vector of the metric and $a_\mu=U^\alpha\nabla_\alpha U_\mu$ is the acceleration of the foliation\footnote{We have defined the dot-product between two vectors $X^\mu$ and $Y^\mu$ with respect to the space-time metric to be $(X\cdot Y)=g_{\mu\nu}X^\mu Y^\mu$.}. The usual notion of a horizon as a casual boundary thus gets resurrected in these theories by the existence of the UH. Even those motions travelling at infinite speed will be forever future-trapped in the inner region once they cross the UH. Due to this, one naively expects that UHs must radiate, have an entropy and, in general, replicate all the features that are usually associated to horizons in Lorentz invariant theories.

Indeed, the thermodynamical properties of UHs have been studied in the recent past \cite{he21,mi15}. In \cite{ber13,Ding:2015fyx} the tunneling picture was tentatively applied to UHs. However, despite showing a proof of concept and anticipating the same result we shall find for the Hawking temperature, we feel that in these early works several technical and physical subtleties -- arising from the recently improved understanding on the shape and dynamical role of the preferred foliation -- were not fully spelled out and require a more in detailed discussion.

In particular, it has been recently shown in \cite{del22} that the two regions bordering the universal horizon must feature a smooth lapse, since otherwise the theory looses its predictability on a fundamental level. This however contradicts the naive expectation that the whole space-time -- exterior and interior -- should be described by the same global foliation or, at the very least, by a discrete set of discontinuous foliations that are, however, globally oriented in the same way. This property manifests itself in the pole structure at the horizon, so that the construction of a well-defined tunneling path is compromised and needs to be revisited. As we shall see in what follows, our renewed analysis of the gravitational tunneling in this framework exposes quite relevant physical lessons concerning the origin and nature of the Hawking radiation in these settings.

In this paper we aim to solve these issues by investigating and formalizing the tunneling picture for UHs in Lorentz-violating theories in a solid way. Our work is organized as follows: First, in Section \ref{sec:relat_tunneling} we summarize and review the standard derivation of Hawking radiation in general relativistic space-times through the tunneling method. Later, in Section \ref{sec:non_re_tunneling} we show how this can be applied to the case of UHs in a straightforward way, by studying arbitrary spherically symmetric solutions and the trajectories of gravitating fields therein, by means of the Hamilton-Jacobi method. Due care will be given to the choice of the privileged time direction in Section \ref{sec:sync_factor}, which influences the temperature of the emitted radiation through the UV behavior of the dispersion relation of the field. Finally, we draw our conclusions in Section \ref{sec:conclusions}.

\section{Relativistic Tunneling}\label{sec:relat_tunneling}
Before we explore the trenches of Lorentz-violating theories, we briefly review the tunneling picture in theories with full invariance under diffeomorphisms. This allows us to introduce the basic concepts that we will later extend to Lorentz-breaking theories. Most importantly, it will helps us clarifying how the tunneling paths are identified, the Hamilton-Jacobi method to describe the tunneling rate, and the consistency criterion for a well-defined Hawking effect. Our review is based on \cite{gia20}, which discusses the relativistic framework generically also for dynamical cases. In the following we restrict ourselves to spherically symmetric and static space-times. However, generalizations are immediate.

\subsection{Gravitational tunneling}

Let us consider a static and spherically symmetric space-time. Its metric can always be written in terms of outgoing (ingoing) rays of light by using Eddington--Finkelstein--Bardeen (EFB) coordinates
\begin{align}\label{efb}
\mbox{d}s^2=-F(r)\mbox{d}\xi^2\pm2 B(r)\mbox{d}\xi\mbox{d}r+r^2\mbox{d}^2S_2,
\end{align}
where d${}^2S_2$ is the line element of the two-dimensional sphere of unit radius, and $\xi=\{u,v\}$ is a placeholder for the light-cone coordinates, $u=t-r^*$ and $v=t+r^*$, with $r^*=\int_r\tfrac{{\rm d}r}{F(r)}$ the tortoise coordinate. Here the functions $F(r)$ and $B(r)$ determine completely the metric. The existence of a horizon is determined by the condition $F(r)=0$, with $F(r)\in\mathcal{C}^2(\mathbb{R})$ at the least. The sign in front of the non-diagonal term is negative (positive) for $\xi\equiv u$ ($\xi\equiv v$).

The horizon, separating trapped and normal regions, acts as a semipermeable surface allowing only penetration from one side. Which paths are classically allowed depend on the nature of the trapped surface. For a black hole, the interior becomes future trapped and only ingoing trajectories can cross the horizon. This can be seen from the null-congruences in Schwarzschild space-time, that generate the lightcone. While the ingoing congruence $l^-=-2\partial_r$ is inwards pointing and well defined across the horizon, the outgoing one $l^+=\partial_v+F(r)\partial_r$ changes its direction from outgoing to effectively ingoing, when travelling from the (exterior) normal region into the (interior) trapped region. This is due to the fact that $F(r)$ changes its sign across the horizon. All curves in the interior are hence classically destined to stay inside.

Similarly, curves outside the horizon cannot classically have departed from within it. Indeed, they show a singular behavior in their momentum and classical action at the horizon. For example, for the radial momentum $k_r$, one finds through the Hamilton--Jacobi method\footnote{Alternatively, one can use the expansion parameters $\theta^V=h^{ab}\mathcal{L}_Vh_{ab}$ of a hypersurface with metric $h_{ab}$ along a vector field $V^\mu$, that are negative along ingoing and positive along outgoing directions. At the horizon, where one direction gets marginally trapped, we get $\theta=0$. For a Schwarzschild black hole, the outgoing direction along $l^+$ becomes trapped, hence, $\theta^+$ changes sign and vanishes when $F(r)=0$. The classical action $S_0$ can be rewritten using $\theta^+$ such that the spatial momentum $k_r\propto1/\theta^+$, thus inducing a pole in the integral \cite{sen15}.} that $k_r\sim F^{-1}(r)$. A small complexification around the horizon $F(r)\to F(r\pm i\varepsilon)$ can be introduced so to resolve this problem, but this also  leads to  a complex action, given that the imaginary part of the classical action
\begin{equation} \label{eq:clas_action_C}
\mathcal{S}_0=-\Omega v+\int \frac{2\Omega}{F(r\pm i\varepsilon)}\mbox{d}r,
\end{equation}
is then determined by the Sokhotski--Plemelj theorem\footnote{\label{ft:footnote3} The Sokhotski--Plemelj theorem is a relation for integrals with a small complexification around a pole
\begin{equation}
\lim_{\varepsilon\to0}\int_a^bdx\ \frac{f(x)}{x\pm i\varepsilon}=\mp i\pi f(0)+\mathcal{P}\left(\int_a^b\frac{f(x)}{x}\right),
\end{equation}
where $\mathcal{P}$ denotes the Cauchy principal value. An alternative path can be taken by complexifying the Hamilton--Jacobi equation by introducing the Feynman $\pm i\varepsilon$-prescription \cite{gia20}.}. Thus, an outgoing trajectory cannot start in the interior of the black hole. The classical one-particle action would become complex along this path, making it classically forbidden \cite{van08}. However, things dramatically change once quantum effects are considered, as  prohibited processes can become quantum-mechanically allowed. In particular, horizon crossing due to gravitational tunneling becomes possible.

Although conceptually different, it is instructive to look first at quantum mechanical tunneling through a finite size barrier. Suppose that we have a bounded one-dimensional $L^1$-integrable classical potential, such that $\lim_{x\rightarrow \pm\infty}V(x)=0$. The probability for a quantum mechanical wave function $\Psi(x)$ with energy $E<{\rm max}\left(V(x)\right)$ to tunnel through it is given by the ratio of the norm of the transmitted wave $\Psi_>$ divided by the norm of the incident wave $\Psi_<$ before hitting the wall
\begin{align}
    |T|^2=\frac{\|\Psi_{>}\|^2}{\|\Psi_{<}\|^2}.
\end{align}

For any shape of the potential, as long as the potential itself is only mildly time-dependent, this can be described by using a WKB-approximation for the wave function
\begin{align}
    \Psi(x)=\Psi_0\exp\left(\tfrac i\hbar\sum_i \mathcal{S}_i(x)\hbar^i\right).
\end{align}

To leading order, and considering the Schr\" odinger equation, the solution for ingoing and outgoing waves is determined by the classical action $\mathcal{S}_0=\pm\int$d$x\sqrt{2m(E-V(x))}$ with $E$ being the energy of the incident wave and $m$ its mass. Thus, one can understand the tunneling process as occurring through a classically forbidden path where the exponent becomes complex -- since inside the potential $E<V(x)$. This has some similarities with the case of gravitational paths, as we see below.

Let us start by noting that gravitational tunneling leads to particle creation by horizons. This has been formulated in \cite{par00,pad02} for the static black hole, and then further refined to a vast variety of examples \cite{van11}. A cartoon picture of the process was already suggested by Hawking \cite{haw75}. Close to but still outside the horizon, a Hawking pair consisting of a positive and a negative energy particle, can be created. The positive energy particle escapes the gravitational well and it is measured by an asymptotic observer, while the negative energy one falls into the black hole (where its existence on-shell is allowed by the spacelike nature of the Killing vector associated to stationarity in time).

Courtesy of the energy budget provided by the black hole, we can describe this from a different perspective. Instead of a pair creation, we interpret the same process as a single particle coming from the interior and tunneling outwards in a quantum mechanical way. This is possible because a) we can always trade an inward-pointing negative Killing energy particle with a positive energy, outward-pointing one; and b) he notion of energy positivity is linked to the Killing vector, which changes its character at the horizon.

This idea can be made explicit by considering a massive scalar field $\phi$ in the space-time described by \eqref{efb} and obeying the Klein-Gordon equation
\begin{equation}\label{kgrel}
\left(\Box-\frac{m^2}{\hbar^2}\right)\phi(x)=0.
\end{equation}
Assuming that the space-time does not change too rapidly, we can use a WKB ansatz for the field
\begin{equation}\label{wkb}
\phi=\phi_0e^{\frac i\hbar\sum_n\hbar^n\mathcal{S}_n(x)}=\phi_0 e^{\frac i\hbar \mathcal{S}_0+\mathcal{O}(\hbar^0)},
\end{equation}
with $\mathcal{S}_0(x)$ the classical action. The constant $\phi_0$ is allowed to have a mild coordinate dependence, but it is usually treated as effectively constant. We can then  plug \eqref{wkb} into \eqref{kgrel} to get the Hamilton--Jacobi equation to lowest order\footnote{In principle we could also consider the sub-leading term $\mathcal S_1$. However, this does not contribute to the tunneling probability.} in $\hbar$
\begin{equation}\label{hje}
\partial_\mu\mathcal{S}_0\partial^\mu\mathcal{S}_0+m^2=0.
\end{equation}

One can formulate a general ansatz for the classical action as $\mathcal{S}_0=\int $d$x^\mu\partial_\mu \mathcal{S}_0$ with d$x^\mu=\sum_K K^\mu (K_\nu {\rm d}x^\nu)/|K|$ where the vectors $K^\mu$ determine the observer, or in other words, the vacuum state. They must be such that they span the full chart of space-time coordinates. Generally, some of the $K^\mu$ are chosen to reflect the symmetries of the setup, i.e. Killing (static space-times), Kodama (spherically symmetric but dynamical), or dual-null vector (general space-times). In the case of a time-like Killing vector $K^\mu\equiv \chi^\mu$, this introduces a covariantly conserved energy, or frequency
\begin{equation}\label{frq}
\Omega:=-\chi^\mu\partial_\mu \mathcal{S}_0,
\end{equation}
that defines a preferred notion of observers. The ansatz of $\mathcal{S}_0$, thus, drastically simplifies to
\begin{align}\label{eq:S0_relativ}
    \mathcal{S}_0=-\Omega v + \int k_r \mbox{d}r,
\end{align}
where we have assumed spherical symmetry. Hence, the trajectory is completely determined by the Killing energy and the function $k_r(r)$, denoting the spatial momentum.

Coming back to our discussion on the quantum-mechanical tunneling, in a system with gravitational tunneling the probability to reach a classically inaccessible region is given through complex paths, which is reminiscent to the case of the one-dimensional potential barrier. The positive energy particle inside a black hole is now interpreted to take a generically complex path across the horizon determined by $\mathcal{S}_0$, as shown in figure \ref{eins}.

The particle first travels along a past-directed null-curve to the horizon which becomes an outgoing future-directed null-path after horizon crossing. As we discussed previously in \eqref{eq:clas_action_C}, this happens when $k_r\sim F^{-1}(r)$, which corresponds to the outgoing trajectory. This also establishes a clear link between the existence of tunneling, and hence of particle production, and the presence of a single pole in the classical action\footnote{Higher order poles lead to a purely real action cf. \cite{gelfand,hormander} for details.}.

\begin{figure}
	\centering
	\includegraphics[width=0.235\textwidth]{pathkh.pdf}
	\caption{The Killing horizon is displayed as a null-line that separates the normal (unshaded) from the trapped region (grayshaded). The tunneling path shows a positive energy particle that starts in the interior on a past-directed outgoing null path, crosses the horizon on a complex path (dashed), and scatters into a future-directed outgoing null path, once it has crossed the horizon. This is equivalent to the negative energy particle tunneling inwards.}
	\label{eins}
\end{figure}

Finally, we can define the tunneling rate as the ratio between the transmitted fraction and the incident wave \cite{sen15}
\begin{equation}\label{rate}
\Gamma=\frac{\|\phi_{>}\|^2}{\|\phi_{<}\|^2}\sim e^{-\frac{2}{\hbar}{\rm Im}(\mathcal{S}_0)},
\end{equation}
which vanishes along classical paths ${\rm Im}(\mathcal{S}_0)=0$, but it is finite otherwise. This is connected to thermodynamics by following \cite{har83,gib93}. Comparing the probability for a detector to absorb a particle $P_{\rm abs}$ with its probability to emit one $P_{\rm em}$ at a fixed energy $\Omega$ we get
\begin{equation}\label{bolli}
\Gamma\equiv\frac{P_{\rm abs}}{P_{\rm em}}=e^{-{\Omega}/{T_{\scaleto{\rm H}{2.5pt}}}},
\end{equation}
thus finding that the detector is in a thermal equilibrium at (horizon) temperature $T_{\rm H}$. The particle absorbed by the detector coincides with the particle that has crossed the horizon. Thus, whenever ${\rm Im}(\mathcal{S}_0)\propto \Omega$, we can read off a horizon temperature from the tunneling rate. However, the thermodynamics is only well-defined whenever the imaginary part is semi-positive definite \cite{gia20}
\begin{equation}\label{heb}
\mbox{Im}(\mathcal{S}_0)>0,
\end{equation}
otherwise the process leads to inconsistencies such as the violation of the probabilistic interpretation. In the subsequent analysis we will show the robustness of this criterion even in Lorentz-violating theories, which underlines the resilience of thermodynamic properties and the generality of the analysis in \cite{gia20}.

\subsection{Black holes in the tunneling picture}\label{sec:bh_rel}

As we discussed before, the classical action for an S-wave (spherically symmetric) in EFB coordinates is given by \eqref{eq:S0_relativ}, where all that remains to be done is to obtain the value of $k_r$. This is done by plugging \eqref{eq:S0_relativ} into \eqref{hje}, so arriving at the Hamilton--Jacobi equation in the EFB vacuum
\begin{equation}
2\Omega k_r-F(r)k^2_r=0.
\end{equation}
Tracing the particle that arrives to the asymptotic region of large radius down to the horizon, we find that its momentum diverges and that the particle thus travels on an approximately light-like trajectory, which allows us to neglect the mass term\footnote{The Hamilton--Jacobi equation is soluble even in the case of keeping the mass term. However, since we will perform later a near horizon approximation, this term will drop out anyways and does not add anything to the analysis but mere complication.} and solve for the radial null momentum $k_r(r)$. One finds two solutions, the first corresponding to the ingoing fields ($k_r=0$), and which is regular across the horizon; and the outgoing momentum \cite{hay09}
\begin{equation}
k_r(r)=\frac{2\Omega}{F(r)},
\end{equation}
which develops the aforementioned pole at the position of the horizon, thus requiring to be analytically continued in $r$ through $F(r)=0$. In our case we choose\footnote{The definition for a quantum process leading to a consistent contour integral, as well as thermodynamics, was derived in \cite{gia20} and is given by $\mbox{Im}(\mathcal{S}_0)>0$. This influences our choice of analytic continuation.} $F(r)\rightarrow F(r-i\varepsilon)$, with $\varepsilon \ll 1$. Alternatively, one can add a Feynman $i\varepsilon$-prescription in the Hamilton-Jacobi equation and solve directly for the imaginary part. As a result, $k_r(r)$ will contribute to Im$(\int$d$r k_r(r))$ such that after some manipulations (cf. \cite{gia20} for details)
\begin{equation}\label{ims0rel}
\mbox{Im}(\mathcal{S}_0)=\lim_{\varepsilon\to0}\mbox{Im}\left(\int_{r_1}^{r_2}\frac{\Omega\mbox{d}r}{ \kappa_{\textsc{kh}}
\left(r-r_{\textsc{kh}}-\frac{i\varepsilon}{\kappa_{\textsc{kh}}
}
\right)}\right)=\frac{\pi\Omega}{\kappa_{\textsc{kh}}},
\end{equation}
where we performed a near-horizon expansion $F(r)\cong\kappa_{\scaleto{\rm KH}{3pt}}(r-r_{\scaleto{\rm KH}{3pt}})+\mathcal{O}\left((r-r_{\scaleto{\rm KH}{3pt}})^2\right)$, and identified the surface gravity of the horizon as $\kappa_{\scaleto{\rm KH}{3pt}}=\tfrac{\partial F}{\partial r}|_{r=r_{\scaleto{\rm KH}{2.5pt}}}$.
Furthermore, we implemented the crossing path of the positive energy mode from the interior starting at $r_1<r_{\scaleto{\rm KH}{3pt}}$ and tunneling to the exterior region to $r_2>r_{\scaleto{\rm KH}{3pt}}$ \cite{par00}. Comparing the tunneling rate \eqref{rate} with the Boltzmann distribution \eqref{bolli} we find
\begin{equation}\label{reltemp}
\Gamma=e^{-{2\pi\Omega}/{\kappa_{\textsc{kh}}}}= e^{-{E}/{T_{\textsc{h}}}},
\end{equation}
yielding the standard result for the horizon temperature $T_H=\frac{\kappa_{\textsc{kh}}}{2\pi}$.

\section{Non-relativistic tunneling} \label{sec:non_re_tunneling}
After summarizing the tunneling approach in relativistic theories, let us now turn to the main topic of this work. In the following, we will extend the tunneling picture to the case of space-times endowed with a UH in Lorentz violating theories. In particular, we will focus on EA gravity, which for our purposes here also includes the low energy limit of Ho\v rava--Lifshitz Gravity. In what follows, we investigate tunneling across UHs, and hereinafter verify the principles formulated in \cite{gia20}, and revisited in the previous section, for Lorentz breaking theories. Their persistence even in absence of Lorentz symmetry fortifies the universality of horizon thermodynamics.

\subsection{Einstein--Aether Gravity}
We consider matter actions coupled to EA gravity, with action
\begin{align}\label{eq:EA-action}
    I_{\rm EA}=\frac{1}{16\pi G}\int \mbox{d}^4x \sqrt{-g} \left(R-{\cal L}_{\rm U} + \lambda(U_\mu U^\mu+1)\right),
\end{align}
where $R$ is the Ricci scalar, $\lambda$ is a Lagrange multiplier implementing a unit norm condition for the aether field $U^\mu$, and ${\cal L}_{\rm U}=K^{\alpha\beta}_{\mu\nu}\nabla_\alpha U^\mu \nabla_\beta U^\nu$, with $K^{\alpha\beta}_{\mu\nu}=c_1 g^{\alpha\beta}g_{\mu\nu}+c_2 \delta^\alpha_\mu \delta^\beta_\nu +c_3 \delta^\alpha_\nu \delta^\beta_\mu + c_4 U^\alpha U^\beta g_{\mu\nu}$ and couplings $c_i\in\mathbb{R}$.

This action is invariant under the FDiff transformations \eqref{eq:FDiff}, which then dictate the form of the matter action. Hereinafter we take the simplest case of a Lifshitz scalar field\footnote{Here we set a possible mass term to vanish. As in the relativistic tunneling, its presence is irrelevant and thus we can neglect it without loss of generality.}
\begin{align}\label{eq:lifshitz}
I_{\rm m}=-\frac12\int {\rm d}^4x \sqrt{-g}\left[\partial^\mu\phi\partial_\mu\phi+\sum_{z=2}^Z\frac{\alpha_{2z}}{\Lambda^{2z-2}}\phi(-\Delta)^z\phi\right],
\end{align}
whose corresponding equation of motion is thus
\begin{align}\label{eq:lifhistz_eom}
    \square \phi -\sum_{z=2}^Z\frac{\alpha_{2z}}{\Lambda^{2z-2}}(-\Delta)^z\phi=0.
\end{align}

Here, all the $\alpha_z$ are dimensionless, and we choose to normalize $\alpha_{2Z}=1$. The scale $\Lambda$ sets the momentum scale at which Lorentz violations become relevant. In 3+1 dimensions, usually $N=3$ is chosen so to enforce power-counting renormalizability of the gravitational action in the case of Ho\v rava--Lifschitz (HL) gravity \cite{Pospelov:2010mp,Anselmi:2007ri,Barvinsky:2015kil}. The d'Alembert as well as the spatial Laplace operators are given by $\square=g^{\mu\nu}\nabla_\mu\nabla_\nu$ and $\Delta=\gamma^{\mu\nu}\nabla_\mu\nabla_\nu$, with $\gamma_{\mu\nu}=g_{\mu\nu}-U_\mu U_\nu$ the metric induced on the leafs orthogonal to the aether. We further restrain ourselves to space-times with UHs, for which it is a necessary condition that the aether is hypersurface orthogonal, taking the form \eqref{eq:aether_ortho}. Examples of vacuum solutions of this kind can be found analytically \cite{berg12} and numerically \cite{Barausse:2011pu} for generic regions of the parameter space of EA gravity. All of them equally correspond to vacuum solutions of Khronometric gravity \cite{Blas:2010hb}. Finally, we will also require the existence of a time-like Killing vector $\chi^\mu$ -- i.e.~both metric and aether are supposed to be Lie dragged by this field --, from which a notion of a conserved Killing energy for a particle of four-momentum $k^\mu$, namely $\Omega=-(\chi\cdot k)$,  can be derived.

Note that the equation of motion \eqref{eq:lifhistz_eom} generically contains higher spatial and time derivatives, since $\gamma^{00}$ and $\gamma^{0i}$ are generically non-vanishing. This fact is quite problematic, because it leads to Ostrogradsky ghosts, signaling a classical runaway instability and a loss of unitarity in the quantum theory. Only when the time direction is identified with the integral lines of $U^\mu$, the equation of motion remains second order in time derivatives. This selects a preferred time direction that every motion must follow within this space-time, and thus a universal preferred frame, up to FDiff transformations of course. In this frame, the scalar field $\phi$ inherits a modified dispersion relation which can be read from \eqref{eq:lifhistz_eom}.

Given a particle four-momentum $k_\mu\phi=-i(\partial_\mu\phi)$, it can always be decomposed in the aether frame as
\begin{align}\label{eq:k_vector}
    k^\mu=\omega U^\mu+k_\rho S^\mu,
\end{align}
where $S^\mu$ is the space-like unit vector orthogonal to $U^\mu$.
Here $\omega$ and $k_\rho$ are the energy and spatial momentum in the aether frame
\begin{align}
    \omega=-(U\cdot k),\quad k_\rho =(S\cdot k).
\end{align}
Hence, in the preferred frame \eqref{eq:lifhistz_eom} implies the dispersion relation
\begin{align}
   \omega^2=k_\rho^2+\sum_{z=1}^Z \alpha_{2z}\frac{k_\rho^{2z}}{\Lambda^{2z-2}}.
\end{align}
Note however that these aether frame quantities are \emph{not constant} along the motion. As said before, it is the Killing energy to be conserved along free motions under the sole influence of metric and aether. Given the above definitions we decompose it as
\begin{align}
  \Omega=  (\chi\cdot k)=\omega(U\cdot \chi)-k_\rho (S\cdot \chi),
  \label{eq:decOm}
\end{align}

\subsection{Spherically symmetric black holes}

In the following we consider arbitrary spherically symmetric and static solutions to the equations of motion obtained from \eqref{eq:EA-action}, written in the form \eqref{efb} with $\xi\equiv v$
\begin{equation}
\mbox{d}s^2=-F(r)\mbox{d}v^2+2B(r)\mbox{d}v\mbox{d}r+r^2\mbox{d}^2S_2,
\end{equation}
and endowed with an aether $U^\mu$ of the form
\begin{equation}
U_\mu \mbox{d}x^\mu=-\frac{1+F(r)A^2(r)}{2 A(r)}\mbox{d}v+B(r)A(r)\mbox{d}r,\quad S_\mu \mbox{d}x^\mu=\frac{1-F(r)A^2(r)}{2 A(r)} \mbox{d}v+B(r)A(r)\mbox{d}r,
\end{equation}
satisfying $U_\mu U^\mu=-1$. As said above, $S^\mu$ is the vector orthogonal to $U^\mu$, which satisfies $S_\mu S^\mu=1$, and it is chosen to be inwards pointing. The functions $F(r)$, $B(r)$, and $A(r)$ must be determined case by case by solving explicitly the equations of motion derived from \eqref{eq:EA-action}. Hereinafter we keep the functions arbitrary, albeit demanding the existence of a Killing as well as a universal horizon. The former is localized by the condition $F(r_{\scaleto{\rm KH}{3pt}})=0$, while the latter is given by
\begin{align}
\left.(U\cdot\chi)\right|_{\scaleto{\rm UH}{3pt}}=-\frac{1+F(r_{\scaleto{\rm UH}{3pt}})A^2(r_{\scaleto{\rm UH}{3pt}})}{2 A(r_{\scaleto{\rm UH}{3pt}})}=0,
\end{align}
with $\chi^\mu=(1,0,0,0)$ the time-like Killing vector of the metric. Incidentally, note that  the lapse $N$ of the foliation described by the surfaces orthogonal to the aether field (cf. below), is $N=-(U\cdot\chi)$ and so it changes sign crossing the UH \cite{del22}. Finally, we also demand asymptotic flatness, which implies that $U^\mu$ becomes parallel to $\chi^\mu$ at large radius, thus fixing the asymptotics of $F(r)$ and $A(r)$ to $F(r)=A(r)=1$ for $r\rightarrow \infty$.

Since the aether defines a physical foliation, it is convenient to change coordinates into the preferred system, that allows us to align the notion of time with the evolution direction of the aether, such that the spatial vector will always be tangent to the foliation. The following transformation allows to change into this system, given by $\{\tau,\rho,\vartheta,\varphi\}$ \cite{cro14}
\begin{equation}\label{trafo}
\tau=v+\int\frac{U_r}{U_v}\mbox{d}r\quad\mbox{and}\quad\rho=\int\frac{S_r}{S_v}\mbox{d}r+v.
\end{equation}

There is an important point to take into account here, since the integral defining $\tau$ takes the form
\begin{equation}
\int\frac{U_r}{U_v}\mbox{d}r=\int \frac{-2B(r)A^2(r)}{1+F(r)A^2(r)}\mbox{d}r,
\end{equation}
it actually diverges at the UH and in its proximity behaves as
\begin{align}
    \int\frac{U_r}{U_v}\mbox{d}r=\frac{B(r_{\scaleto{\rm UH}{3pt}})A(r_{\scaleto{\rm UH}{3pt}}) }{\left. \partial_r N\right|_{\scaleto{\rm UH}{3pt}} } \ln\left|r-r_{\scaleto{\rm UH}{3pt}}\right|+{\cal O}(r-r_{\scaleto{\rm UH}{3pt}}),
\end{align}
where we have expanded the lapse as $N=\left. \partial_r N\right|_{\scaleto{\rm UH}{3pt}} (r-r_{\scaleto{\rm UH}{3pt}})$, and taken into account that $A(r)$ and $B(r)$ are regular everywhere. In these coordinates, the position of the UH thus corresponds to $\tau\rightarrow \infty$, with the foliation exhibiting an infinite accumulation there. This matches the behavior of the light-cone coordinate $u$ in standard approaches to Hawking radiation (cf. \cite{Jacobson:2003vx} and references within), and hence serves as a smoking gun to expect a similar property from the UH. Moreover, we can also note that the logarithm requires an analytic continuation through $r=r_{\scaleto{\rm UH}{3pt}}$. This will introduce an imaginary part that will eventually resurface in the solution to the Hamilton--Jacobi equation for the field in \eqref{eq:lifshitz}, in similar manner to the relativistic case discussed in section \ref{sec:bh_rel}.

In foliation adapted coordinates the metric takes the ADM form
\begin{equation}\label{metricpc}
\mbox{d}s^2=-(N^2-N^i   N_i)\mbox{d}\tau^2+2N_\rho \mbox{d}\tau\mbox{d}\rho+\gamma_{ij}\mbox{d}x^i \mbox{d}x^j,
\end{equation}
where the Latin indices now run over the spatial directions $\{\rho, \vartheta,\varphi\}$ only. The lapse $N$, shift vector $N^i$, and spatial metric $\gamma_{ij}$, are given by
\begin{align}\label{eq:lapse}
    &N=\frac{1+F(r)A^2(r)}{2 A(r)},\\
    &N^i=0,\\
    &\gamma_{ij}\mbox{d}x^i \mbox{d}x^j=\frac{(1-F(r)A^2(r))^2}{4 A^2(r)}\mbox{d}\rho^2+r^2\mbox{d}^2S_2,
\end{align}
where $r$ must be understood as $r(\tau,\rho)$ at all times. Notice that the UH is now determined by the condition $g_{\tau\tau}=0$, which again points towards a full analogy with the relativistic case. This criterion is fully consistent with the covariant criterion given in \cite{di09}.

We thus proceed as in the relativistic case, by constructing a solution for the Lifshitz scalar field $\phi$ by means of the WKB approximation, obtaining the corresponding  Hamilton--Jacobi equation. To this aim, we need the inverse metric
\begin{equation}\label{eq:inverse_metric}
g^{-1}=\mbox{diag}\left(-\frac{1}{N^2},\gamma^{\rho\rho},\frac{1}{r^2},\frac{1}{r^2\sin^2(\vartheta)}\right).
\end{equation}
All divergences when approaching the UH are encoded in the behavior of the lapse, so mere observation is enough to spot that at the UH, only $g^{00}$ blows up, while the other components remain finite. Note also that close to the Killing horizon all coefficients are regular. This will be important in stressing the difference between the Killing and the universal horizon in the setting at hand.

\subsection{WKB approximation and Hamilton--Jacobi equation}
In order to derive the Hamilton--Jacobi equation we first need to formulate an ansatz for our field through the WKB approximation \eqref{wkb}.
 Again, in the semiclassical limit, the dominant contribution is given by the classical action $\mathcal{S}_0$.  Using the ansatz\footnote{From now on, we work in natural units such that although formally $\hbar$ plays the role of a smallness parameter that is later identified with the physical $\hbar$, we will set it to one.}
\begin{equation}\label{s0dec}
\mathcal{S}_0=\int\mbox{d}x^a\partial_a\mathcal{S}_0=-\int \omega\mbox{d}\tau+\int k_\rho\mbox{d}\rho,
\end{equation}
where we have assumed spherically symmetric waves. Using it, we can obtain the corresponding leading order of the Hamilton--Jacobi equation by inserting the WKB approximated field in \eqref{eq:lifhistz_eom} to get
\begin{equation}\label{hamjac}
-\frac{\omega^2}{N^2}+\gamma^{\rho\rho}k_\rho^2+\frac{(\gamma^{\rho\rho})^Z}{\Lambda^{2Z-2}}\left(k_\rho^{2Z}+G(k_\rho,\nabla k_\rho;\Lambda)\right)=0.
\end{equation}
The operator $G(k_\rho,\nabla k_\rho;\Lambda)$ contains subleading terms that depend on a combination of $k_\rho$ and derivatives thereof, but which are suppressed at the UH. Since the leading contribution is given by the highest power, we singled this term out.

To understand better the behavior of the quantum fields across the horizon, we investigate the structure of \eqref{hamjac} by following the discussion in \cite{mi15,he21}. We distinguish two regimes -- the \emph{soft regime}, corresponding to $\Lambda\to\infty$, and the \emph{hard regime} $\Lambda\to 0$. The former corresponds to rays that cross the universal horizon with finite momentum, while in the latter they have been red-shifted by climbing up the gravitational well, such that their momentum diverges when traced back to the universal horizon, lingering there eternally. A detailed discussion on the behavior of the modes can be found in \cite{mi15,he21}, and it is summarized in Figure \ref{pmodes}.
\begin{figure}[]
	\centering
	\includegraphics[width=0.4\textwidth]{modes.pdf}
	\caption{Structure of the modes. Close to the UH and outside of it we have four modes -- ingoing $\phi_{\rm in}$, outgoing $\phi_{\rm out}$, and the two tails of the trapped one $\phi_{\rm tr}$. However, the latter is indistinguishable from the former in the region $r\sim r_{\rm UH}$. Note that the mode $\phi_{\rm p}$ is the Hawking partner to the red mode. The UH is denoted by the purple line, while the Killing horizon is showed as dotted line. The zig-zag line on the left depicts the singularity.} \label{pmodes}
\end{figure}

At any moment, outside the universal horizon but still nearby, we have four solutions. The ingoing blue mode $\phi_{\rm in}$ (soft), the outgoing red mode $\phi_{\rm out}$ (hard), and the \emph{trapped} orange mode $\phi_{\rm tr}$, which departs the UH as hard, but eventually turns back, crossing it as a soft mode. Although this mode is clearly different form the other two, it mimics their behaviors when we zoom close to the UH. In order to single out the orange mode, we would need to trace it across the region enclosed between the two horizons. For our purposes here, it is thus enough to talk about soft and hard modes, since the trapped one degenerates locally with the ingoing and outgoing rays.

In the interior of the UH we find a similar structure, although with the roles of $\phi_{\rm tr}$ and $\phi_{\rm in}$ exchanged. While the former remains soft, the latter becomes hard when approaching the UH from its interior. In this region we also find a purely hard mode $\phi_{\rm p}$, corresponding to the CP reversed --- and thus also T reversed, provided that CPT is conserved in our framework --- partner of $\phi_{\rm out}$. This mode is equally diverging as the latter in the UH, but travels from the singularity towards the horizon\footnote{Note however that the direction of the evolution of the rays is tied to the notion of time advance as measured by an observer. The description depicted here corresponds to an observer sitting \emph{outside} the UH. A different observer sitting on its interior would instead measure a reversed time flow -- because the sign of their lapse is flipped, cf. \cite{del22} -- and would assign a reversed evolution to the modes.}.

Once this is settled, and in order to understand the contribution of all these modes to the tunnelling probability, we find the ``soft" and ``hard" limits of the Hamilton--Jacobi equation to be
\begin{eqnarray}
-\frac{\omega^2}{N^2}+\gamma^{\rho\rho}k_\rho^2+m^2+\mathcal{O}\left(\frac{1}{\Lambda}\right)&=&0,\label{shje}\qquad\mbox{Soft}\\
\frac{1}{\Lambda^{2Z-2}}\left(k_\rho^{2Z}+G(k_\rho,\nabla k_\rho;\Lambda)+\mathcal{O}(\Lambda)\right)&=&0. \qquad\mbox{Hard}\label{hhje}
\end{eqnarray}

As previously mentioned however, and differently from the Killing energy $\Omega$, neither $\omega$ nor $k_\rho$ are constants of motion, which obstructs the efforts of solving the previous equations explicitly. However, we can make use of the relation Eq.~\eqref{eq:decOm} among these quantities to write
\begin{equation}\label{kleinomega}
\omega=-\frac{k_\rho (S\cdot\chi)+\Omega}{(U\cdot\chi)}.
\end{equation}

We also need to solve $k_\rho$ in terms of $\Omega$ and its associated spatial momentum $k_r=(k\cdot \zeta)$, where $\zeta^\mu$ is the vector orthogonal to $\chi^\mu$ and satisfying $\zeta^2=-\chi^2$, taken inwards pointing. Contracting now \eqref{eq:k_vector} with $\zeta^\mu$ and solving for $k_\rho$ we get
\begin{align}\label{krho}
    k_\rho=-\Theta \left(\Omega+\frac{(U\cdot \chi)}{(U\cdot \zeta)}k_r\right),
\end{align}
where we have defined the position dependent angle to be
\begin{align}
    \Theta=\left[(S\cdot\chi)-\frac{(S\cdot\zeta)(U\cdot\chi)}{(U\cdot \zeta)}\right]^{-1}.
\end{align}
Asymptotically $\Theta$ vanishes, while close to the UH, we find instead
\begin{align}\label{eq:theta_limit}
    \left.\Theta\right|_{\rm UH}=(S\cdot\chi)^{-1}.
\end{align}
By substituting these quantities into the Hamilton--Jacobi equation, we have all the required tools to study the contribution of all the different modes to the tunneling probabilities.

\subsubsection{Soft Hamilton--Jacobi equation}
Let us start by considering the soft Hamilton--Jacobi equation \eqref{shje} with the substitutions \eqref{kleinomega} and \eqref{krho}. This introduces additional dependencies on the lapse function, as we can see by noting that $(U\cdot\chi)=-N$. Additionally, we expect these modes to cross the universal horizon with finite momentum. Taking this into account and working close to the universal horizon, so that $N\sim \left. \partial_r N\right|_{\scaleto{\rm UH}{3pt}}(r-r_{\scaleto{\rm UH}{3pt}})$, the leading term in \eqref{shje} becomes
\begin{align}
    \frac{\left(\Omega+k_\rho (S\cdot\chi)\right)^2}{N^4}+{\cal O}\left(\frac{1}{N^2}\right)=0,
\end{align}
with solution
\begin{align}
k_\rho=-\frac{\Omega}{(S\cdot\chi)}.
\end{align}
In terms of $k_r$ this turns out to be
\begin{align}
    k_r^{\rm s}=\Omega\left[1-\Theta (S\cdot \chi)\right]\frac{(U\cdot \zeta)}{(U\cdot \chi)}=-\Omega\left[1-\Theta (S\cdot \chi)\right]\frac{(U\cdot \zeta)}{N}.
\end{align}
Here we have indicated the soft (and later the hard) modes with the superscript `s' (and `h' respectively).

What is left is to evaluate $\mathcal{S}_0$ to obtain its imaginary part, if any. In order to do this, we integrate the solution above through a contour integral connecting the exterior and interior of the universal horizon. Since the soft modes are regular along the UH, they have support along $r=r_{\scaleto{\rm UH}{3pt}}$ and the integration is straightforward
\begin{align}
{\rm Im}\left(\mathcal{S}^{\rm s}_0\right)&={\rm Im}\left(\oint_c k^{\rm s}_{r}(r)\mbox{d}r\right)=-\frac{\Omega \left[1-\left.\Theta\right|_{\scaleto{\rm UH}{3pt}} (S\cdot \chi)_{\scaleto{\rm UH}{3pt}}\right](U\cdot \zeta)_{\scaleto{\rm UH}{3pt}}}{\left. \partial_r N\right|_{\scaleto{\rm UH}{3pt}}}{\rm Im}\left(\oint_c\frac{\mbox{d}r}{r-r_{\scaleto{\rm UH}{3pt}}}\right),
\end{align}
where we have used $(U\cdot \chi)=-N$. Using again the Sokhotski--Plemelj theorem (cf. footnote \ref{ft:footnote3}) to resolve the pole in the lapse, this evaluates to
\begin{align}
   {\rm Im}\left(\mathcal{S}^{\rm s}_0\right) =\frac{\pi \Omega \left[1-\left.\Theta\right|_{\scaleto{\rm UH}{3pt}} (S\cdot \chi)_{\scaleto{\rm UH}{3pt}}\right](U\cdot \zeta)_{\scaleto{\rm UH}{3pt}}}{\left. \partial_r N\right|_{\scaleto{\rm UH}{3pt}}}=0,
\end{align}
which vanishes due to \eqref{eq:theta_limit}. This implies that the modes following these trajectories do not escape the gravitational well. This is in agreement with their description as purely incoming, and thus classically allowed, through the UH.

\subsubsection{Hard Hamilton--Jacobi equation}
Now we take the version \eqref{hhje} of the Hamilton--Jacobi equation and proceed similarly to the previous case. Note that although it does not look dynamical, it yields a dynamical equation for the free-falling observer once the relation \eqref{krho} is imposed. The obvious solution to \eqref{hhje} is $k_\rho=0$, which in turns implies
\begin{equation}\label{hardmodes}
    k_\rho=0\quad \longrightarrow\quad  k^{\rm h}_r=- \Omega\times\frac{(U\cdot \zeta)}{(U\cdot\chi)}.
\end{equation}

As before, we again integrate $k^{\rm h}_{r}$ close to the universal horizon through a contour integral (cf. footnote \ref{ft:footnote3}). However, this time we find that the hard modes do not have support through the UH, so the tunneling path has to be made out of combining an interior hard mode -- either $\phi_{\rm p}$ or $\phi_{\rm in}$ -- and an exterior one -- $\phi_{\rm out}$ or $\phi_{\rm tr}$. In practice, however, this simply corresponds to analytically continuing one of the modes through the UH, while ensuring that $\partial_r N$ is continuous -- cf. \cite{del22}. Moreover, by combining the modes in this way, we find a precise reproduction of Hawking's cartoon picture of gravitational tunneling.

Performing the integral, this time we find a non-trivial contribution to the imaginary part
\begin{equation}\label{eq:result_hard}
{\rm Im}\left(\mathcal{S}^{\rm h}_0\right)={\rm Im}\left(\oint k^{\rm h}_{r}(r)\mbox{d}r\right)=\Omega \times \frac{ (U\cdot \zeta)_{\scaleto{\rm UH}{3pt}}}{\left. \partial_r N\right|_{\scaleto{\rm UH}{3pt}}}{\rm Im}\left(\oint\frac{\mbox{d}r}{r-r_{\scaleto{\rm UH}{3pt}}}\right)=\pi \Omega\times \frac{(U\cdot \zeta)_{\scaleto{\rm UH}{3pt}}}{\left. \partial_r N\right|_{\scaleto{\rm UH}{3pt}}}=\frac{\pi \Omega  \sqrt{-F(r_{\scaleto{\rm UH}{3pt}})}}{\left. \partial_r N\right|_{\scaleto{\rm UH}{3pt}}},
\end{equation}
where we have used \eqref{eq:theta_limit} and evaluated $(U\cdot\zeta)_{\scaleto{\rm UH}{3pt}}= \sqrt{-F(r_{\scaleto{\rm UH}{3pt}})}$.

In contrast to the ingoing mode, we obtain here instead a non-vanishing imaginary contribution to $\mathcal{S}^{\rm h}_0$. This highlights the fact that the outgoing mode corresponds to trajectories whose momentum diverges at the universal horizon, developing a pole which demands an analytical continuation. Classically, these are unable to cross the surface at $r=r_{\scaleto{\rm UH}{3pt}}$, but quantum mechanical tunneling makes this possible.

Finally, as in the relativistic case, we can obtain the temperature of the radiation composed by these modes through a comparison with the Boltzmann distribution
\begin{equation}
\Gamma=\frac{\|\phi_>\|}{\|\phi_<\|}\propto e^{-2{\rm Im}(\mathcal{S}_0)}\equiv e^{-\Omega/ {T_{\scaleto{\rm H}{2.5pt}}}}
\end{equation}
so that
\begin{equation}\label{eq:temp}
T_{\rm H}=\frac{\left. \partial_r N\right|_{\scaleto{\rm UH}{3pt}}}{2\pi\sqrt{-F(r_{\scaleto{\rm UH}{3pt}})}}.
\end{equation}

Note that our result for
${\rm Im}(\mathcal{S}_0)$ also fulfills the property \eqref{heb} which shows that this setup admits well-defined thermodynamics. Note in particular that in the case of UHs discussed here, the existence of an analytic continuation of $N^{-1}$ through $r=r_{\scaleto{\rm UH}{3pt}}$ requires $\left.\partial_r N\right|_{\scaleto{\rm UH}{3pt}}$ to be continuous across the UH, thus connecting the consistency of the tunneling approach with the results of \cite{del22}.

Finally, we just need to evaluate the value of $\left.\partial_r N\right|_{\scaleto{\rm UH}{3pt}}$. In order to do that, we construct the acceleration of the aether $a^\mu=U^\nu\nabla_\nu U^\mu$, and contract it with the Killing vector, getting
\begin{align}
    (\chi\cdot a)=U^\nu \partial_\nu (U\cdot\chi)=-U^r \partial_r N,
\end{align}
where we have used that the metric and aether are static and that $\chi^\mu$ is a Killing vector. Evaluating it close to the UH we find
\begin{align}
    - \frac{\left.\partial_r N\right|_{\scaleto{\rm UH}{3pt}}}{\sqrt{-F(r_{\scaleto{\rm UH}{3pt}})} } = (\chi\cdot a)_{\scaleto{\rm UH}{3pt}} = -2 \kappa_{\scaleto{\rm UH}{3pt}},
\end{align}
where $\kappa_{\scaleto{\rm UH}{3pt}}$ is the surface gravity evaluated at the UH. Substituting this onto \eqref{eq:temp} we finally arrive at the horizon temperature
\begin{align}
    T_{\rm H}=\frac{\kappa_{\scaleto{\rm UH}{3pt}}}{\pi}.
\end{align}
This result agrees with that obtained in \cite{del22} by other methods.

\subsubsection{Radiation from the Killing horizon?}\label{sec:killing_radiation}
The analysis performed in this section hints towards the conclusion that the UH controls the thermodynamic properties of the solution, rather than the Killing horizon. To confirm this, we can observe what happens at the latter. In this region of space-time, neither the elements of the inverse metric \eqref{eq:inverse_metric}, nor the relations \eqref{kleinomega} and \eqref{krho} between aether frame and EFB frame quantities diverge at the Killing horizon. Due to this, the momentum of the modes will also be regular when crossing the Killing horizon, signaling that it is not a causal barrier anymore and can be exited in a finite -- perhaps long -- time. Thus, the integral $\oint k_r(r)$d$r$ along a path crossing the Killing horizon will be strictly real for all the modes, no matter their character. Since no mode develops a pole at the Killing horizon, there is no need for an analytic continuation. This shows that the UH is the sole responsible for the thermodynamical properties of the system, being also the only true causal boundary within the spacetime. The role of the Killing horizon must be, at most, to introduce a graybody factor that might distort the shape of the distribution measured by an observer sitting at large radius, compared to the one emitted close to the universal horizon. However, this analysis is beyond the scope of this work.

\section{Time reparametrization invariance and UV sensitivity of the temperature}\label{sec:sync_factor}

In the previous section we obtained the temperature of the radiation emitted by the UH via the Hamilton--Jacobi method and the ansatz \eqref{s0dec}. The reader might notice that in writing this ansatz we chose a particular foliation time $\tau$, which identifies the lapse function of the foliation by $N=-(U\cdot\chi)$. However, this construction is far from unique. Indeed, as discussed in the introduction, EA configurations are invariant under time reparametrizations
\begin{align}\label{eq:time_rep}
    \tau\rightarrow \hat\tau(\tau),
\end{align}
which preserve the form of $U^\mu$, but not the choice of the lapse function. Under a transformation \eqref{eq:time_rep} the lapse function transforms as
\begin{align}\label{eq:lapse_relation}
    \hat{N}=\frac{{\rm d}\tau}{{\rm d}\hat \tau} N,
\end{align}
which has a strong effect on the value of the temperature \eqref{eq:temp}. This is not a surprise of course. Taking into account that the temperature is a measure of the mean energy per particle, and that time and energy are conjugated variables, it is clear to see how a redefinition of the former affects the latter. What thus fixes the right choice of foliation time? In order to answer this question, we must take into account that the rays that arrive to the asymptotic region of large radius have climbed the gravitational well by following free-falling trajectories, understood here as those ruled by their equations of motion. Thus, we synchronize the clocks in the asymptotic region with that dictated by the rays. This is tantamount to select the proper vacuum state, and the right notion of energy for the field.

In order to make this connection explicit, we aim to compute the four-velocity of the ray close to the UH, following \cite{cro14}. We thus recall \eqref{kleinomega} and solve it together with the dispersion relation $\omega(k_\rho)$ to obtain the value of $k_\rho$. For this, it suffices to consider only the UV limit of the dispersion relation, since the modes that arrive at the asymptotic region correspond to hard blue-shifted modes close to the UH. Hence, we take $\omega= k_\rho^{Z}/\Lambda^{Z-1}$, and perform a WKB expansion for large momentum $k_\rho\sim \Lambda$. At leading order we get
\begin{align}\label{lambdaentw}
    k_\rho^{Z-1} \Lambda^{1-Z}= \frac{1-F(r)A^2(r)}{1+F(r)A^2(r)}.
\end{align}

From here we compute the group velocity of the field as
\begin{align}\label{ck}
    c(k_\rho)=\frac{{\rm d}\omega}{{\rm d} k_\rho}=Z\times \frac{1-F(r)A^2(r)}{1+F(r)A^2(r)},
\end{align}
and we use it to build the four-velocity
\begin{align}
    V^\mu=U^\mu + c(k_\rho) S^\mu.
\end{align}

The vector $V^\mu$ is then employed to extract the trajectory of the field. In particular we are interested in ${\rm d} v/{\rm d}r=V^0/V^1$, which close to the UH reads
\begin{align}\label{strahl}
\left.\frac{{\rm d} v}{{\rm d}r}\right|_{\scaleto{\rm UH}{3pt}}=-\frac{Z}{Z-1} \frac{B(r_{\scaleto{\rm UH}{3pt}})}{\sqrt{-F(r_{\scaleto{\rm UH}{3pt}})}}\frac{1}{N}+{\cal O}(N).
\end{align}

We can compare this with the relation ${\rm d}\tau/{\rm d}r$, which can be obtained by noting that close to the UH, $U_0\sim 0$ so that
\begin{align}
    -N {\rm d}\tau= \left. U_\mu {\rm d}x^\mu\right|_{\scaleto{\rm UH}{3pt}}=U_r {\rm d}r ,
\end{align}
and therefore
\begin{align}
   \left. \frac{{\rm d}\tau}{{\rm d}r}\right|_{\scaleto{\rm UH}{3pt}}=-\frac{U_r}{N}=-\frac{B(r_{\scaleto{\rm UH}{3pt}})}{\sqrt{-F(r_{\scaleto{\rm UH}{3pt}})}}\frac 1 N.
\end{align}

Thus, we observe that the clock of the rays aligns with the foliation clock close to the UH, with a proportionality factor
\begin{align}
    \frac{{\rm d}\tau}{{\rm d}v}=\frac{Z-1}{Z},
\end{align}
which in turn implies that the ray sees a lapse
\begin{align}
    \hat N=\frac{Z-1}{Z} N,
\end{align}
close to the UH. Propagating this synchronization factor throughout our computation, we thus find
\begin{align}\label{eq:true_temp}
    \hat{T}_{\rm H}=\frac{Z-1}{Z}\frac{\kappa_{\scaleto{\rm UH}{3pt}}}{\pi},
\end{align}
which again agrees with previous results in the literature \cite{del22,he21}.

\section{Conclusions and discussion}\label{sec:conclusions}

In this work, we have studied the application of the gravitational tunneling method to derive the distribution of radiation emitted by a universal horizon in Einstein-Aether and Ho\v rava--Lifshitz Gravity. We have done it by focusing on spherically symmetric black hole solutions, endowed with a gravitating scalar field satisfying the symmetries of the background and with an anisotropic scaling between time and spatial directions. Our findings suggest that the collection of knowledge built from understanding the role of quantum fields close to the event horizon of a general relativistic black hole translates almost straightforwardly to the case under study here.

In particular, we have found that there exist non-vanishing tunneling probabilities associated to classically forbidden trajectories escaping the interior of the UH, which acts here as a universal causal boundary, a role which is again akin to the one of the event horizon in General Relativity. The classical action for these trajectories, which peel off the UH when traced back in time from the asymptotic region where our observer is placed, develops a non-vanishing imaginary part, inherited from a pole in the spatial momentum, which leads to a well-defined temperature once the tunneling probability is compared to a Boltzmann distribution, given by \eqref{eq:true_temp}. This can be interpreted as the existence of Hawking radiation emitted by the UH in a thermal ensemble made of quanta of the scalar field. The advantage of the tunneling method here with respect to other approaches is apparent. Through analyzing the behavior of the modes of the scalar fields close to the UH we have been able to obtain the value of the temperature for all possible static and spherically symmetric black hole solutions endowed with a UH. This is in agreement with recent works where this issue is analyzed from other directions.

An important property to highlight here is that the temperature shows a dependence on the dispersion relation of the scalar field through the pre-factor $(Z-1)/Z$, where $2Z$ is the highest power of the spatial preferred momentum in the square of the dispersion relation -- i.e. $\omega^2\sim k_\rho^{2Z}$ at large momentum. This factor has an important significance for model building and for the consistency of the model discussed here. Unless all matter species coupled to gravity share the same UV scaling in their dispersion relation, the temperatures of the different thermal ensembles emitted by the black hole will differ. In such a case, it is not difficult to envision how to construct a perpetuum mobile of the second kind, thus violating the laws of thermodynamics \cite{Dubovsky:2006vk}. Therefore, our results here suggest that a universal high energy behavior of the dispersion relation should exist in a consistent theory. This is in line with findings and discussions on the renormalization group flow properties of models coupling Ho\v rava--Lifshitz Gravity to matter actions, where indeed a universal UV scaling is present and controlled by the gravitational action \cite{Pospelov:2010mp}.

Some open questions linger though. In particular, and although we have determined the spectrum of radiation emitted by the UH, there is still a lack of profound understanding on the role of the Killing horizon in the process. Although no radiation is emitted by this surface in the present setting, as we have discussed in subsection \ref{sec:killing_radiation}, it might still have an important effect on the distribution of modes measured by an observer sitting at large radius, in the exterior of all horizons. Different rays of different energies will climb the gravitational well at a different rate, due to the momentum dependence of the dispersion relation and group velocity. In particular, rays with small momentum will linger for a long time close to the Killing horizon before being able to escape from it, accumulating energy and shifting their frequency. One could then expect a distortion of the spectrum measured by an asymptotic observer when compared to the emitted one
\begin{align}
    \rho_{\rm measured}(\Omega)=\int \frac{{\rm d}\Omega'}{2\pi}\  {\cal G}(\Omega,\Omega')\times \rho_{\rm emitted}(\Omega'),
\end{align}
where the graybody kernel ${\cal G}(\Omega?,\Omega)$ encodes the effects introduced by the Killing horizon. Its computation is beyond the scope of this work, but is poses an exciting problem for future research. Is the measured spectrum still the one of a thermal ensemble with temperature $\hat T_{\rm UH}$? Is it shifted so that it looks like a thermal emission with $T_{\rm H}$ instead, thus washing off the presence of Lorentz violations? Or is it something in between?

Finally, it is worth to mention that all of our results here are obtained by assuming a fixed background with no back-reaction onto the geometry due to the presence of the matter fields. In particular, the black holes discussed here are solutions only to the \emph{low energy} Lagrangian of HL Gravity -- usually denoted by ${\cal L}_2$ -- and described here in terms of EA gravity, which serves as an effective field theory description of Lorentz violations at the second derivative level.
However, the more realistic -- and perhaps UV complete \cite{Barvinsky:2017kob,Barvinsky:2019rwn,Barvinsky:2021ubv} -- action of HL Gravity contains also terms with four and six derivatives -- in $3+1$ dimensions. Taking into account the fact that the modes investigated here, and ultimately responsible of the Hawking radiation, are blue-shifted to high energies, one could wonder if, in a general case, they would trigger perturbations of the geometry with comparable momentum, thus exiting the applicability of the low energy theory ${\cal L}_2$. This would demand the presence of the higher order terms in order to account for their dynamics. Stability of the UH, or even of the full space-time, in the presence of these perturbations is so far an open problem in the field; and one that we should keep in mind when discussing scenarios like the one in this work.

\acknowledgements
M. H-V. wants to thank the APP department at SISSA for their hospitality during the final stages of this work. The work of F. D. P., S. L., and M. S. has been supported by the Italian Ministry of
Education and Scientific Research (MIUR) under the grant PRIN MIUR 2017-MB8AEZ. The work of M. H-V. has been supported by the Spanish State Research Agency MCIN/AEI/10.13039/501100011033 and by the EU NextGenerationEU/PRTR funds, under grant IJC2020-045126-I. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya.