\section{Introduction}
Quantum droplets are a prime example that the mean-field (MF) description may fail even for weakly interacting Bose gas ~\cite{petrov_quantum_2015,lima_beyond_2012,lima_quantum_2011,bottcher_new_2021}. The MF approach predicts an unstable weakly interacting system, where the corresponding attractive and repulsive contributions nearly cancel each other. This gives way to the enhanced role of zero-point energy fluctuations, usually named the Lee-Huang-Yang (LHY) term~\cite{lee_eigenvalues_1957}, stabilising the emerging droplet. Quantum liquids, however, are more robust in lower dimensions due to the increased quantum fluctuations~\cite{Petrov2016,Santos2017,Malomed2018,Zin2018,Ilg2018,Rakshit2019,Edmonds2020,morera_quantum_2020,morera_universal_2021,guijarro_ultradilute_2022,Ferlaino2022}. 

Particularly, in 1D, quantum droplets have been recently theoretically studied for strongly correlated systems with non-local interactions both in continuum and in the presence of an optical lattice~\cite{oldziejewski_strongly_2020,de_palo_polarization_2021,de_palo_formation_2022,morera_superexchange_2022}. In the former case, a quantum liquid obeys the so-called Lieb-Liniger Gross-Pitaevski~\cite{kopycinski_beyond_2022} equation (LLGPE) that, in the limit of weakly interacting gas, reduces to a standard Gross-Pitaevski equation with the 1D LHY term for short-range interaction. Importantly, in all interaction regimes in 1D, quantum droplets exhibit qualitatively the same ground-state features such as being ultradilute self-bound objects (negative energy for an untrapped state) and having constant chemical potential marked by the flat-top density profile~\cite{bottcher_new_2021}. 


A hallmark of one-dimensional quantum systems with attractive forces is a transition between droplets and solitons \cite{edler_quantum_2017,tylutki_collective_2020,oldziejewski_strongly_2020}. Both are bound states, present for an arbitrary number of particles for different interaction regimes. Namely, a droplet solution appears when short-range repulsion prevails over long-range attraction, contrary to the bright soliton case. Although quantum droplets in 1D exhibit a characteristic flat-top density profile for greater numbers of particles, for smaller systems, one can only distinguish between a bright soliton and a droplet by studying their excitations spectrum. 
This is motivated by the fact, that in a one-dimensional Bose system at least in case without LHY corrections, 
there are no collective modes in a bright soliton -- there are continuum modes solely~\cite{Castin2001}. In contrast, 
the excitation spectrum of a droplet supports small-amplitude collective excitations~\cite{tylutki_collective_2020}.
Moreover, the excitation spectrum differs between 1D and 3D droplets, for the latter supports both bulk and surface modes \cite{Bulgac2002,petrov_quantum_2015}. 

Little is known, however, how the above picture changes by introducing strong correlations between the particles. For instance, for repulsive contact interaction, one can experimentally show that while tuning from the weakly to the strongly interacting regime via a magnetic Feshbach resonance, the dispersion changes dramatically by emerging holelike excitations absent for higher dimension and weaker interactions~\cite{meinert2015}. The recent progress in preparing a bosonic one-dimensional quantum gas of dysprosium~\cite{Kao296} warrants further studies on the exact nature of excitations in one-dimensional systems with simultaneously present attraction and repulsion.


In this work, we investigate the excitation spectrum of a one-dimensional strongly correlated dipolar Bose gas. The complementary problem of solitonic excitations in such system was studied in~\cite{Kopycinski2022b}. The system is described with an effective, hydrodynamic approach introduced in earlier works \cite{oldziejewski_strongly_2020,kopycinski_beyond_2022}. The problem poses a serious challenge in solving because of the presence of non-local dipolar interaction and a quite intricate form of short-range contribution. We thus first focus on a analytical regime of fermionised bosons and vanishing range of dipolar forces. We find a qualitative agreement with previous results about weakly interacting Bose-Bose mixtures~\cite{tylutki_collective_2020} further confirming robustness of quantum droplets in lower dimensions. Next, we generalise our observations to the full problem.

Let us outline the most important findings of our paper. We analyze the excitation spectrum of dipolar quantum droplets by solving the Bogoliubov-de Gennes (BdG) equations. This is done with an emphasis on regimes, where this complex problem simplifies and analytical treatment is at hand. In particular, in the regime of strong contact interactions and short range of dipolar forces, the density profile of a droplet can be found analytically. We thoroughly study this regime obtaining the expression for the droplet energy and characterize the excitation spectrum. What is more, insights from that analysis motivate a very simple rectangular ansatz for the density profiles of droplets. Within approximation of such ansatz, we derive simple and useful formulas for the droplet width and energy. Moreover, we obtain expressions for the excitation spectrum and the number of phononic bound modes that a droplet can exhibit. The ansatz works also outside the aforementioned regime of strong contact interactions and short-ranged dipolar interactions as long as the droplet's shape is dominated by the bulk, flat-top part. In addition to the study of excitations within the framework of BdG equations, we numerically analyze the response to the initial perturbation in the full, nonlinear dynamics.

The paper is organized as follows. In Sec.~\ref{sec:model} we introduce the most important features of the system under study. We briefly describe the effective dipolar interactions in one dimension and the hydrodynamic description of strongly interacting Bose gas. This is followed by the analysis of the ground state of the system. In particular, we specify the regime, for which the density profile of the droplet can be found exactly and introduce the rectangular ansatz for the density profiles. After that, in Sec.~\ref{sec:exc} we present the results for the excitation spectrum obtained from the solution of BdG equations. In Sec.~\ref{sec:nonlin} we broaden this analysis by studying the response to initial perturbation in dynamics given by the nonlinear, hydrodynamic equation. Finally, in Sec.~\ref{sec:summ} we give a summary of results presented in the paper.

\begin{figure}
    \centering
    \includegraphics[width=0.7\textwidth]{graphical-abstract.pdf}
    \caption{We study a gas of atoms trapped in a quasi-1D tube, with dipole moments polarized. Here, we present head-to-tail configuration that leads to attractive interactions between dipoles. Exciting a quantum droplet state may result in an emission of the matter waves. }
    \label{fig:graphical-abstract}
\end{figure}

\section{Model}\label{sec:model}
We consider a system of $N$ bosonic atoms with mass $m$ confined to a quasi-one-dimensional geometry, i.e. we assume a very tight harmonic trap with frequency $\omega_\perp$ in the $y$ and $z$ directions, but particles are unconfined in the 
$x$ direction.
The atoms interact via repulsive contact potential and attractive dipolar forces. As shown in Fig.~\ref{fig:graphical-abstract}, the atomic dipole moments, $\bm{\mu}_D$ are
aligned with an external magnetic field so that the interaction is attractive.
We assume that, the trap in the perpendicular $y$ and $z$ directions is so steep, that the motion in these directions is frozen, i.e. atoms occupy the ground state Gaussians. In this situation, one can integrate out the spatial $y$ and $z$ degrees of freedom to obtain an effective model describing the system restricted to $x$ direction only.

This procedure leads to the effective 1D inter-particle interaction potential
\begin{equation}
    V(x-x')=g \delta(x-x')-g_\text{dd} V_\text{dd}^\sigma(x-x'), \qquad g,g_\text{dd}>0.
\end{equation}
The first term expresses the short-range interaction, with the interaction coefficient  $g$. The second term, coming from the long-range dipolar potential takes the form  \cite{deuretzbacher_ground-state_2010}
 \begin{equation}
   V_\text{dd}^\sigma(x-x')=\frac{v(\frac{x-x'}{\sigma})}{\sigma}, \quad v(u)=\frac{1}{4}\bigg( -2|u|+\sqrt{2 \pi} (1+u^2)e^{u^2/2} \text{Erfc}(|u|/\sqrt{2})\bigg),
 \end{equation}
 where 
$g_\text{dd}=\frac{\mu_0 \mu_D^2(3\cos^2 \theta-1)}{4\pi\sigma^2}$, with $\theta$ denoting the angle between orientation of the dipoles and $x$ axis.
The full dipolar potential reduced to 1D, i.e. $V_\text{dd}^\sigma(x)$ is finite at the origin and it has the characteristic effective range that depends on the oscillatory length of the  transverse confinement $\sigma=\sqrt{\frac{\hbar}{m \omega_\perp}}$. 
In particular, when transverse trap is getting tighter with $\omega_\perp\to\infty$, then $\sigma \to 0$, and the quasi-1D dipolar potential tends to the delta function $V_\text{dd}^\sigma(x) \to \delta(x)$.

In this paper, we wish to study the regime of strong contact interactions. That is why we cannot simply adapt the standard mean-field approach typically used for the description of 1D dipolar quantum droplets, i.e. GP equation supplemented with LHY term \cite{bottcher_new_2021}. Instead, we turn to the effective, hydrodynamic equation described earlier in \cite{kopycinski_beyond_2022,oldziejewski_strongly_2020}, which is based on the density field $\rho(x)$ and the velocity field $v(x)$. Briefly, we use local density approximation and locally, the energy of contact interactions is expressed with the exact ground state energy function $E_\text{LL}(\gamma)$ from the Lieb-Liniger model \cite{lieb_exact_1963,lieb_exact_1963-1}
\begin{equation}
    E_\text{LL}(\gamma)=\frac{\hbar^2}{2m}N \rho^2 e_\text{LL}(\gamma)
    \label{eq:energy-Lieb-Liniger}
\end{equation}
where $\gamma=\frac{m}{\hbar^2}\frac{g}{\rho}$ with $\rho$ denoting density of the gas. In order to account for the density slowly varying in space, we promote the parameter $\gamma$ to a function of position $\gamma \to \gamma(x)=\frac{m}{\hbar^2}\frac{g}{\rho(x)}$. The function $e_\text{LL}(\gamma)$ present in \eqref{eq:energy-Lieb-Liniger} has no known compact analytical form but there exist accurate and handy approximations to it \cite{lang_ground-state_2017,ristivojevic_conjectures_2019,prolhac_ground_2017}. Energy associated with the dipolar interactions is approximated in the standard mean-field way. Precisely, we assume that the dipolar interaction couples densities in different positions and neglect contributions from higher order correlations between atoms. Such an approximation is expected to work when the range of the dipolar interactions $\sigma$ is much larger than the average interparticle spacing defining the fermionization length scale~\cite{oldziejewski_strongly_2020}. The full energy functional reads
\begin{equation}
\begin{aligned}
    E[\phi]&= N\int dx \Big[\frac{\hbar^2}{2m}\left|\frac{\partial \phi}{\partial x}\right|^2+\frac{N^2\hbar^2}{2m}|\phi|^6e_{\rm LL} \bigg(\frac{gm}{\hbar^2 N|\phi|^2}\bigg) + \\
    &-\frac{1}{2}N g_\text{dd} \int dx'V_\text{DD}^\sigma(x-x')|\phi(x)|^2|\phi(x')|^2  \Big].
	\label{eq:energy-LL-GPE-full}
\end{aligned}
\end{equation}
The normalized complex field $\phi(x)$ is directly related to two hydrodynamic fields $\rho(x)$ and $v(x)$ via relation
\begin{equation}
    \phi (x) =\sqrt{\rho(x)/N} e^{i \varphi(x)},
\end{equation}
where $\hbar \partial_x \varphi(x)=m\, v(x)$.
Applying variational principle to the energy functional \eqref{eq:energy-LL-GPE-full} leads to
 \begin{align}
 \label{eq:LLGPdyn}
       i \frac{ \partial\phi(x,t)}{\partial t}&=-\frac{1}{2} \frac{\partial^2 \phi(x,t)}{\partial x^2}+\mu_{\rm LL}\left[N|\phi(x,t)|^2\right]\,\phi(x,t)\nonumber\\
       &-g_\text{dd}N \int dx' V_\text{dd}^\sigma(x-x')|\phi(x', t)|^2\phi(x, t),    
 \end{align}
where, for brevity, we use $\hbar=m=1$. In the equation above we have introduced a functional defining local chemical potential in the Lieb-Liniger model  \cite{lieb_exact_1963}
\begin{equation}
         \mu_{\rm LL}\left[N|\phi|^2\right]=\frac{1}{2}N^2|\phi|^4\left(3e_{\rm LL}\left(\frac{g}{N|\phi|^2}\right)-\frac{g}{N|\phi|^2}e'_{\rm LL}\left(\frac{g}{N|\phi|^2}\right)\right).
  \end{equation}
The stationary version of Eq.~\eqref{eq:LLGPdyn} reads
 \begin{equation}
 \label{eq:LLGPstat}
       \mu\, \phi_0(x)=-\frac{1}{2} \frac{\partial^2 \phi_0(x)}{\partial x^2}+\mu_{\rm LL}\left[N|\phi_0 (x)|^2\right]\,\phi_0(x)-g_\text{dd}N \int dx' V_\text{dd}^\sigma(x-x')|\phi_0(x')|^2\phi_0(x),
\end{equation}
where $\mu$ in turn is the chemical potential of the stationary solution $\phi_0$.
 Let us note that a simplified version of Eq.~\eqref{eq:LLGPstat}, but with an approximate expression for $e_\text{LL}(\gamma)$, was studied before in \cite{oldziejewski_strongly_2020}. Importantly,
 for strong short range interaction and sufficiently weak dipolar potential, the ground state of Eq.~\eqref{eq:energy-LL-GPE-full} is a quantum droplet, with characteristic  flat-top density profile~\cite{oldziejewski_strongly_2020}.
 
In general, Eq.~\eqref{eq:LLGPstat} does not have an analytical solution and one must use numerical methods such as imaginary time evolution to solve it. However, there are some special cases, where the equation simplifies and becomes exactly solvable, yielding closed analytical expressions for density profiles of quantum droplets. In this paper, we will mainly focus on these regimes, as the analytical solution allows for a detailed analysis of elementary excitations and comparison with approximate ansatz solutions.
\subsection{Analytical regime: $g \to \infty$ and $\sigma \to 0$}
In the limit of infinitely strong contact repulsion $g \to \infty$ (fermionization regime with $e_\text{LL}(\gamma \to \infty)=\pi^2/3$) and short-range dipolar interactions $\sigma \to 0$ 
the Eq.~\eqref{eq:LLGPdyn} can be significantly simplified and one finds its exact analytical solutions. Therefore, hereafter we refer to this regime as {\it analytical regime}. In this regime  Eq.~\eqref{eq:LLGPdyn} reduces to \footnote{To keep $g_\text{dd}$ finite upon decreasing $\sigma$, one should simultaneusly tune $\theta$ to make $3 \cos^2 \theta-1$ sufficiently small.}
\begin{equation}
 \label{eq:LLGP_analytic}
      \mu \phi=-\frac{1}{2} \frac{\partial^2 \phi}{\partial x^2}+\frac{\pi^2N^2}{2}|\phi|^4\phi-g_\text{dd}N|\phi|^2\phi,
\end{equation}
which was studied also in  \cite{baizakov_solitons_2009}.
In this special limit the energy functional simplifies to
\begin{equation}
\begin{aligned}
    E[\phi]&= N\int dx \Big[\frac{1}{2}\left|\partial_x \phi\right|^2+\frac{N^2}{2}\frac{\pi^2}{3}|\phi|^6 -\frac{1}{2}N g_\text{dd} |\phi(x)|^4  \Big].
	\label{eq:energy-LL-GPE-analytic}
\end{aligned}
\end{equation}
Note, that the resulting equation contains two competing nonlinear terms, that are different than the nonlinearities appearing for weak interactions only \cite{lima_beyond_2012}.

\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{fig1.pdf}
    \caption{Panel (a): Analytical solutions \eqref{eq:LLGP_analytic_sol} in the fermionization regime, with range of dipolar interactions $\sigma =0$ for different particle numbers $N$ and fixed $g_\text{dd}=30$. The droplet becomes wider with increasing $N$. Moreover, the density profiles become dominated by the bulk, flat-top part and approach a rectangular shape. Panel (b): Fermionization regime with fixed $g_\text{dd}=280$. Numerical solutions of Eq.~\eqref{eq:LLGPstat} with finite dipolar range $\sigma=0.25$ (solid lines) together with $\sigma=0$ analytical solution \eqref{eq:LLGP_analytic_sol} to the Eq.~\eqref{eq:LLGP_analytic} (dashed lines). The differences are small and visible only close to the edges of the droplet. Parameter $\kappa \approx 10.6$ fulfills the condition defining the range of applicability of \eqref{eq:LLGP_analytic}. The atom number is so large that analytical solutions \eqref{eq:LLGP_analytic_sol} are very close to the rectangular ansatz \eqref{eq:rect_ansatz}.}
    \label{fig1}
\end{figure}

The ground state of Eq. ~\eqref{eq:energy-LL-GPE-analytic} has been derived in the context of the nonlinear optics \cite{Pushkarov1979Nov}, and it reads
\begin{equation}
\label{eq:LLGP_analytic_sol}
    \phi^{\infty} (x)=\sqrt{\frac{\sqrt{3}g_\text{dd}}{2\pi \eta}} \frac{\tanh (\eta)}{\sqrt{1+\text{sech} (\eta) \,\cosh (x/a)}},
\end{equation}
where 
\begin{equation}
    \eta =\frac{\pi N}{\sqrt{3}} \qquad a = \frac{\pi}{\sqrt{3}g_\text{dd}\tanh (\eta)}.
\end{equation}
Plots of the solution \eqref{eq:LLGP_analytic_sol} for some exemplary parameters are presented in Fig.~\ref{fig1}(a).
We also derived the total energy \eqref{eq:energy-LL-GPE-analytic} of the ground state solution \eqref{eq:LLGP_analytic_sol}, in terms of its three contributions,  corresponding to the kinetic energy of the envelope $E_{\rm kin}$, the Tonks-Girardeau energy $E_{\rm TG}$ and the dipolar $E_{\rm dip}$ energy. The explicit formulas are 
\begin{equation}
\label{eq:ene1_analytic}
    E^{\infty} _{\rm kin}=\frac{3\sqrt{3}}{16 \pi^3}g_\text{dd}^2 \tanh \eta \Bigg( -\frac{\eta}{\cosh \eta \sinh \eta} +1 \Bigg)
\end{equation}
\begin{equation}
    E^{\infty} _{\rm TG}=\frac{3\sqrt{3}}{16\pi^3} g_\text{dd}^2 \tanh \eta \Bigg( (2+\text{sech}^2 \eta) \eta \coth \eta -3 \Bigg)
\end{equation}
\begin{equation}
    E^{\infty} _{\rm dip}=-\frac{12\sqrt{3}}{16\pi^3} \pi^3\tanh \eta \Bigg( \eta \coth \eta  -1 \Bigg).
\end{equation}
In the limit of large number of atoms one obtains  a simple expression for the total energy (see Fig.~\ref{fig2} for comparison).
\begin{equation}
    E^{\infty} = E^{\infty} _{\rm kin} + E^{\infty} _{\rm TG} + E^{\infty} _{\rm dip} \approx -\frac{3g_\text{dd}^2}{8 \pi^2}N +\frac{3\sqrt{3}g_\text{dd}^2}{8 \pi^3}.
\label{eq:ene_analytic}
\end{equation}
Note that the $N$-independent term in one dimension can be interpreted as surface tension energy. It contains contribution from all types of energies. Its precise value may be useful for study of the phenomena related to fragmentation of droplets~\cite{de_palo_formation_2022}, as it gives the amount of energy required to split a droplet into two smaller ones.

Before further analysis, let us comment on assumptions leading to Eq.~\eqref{eq:LLGP_analytic}. While it is clear that infinitely strong contact interactions are physically reasonable~\cite{kinoshita_observation_2004,paredes_tonksgirardeau_2004}, the zero-range limit of dipolar interactions $\sigma \to 0$ may appear suspicious. Indeed, one of the assumptions leading to the approximate, hydrodynamic equation \eqref{eq:LLGPstat} corresponds to a range of the dipolar interactions much larger than the average interparticle spacing $d$.
However, interparticle spacing does not enter Eq.~\eqref{eq:LLGPstat}  explicitly. 
Therefore, one can find a regime with finite $\sigma$, where the characteristic width of the droplet $W$ is much larger than the range of the dipolar interactions $\sigma$. At the same time, this regime is characterized by $\sigma$ much larger than the average interparticle spacing approximately given by $W/N$. We observe that in such a regime, given by the condition $W \gg \sigma \gg W/N$, one can perform the limit $\sigma \to 0$ in Eq.~\eqref{eq:LLGPstat} obtaining roughly the same density profile as compared to the one with finite dipolar range $\sigma$. We confirm that on Fig.~\ref{fig1}(b). 
Later in the text, we estimate width of the droplets as $W=2\pi^2N/3g_\text{dd}$, hence we can write a condition for reliability of Eq.~\eqref{eq:energy-LL-GPE-full} in terms of parameters of our model in the following way $N \gg \kappa \gg 1$, where $\kappa:= 3 g_\text{dd} \sigma/2 \pi^2$.




 
 The plots presented in Fig.~\ref{fig1}(b) suggests yet another approximation, which will be discussed in the next section.
 



\begin{figure}[h]
    \centering
    \includegraphics[width=0.8\textwidth]{fig2.pdf}
    \caption{Analytical regime: energy differences between exact energy $E^{\infty}$ and its large $N$ asymptotics (the right hand side of Eq.~\eqref{eq:ene_analytic}) and the energy \eqref{eq:exc_rect} of the rectangular ansatz approximation with fixed $g_\text{dd}=1$. The differences are negligible for large enough $N$. Inset: Comparison between the exact energy $E^{\infty}$ together with large $N$ asymptotics and rectangular ansatz formula \eqref{eq:exc_rect}. Dipolar attraction is fixed at $g_\text{dd}=1$.}
    \label{fig2}
\end{figure}
\subsection{Rectangular ansatz for the ground state density  \label{subsec:rect-ansatz}}
Looking at Fig.~\ref{fig1}(a) one can notice that when $N$ becomes larger and larger, the solution \eqref{eq:LLGP_analytic_sol} starts to resemble a rectangle. This motivates us to introduce the following,  very simple ansatz
\begin{equation}
\label{eq:rect_ansatz}
    \phi_W(x) =\frac{1}{\sqrt{W}} \text{rect}(x/W),
\end{equation}
where the width $W$ is the only parameter of the ansatz. Although we motivate it mainly by the shape of solutions in the analytical regime, we present results for more general case $g \neq \infty$ and $\sigma \neq 0$. In principle, flat-top shapes resembling rectangle may appear outside the range of parameters for which \eqref{eq:LLGPstat} is exactly solvable. We plug our ansatz to the energy functional Eq.~\eqref{eq:LLGPstat} and minimize energy with respect to $W$. Additionally, we neglect term involving $\frac{d\phi}{dx}$, in the same spirit as neglecting the kinetic energy in the Thomas-Fermi approximation. Namely, in the number of parameters we checked that the kinetic energy is much smaller  than interaction energy,  stemming from the other two terms. This can be confirmed in the analytical regime. The asymptotic form of the kinetic energy \eqref{eq:ene1_analytic} becomes $N$-independent there, and thus its contribution to the whole energy \eqref{eq:ene_analytic} is negligible for large $N$. Moreover, to simplify the integral involving dipolar potential, we assume that $W \gg \sigma$. 
Under these assumptions, we look for the width $W$ minimizing energy
\begin{equation}
    \frac{dE}{dW}=\frac{1}{2}N^3\Bigg[-\frac{2}{W^3}e_\text{LL}\Big(\frac{gW}{N}\Big)+\frac{1}{W^3}\frac{gW}{N}e_\text{LL}'\Big(\frac{gW}{N}\Big)\Bigg]+g_\text{dd}\frac{N^2}{2W^2}=0.
\end{equation}
The equation above may be formulated in a  compact form
\begin{equation}
\label{eq:minimum}
    \frac{1}{\gamma_W}\big(2e_\text{LL}(\gamma_W)-\gamma_W e_\text{LL}'(\gamma_W)\big)=g_\text{dd}/g,
\end{equation}
 with optimal $W$ hidden in $\gamma_W:=gW/N$. From properties of the $e_\text{LL}(\gamma)$ function we observe that this equation has a solution only if $g_\text{dd}/g<1$, and therefore we expect flat-top solutions only if contact repulsive interactions are stronger than dipolar attraction $g>g_\text{dd}$. 

 This is in agreement with previous observations~\cite{oldziejewski_strongly_2020,edler_quantum_2017} suggesting soliton-droplet crossover near $g \approx g_\text{dd}$.
Additionally, we may calculate the corresponding chemical potential and energy in the rectangular approximation
\begin{equation}
    \mu=\mu_\text{LL}[N/W]-g_\text{dd}\frac{N}{W}, \qquad E=\frac{N^3}{2W^2}e_\text{LL}(\gamma_W)-\frac{N^2 g_\text{dd}}{2W}.
    \label{eq:muEneBox}
\end{equation}
In the limit $g \to \infty$ equation \eqref{eq:minimum} becomes exactly solvable, leading to the explicit formulas for the droplet width 
\begin{equation}
    W=\frac{2 \pi^2 N}{3g_\text{dd}},
\end{equation}
the chemical potential and the energy
\begin{equation}
\label{eq:muEneBox_analytic}
    \mu_W =-\frac{3g_\text{dd}^2}{8 \pi^2}, \qquad E_W=-\frac{3g_\text{dd}^2N}{8 \pi^2}.
\end{equation}
In the analytical regime we can directly compare rectangular ansatz energy to the exact expressions. From comparison of Eqs.~\eqref{eq:ene_analytic} and \eqref{eq:muEneBox_analytic}, we see that rectangular ansatz gives correct value of the term proportional to $N$ but neglects $N-$ independent constant corresponding to surface tension energy. This is in agreement with our initial assumption of large $N$. Additional results in this approach, for energy differences,  are presented in Fig.~\ref{fig2}.

The rectangular ansatz leads to simple analytic formulas characterizing parameters of a droplet. In the remaining part of the paper we use  it further to derive and understand the elementary excitations in the system.
\section{Elementary excitations}\label{sec:exc}
The main goal of this paper is to study excitations of 1D dipolar droplets. Elementary excitations govern their low-energy dynamics, characterize the response to small perturbations and can be used to study the low-temperature thermodynamics. To access elementary excitations and energies, we linearize our equations of motion \eqref{eq:LLGPdyn}. This procedure, on a formal level, is equivalent to the standard Bogoliubov-de~Gennes (BdG) framework applied to nonlinear Schrodinger-like equations \cite{kaiser_bose-einstein_2001}. In general, we consider a stationary solution $\phi_0(x)$ to the Eq.~\eqref{eq:LLGPdyn} and wish to characterize a response to some small perturbation of that field.  We consider the following, standard ansatz for the perturbed time-dependent field
\begin{equation}
\label{eq:linansatz}
    \phi(x,t)=\Big(\phi_\text{0}(x)+\delta \phi(x,t)\Big)e^{-i\mu t},
\end{equation}
where $\mu$ is the chemical potential of stationary solution $\phi_0(x)$ and the small dynamical perturbation is assumed in the form \cite{pitaevskii_bose-einstein_2016,kaiser_bose-einstein_2001}
\begin{equation}
\label{eq:deltaansatz}
    \delta \phi(x,t)= u(x)e^{-i\epsilon t}+v^*(x)e^{i\epsilon t}
\end{equation}
where the functions $u(x)$ and $v(x)$ characterize the shape of the perturbation and $\epsilon$ denotes the excitation energy. Derivation of BdG equations in our case amounts to substituting the perturbed state, as given in Eq.~\eqref{eq:linansatz}, to the LL-GPE, given in Eq.~\eqref{eq:LLGPdyn} and linearizing it. We get
\begin{subequations}
\begin{align}
    &\epsilon u(x)= \Big( -\frac{1}{2}\partial_x^2 -\mu +f_1(x)+\mathcal{J}_1+\mathcal{J}_2 \Big) u(x)+\Big(f_2(x)+\mathcal{J}_2\Big) v(x),
    \label{eq:bdguv1}
\\
    -&\epsilon v(x)=\Big( -\frac{1}{2}\partial_x^2 -\mu +f_1(x)+\mathcal{J}_1+\mathcal{J}_2 \Big) v(x)+\Big(f_2(x)+\mathcal{J}_2\Big) u(x),
    \label{eq:bdguv2}
\end{align}
\end{subequations}
where we defined two functions related to the LL ground state
\begin{subequations}
\begin{align}
    &f_1(x):=\frac{1}{2}N^2\phi_0^4(x)\Big(9e_\text{LL}\big(\gamma(x)\big)-5\gamma(x)e_\text{LL}'\big(\gamma(x)\big)+\gamma^2(x)e_\text{LL}''\big(\gamma(x)\big)\Big),
\\
&f_2(x):=\frac{1}{2}N^2\phi_0^4(x)\Big(6e_\text{LL}\big(\gamma(x)\big)-4\gamma(x)e_\text{LL}'\big(\gamma(x)\big)+\gamma^2(x)e_\text{LL}''\big(\gamma(x)\big)\Big),
\end{align}
\end{subequations}
and two operators related to non-local interaction:
\begin{subequations}
\begin{align}
    &\mathcal{J}_1h(x) :=-g_\text{dd}N\int dx'\,V_\text{dd}^\sigma(x-x')\phi_0^2(x') h(x),
\\
&\mathcal{J}_2h(x):=-g_\text{dd}N\int dx'\,V_\text{dd}^\sigma(x-x')\phi_0(x')\phi_0(x) h(x').
\end{align}
\end{subequations}
Following the lines of Ref.~\cite{Ronen2006}, we introduce new variables $r(x):=u(x)+v(x)$ and $s(x):=v(x)-u(x)$ transforming our equations to
\begin{equation}
    \begin{pmatrix}
    0 & \hat{\mathcal{A}} \\
    \hat{\mathcal{B}} & 0
    \end{pmatrix}\begin{pmatrix}
  r\\
  s
  \end{pmatrix}=\epsilon \begin{pmatrix}
  r\\
  s
  \end{pmatrix}.
  \label{eq:rsrelation}
\end{equation}
with 
\begin{subequations}
\begin{align}
    &\hat{\mathcal{A}}:=-\frac{1}{2}\partial_x^2 - \mu + f_1(x)-f_2(x)+\mathcal{J}_1,
\\
    &\hat{\mathcal{B}}:=-\frac{1}{2}\partial_x^2 - \mu + f_1(x)+f_2(x)+\mathcal{J}_1+2\mathcal{J}_2.
\end{align}
\end{subequations}
We may further simplify our problem by acting with the matrix above once again on the Eq. \eqref{eq:rsrelation}. We get
\begin{subequations}
\begin{align}
    &\hat{\mathcal{A}}\hat{\mathcal{B}}r(x)=\epsilon^2r(x),
    \label{eq:bdgfinal1}
\\
    &\hat{\mathcal{B}} \hat{\mathcal{A}}s(x)=\epsilon^2s(x).
    \label{eq:bdgfinal2}
\end{align}
\end{subequations}
Note that it is sufficient to solve the eigenproblem for just one equation \eqref{eq:bdgfinal1} or \eqref{eq:bdgfinal2}. The remaining eigenvectors (either $r(x)$ or $s(x)$) can be found from the Eq. \eqref{eq:rsrelation} relating $r(x)$ and $s(x)$ functions. The normalization of $u(x)$ and $v(x)$ is given by the standard condition
\begin{equation}
\label{eq:uvnorm}
    \int dx \Big(|u(x)|^2-|v(x)|^2\Big)=\frac{1}{N}
\end{equation}
that for modes with non-zero energy may be expressed with a condition for $r(x)$ function as $\epsilon^{-1}\int dx  \,r(x) \hat{\mathcal{B}} r(x) =\frac{1}{N}$.


Note that using Eq.~\eqref{eq:linansatz}, we may write down the approximate time evolution of a density profile
\begin{equation}
\label{eq:densityPert}
    \rho(x,t)=|\phi(x,t)|^2\approx|\phi_0(x)|^2+2\phi_0(x)\Big(u(x)+v(x)\Big)\cos(\epsilon t).
\end{equation}
Therefore, function $r(x)=u(x)+v(x)$ is directly related to the shape of the density perturbation.

The BdG equations \eqref{eq:bdgfinal1} and \eqref{eq:bdgfinal2} are solved in the momentum space. We discretize the space and work on finite matrices. In the analytical regime, we find that the low-lying energies converge when one takes sufficiently large number of the numerical grid points. However, outside of the analytical  regime, first one has to find numerically the stationary solution $\phi_0$ of Eq.~\eqref{eq:LLGPstat} to later use it in our solver for the BdG equations. 
Numerical solution of \eqref{eq:LLGPstat} has some inevitable errors, which later strongly amplify in the numerics related to the BdG equations, strongly affect reliability of the excitation energies. We did not fully overcome these problems, therefore we focus mainly on the excitations in the analytical regime and try to draw some conclusions for the more general case. The general case is studied then numerically by studying the evolution of a perturbed state in the frame of the time-dependent LLGPE. 
\subsection{Analytical regime: $g \to \infty$ and $\sigma \to 0$}
In the analytical regime,  $g \to \infty$ and $\sigma \to 0$, we use the exact solution for the density profile $\phi_0$ to find numerically its excitations.

We illustrate our result in Fig.~\ref{fig:fig3} for a specific parameters corresponding to the analytical regime.
In Fig.~\ref{fig:fig3}(a) we see several lowest-lying excitation energies for an exemplary quantum droplet, whose density profile is presented in the inset of that figure. We have two zero-energy modes related to the breaking of the translational and phase symmetries. Next, we have a couple of modes that we call "bound" as their energy is below chemical potential and their profiles vanish outside of the droplet. These modes are either symmetric or antisymmetric functions of the position (see Figs.~\ref{fig:fig3}(b,c)). Importantly, apart from a small region near the edges of the droplet, the modes strongly resemble consecutive standing waves of the infinite box potential with the width given by $W$, implying that the edges of droplets effectively impose open boundary conditions. Above the absolute value of chemical potential we have scattering modes, with non-zero probability of finding a particle outside the quantum droplet (cf. Fig.~\ref{fig:fig3}(d)). The discrete character of the spectrum of these modes is inherited from the imposed periodic boundary conditions. In general, the scattering modes have continuous spectrum. We checked that by accessing different momenta via varying the size of the periodic box.
\begin{figure}[h]
    \centering
    \includegraphics[width=0.75\textwidth]{fig3.png}
    \caption{Panel (a): excitation energies found from the solution of BdG equations for $N=30$ and $g_\text{dd}=100$. There are two zero-energy modes and eight bound modes with energies below chemical potential. The rest can be classified as scattering modes. In the inset we present density profile of the stationary solution. Additionally, we plot excitation energies found in the rectangular ansatz that match well with the energies of bound modes. Panels (b) and (c): shapes of exemplary bound modes (antisymmetric and symmetric, respectively). Panel (d): an example of a scattering mode. Note the non-zero probability of finding particle outside the droplet. }
    \label{fig:fig3}
\end{figure}

 This structure of elementary modes turns out to be generic in the considered regime of $g\to \infty$ and $\sigma \to 0$. In the Fig.~\ref{fig:fig4} we present how excitation spectrum of the bound modes changes when system parameters $N$ and $g_\text{dd}$ are varied. As we see in Fig.~\ref{fig:fig4}(a), the number of modes grows with the particle number $N$. Moreover, from Fig.~\ref{fig:fig4}(b) we see that for fixed $N$, number of modes is independent of $g_\text{dd}$ and the excitation energies grow with $g_\text{dd}$. Interestingly, these observations can be elucidated with an approximate simple analytical approach, based on the rectangular ansatz introduced earlier.
 \subsubsection{Rectangular ansatz approach to the excitations}
 In the Sec. \ref{subsec:rect-ansatz} we show, that width, energy and density of a droplet may be approximated by the results based on a simple rectangular ansatz \eqref{eq:rect_ansatz}..

 We utilize this approach further and solve BdG equations taking $\phi_0(x)=\phi_W(x)$. In such case BdG equations may be written as
 \begin{subequations}
\begin{align}
    &\epsilon u(x)= \Big( -\frac{1}{2}\partial_x^2 +\pi^2 (N/W)^2-g_\text{dd} N/W\Big) u(x)+\Big(\pi^2(N/W)^2-g_\text{dd}N/W\Big) v(x),
\\
    -&\epsilon v(x)=\Big( -\frac{1}{2}\partial_x^2 +\pi^2 (N/W)^2-g_\text{dd} N/W \Big) v(x)+\Big(\pi^2(N/W)^2-g_\text{dd}N/W\Big) u(x).
\end{align}
\end{subequations}
We are looking for solutions in the form of standing waves $u(x)=u \cos (px)$ and $v(x)=v\cos (px)$ with $p=\frac{\pi}{W}, \frac{2\pi}{W},\ldots$. Such values of allowed momenta lead to vanishing of the modes at the boundary of the droplet, as expected. These equations can be readily solved giving
\begin{equation}
\label{eq:exc_rect}
    \epsilon(p)=\sqrt{\Big((\pi N/W)^2-g_\text{dd}N/W\Big)p^2+\Big(\frac{p^2}{2}\Big)^2}.
\end{equation}
Using $W=\frac{2 \pi^2 N}{3g_\text{dd}}$ we obtain explicit excitation energy of $m$-th mode with momentum $p=\frac{m \pi}{W}$
 \begin{equation}
    \epsilon_m=\frac{m}{N} \Big( \frac{3}{2 \pi} \Big)^2\sqrt{\frac{1}{3}+\frac{m^2}{4N^2} }\, g_\text{dd}^2.
\label{eq:excbox}
\end{equation}
Hence, excitation energies for fixed $N$ grow quadratically with $g_\text{dd}$ (cf. Fig.~\ref{fig:fig5}(b)). For large $N$ and fixed $g_\text{dd}$, energies decay 
as $\propto N^{-1}$ (cf. Fig.~\ref{fig:fig5}(a)) and are proportional to $m$, resembling phononic dispersion relation.

In principle, the spectrum given in Eq.~\eqref{eq:excbox} consists of infinitely many eigenenergies. We phenomenologically cut it on the level $\epsilon_{m_\text{max}}=-\mu_W = \frac{3g_\text{dd}^2}{8 \pi^2}$. This gives  condition for the maximal number of modes inside a droplet $m_\text{max}$:
\begin{equation}
    m_\text{max}=\sqrt{\frac{1}{3}(\sqrt{5}-2)}N\approx0.28 N.
\end{equation}
The number of  bound modes grows linearly with the particle number $N$ and it is independent on $g_\text{dd}$ in accordance with Fig.~\ref{fig:fig4}(b).

Results from rectangular ansatz approach are compared to the full solution of the BdG equations in Figs.\ref{fig:fig4}(a,b) and Fig.~\ref{fig:fig3}(a). We see a good agreement, especially for the lowest-lying excitations with the eigenenergies far from $-\mu$. The framework built on the rectangular ansatz allows for particularly simple analytical description of droplet excitations that capture the most important features. In particular, excitation modes may be approximately seen as a cosine standing waves of a box with a width given by $W$ and the associated energy of $m$-th mode is determined from \eqref{eq:excbox}.
 \begin{figure}[h!]
     \centering
     \includegraphics[width=0.75\textwidth]{fig4.png}
     \caption{Panel (a): Excitation spectrum of bound modes for a droplet with $g_\text{dd}=5$ as a function of particle number $N$. We see that as $N$ grows, the number of bound modes increases and their energies get smaller. Panel (b): Excitation spectrum of bound modes for $N=30$ particles as a function of $g_\text{dd}$. The energies are growing functions of $g_\text{dd}$. The number of the bound modes does not depend on $g_\text{dd}$. In both panels dashed lines mark the results obtained within the framework of rectangular ansatz, Eq.~\eqref{eq:excbox}.} 
     \label{fig:fig4}
 \end{figure}

At this point, it is worth to address two issues. First, analogous excitations were already studied in the droplets in a weakly interacting Bose-Bose mixtures~\cite{tylutki_collective_2020}, that results from the Gross-Pitaevskii equations extended by the LHY terms~\cite{petrov_ultradilute_2016}. That results are qualitatively similar to our findings. In both setups, droplets have bound modes displaying phononic dispersion relation and scattering modes. Moreover, the number of bound modes is finite and scales proportionally to $N$ for large droplets. To understand these analogies, one can note that both equations are similar in the sense that they are local equations with competing nonlinearities. In fact, the only difference lies in the different power-law dependence of the terms corresponding to the repulsion and attraction in the system. This analogy in principle may be broken in the regime, where the range of interactions is comparable to the width of the droplet and Eq.~\eqref{eq:LLGPstat} becomes truly non-local. 

Solving the BdG equations for general case of finite $g$ and $\sigma$ is difficult. On the other hand, we have observed a good agreement between numerical solutions of BdG equations and analytical results based on the simple rectangular ansatz. The latter can be applied for any $g$ and $\sigma$ and practically require only flat-top density profile. Therefore, in the next section we briefly generalize results from this section to parameters outside analytical regime.
\subsection{General case in the flat-top regime}
The BdG equations \eqref{eq:bdguv1} and \eqref{eq:bdguv2} used for the case of the rectangular ansatz, but without imposing any restrictions on $g$ and $\sigma$ take the form
\begin{subequations}
\begin{align}
    &\epsilon u(x)= \Big( -\frac{1}{2}\partial_x^2 +v_\text{LL}^2[N/L]+\mathcal{J} \Big) u(x)+\Big(v_\text{LL}^2[N/L]+\mathcal{J}\Big) v(x),
    \label{eq:unibdg1}
\\
    -&\epsilon v(x)= \Big( -\frac{1}{2}\partial_x^2 +v_\text{LL}^2[N/L]+\mathcal{J} \Big) v(x)+\Big(v_\text{LL}^2[N/L]+\mathcal{J}\Big) u(x),
    \label{eq:unibdg2}
\end{align}
\end{subequations}
where we defined the following operator
\begin{equation}
    \mathcal{J} h(x) = -\frac{g_\text{dd} N}{W} \int dx' V_\text{dd}^\sigma(x-x') h(x'),
\end{equation}
and the Lieb-Liniger speed of sound \cite{lieb_exact_1963,ristivojevic_excitation_2014}
\begin{equation}
v_\text{LL}[N/W]=\frac{N}{W}\sqrt{ 3e_{\rm LL}(\gamma_W)-2\gamma_W e_{\rm LL}'(\gamma_W)+\frac{1}{2}\gamma_W^2 e_{\rm LL}''(\gamma_W)}.
\end{equation}
Let us remind that the width of the optimal ansatz rectangle $W$ is determined here numerically from the Eq.~\eqref{eq:minimum}. We find solutions of these equations in the form of standing waves and obtain their spectrum:
\begin{equation}
 \epsilon(p)=\sqrt{\Big(v_\text{LL}^2[N/W]-\frac{g_\text{dd} N}{W}\hat{V}_\text{dd}^\sigma(p)\Big)p^2+\Big(\frac{p^2}{2}\Big)^2}.
 \label{eq:exc_rect_gen}
\end{equation}
where $p=\frac{\pi}{W}, \frac{2\pi}{W},\ldots$ and 
$\hat{V}_\text{dd}^\sigma(p)$ denotes the Fourier transform of the dipolar potential. The number of bound modes can be determined from the transcendental equation $\epsilon(\pi m_\text{max} /W)=-\mu=-\mu_\text{LL}[N/W]+g_\text{dd}N/W$. By taking the limits $g \to \infty$ and $\sigma \to 0$ on the Eq.~\eqref{eq:exc_rect_gen} we obtain, as a special case,  the formula \eqref{eq:exc_rect}. 

Let us emphasize once again that the results for nonzero $\sigma$ and finite $g$ rely on the assumption, that the stationary density profile is well approximated by a rectangle. In general, this does not have to be the case, as the Eq.~\eqref{eq:LLGPstat} has solutions that are closer to Gaussian-like solitonic shapes \cite{oldziejewski_strongly_2020}. In such situations, excitation spectrum will differ significantly and one should attempt to solve the full BdG equations.



\section{Nonlinear response to a perturbation}\label{sec:nonlin}
In this section we wish to go beyond the linearized version of the dynamics given by Eq.~\eqref{eq:LLGPdyn}. Precisely, to understand better the low energy dynamics of the droplets we numerically solve the full time-dependent equation \eqref{eq:LLGPdyn}. As initial states we take  droplets perturbed with the modes found from BdG equations, analyzed in detail in the previous section. In particular, we want to check the differences between linearized dynamics and the full, nonlinear one given by \eqref{eq:LLGPdyn}.

The results are presented in Fig. 
\ref{fig:fig5}. In the first two cases, corresponding to Fig. \ref{fig:fig5} (a,b), we consider a droplet with the same parameters as in Fig. \ref{fig:fig3}, namely $N=30$, $g_\text{dd}=100$ and $\sigma=0$. In Fig.~\ref{fig:fig5}(a) we study the dynamics of a bound mode with $m=3$ with the excitation energy given by \eqref{eq:excbox} and the shape proportional to cosine standing wave $u(x),v(x) \sim \cos \frac{3 \pi x}{W}$.

\begin{figure}[h!]
    \centering
    \includegraphics[width=0.75\textwidth]{fig5.png}
    \caption{Dynamics of the perturbed droplet with $N=30$, $g_\text{dd}=100$ and $
   \sigma=0$ (first row) or $
    \sigma=0.05$ (second row). First column: time evolution of bound cosine mode with $m=3$, revival time of the wave matches the period $T_\text{exc}$ related to the excitation energy  given by \eqref{eq:excbox} for the first row and \eqref{eq:exc_rect_gen} for the second row. Second column: dynamics of $m=14$ cosine mode, with the excitation energy above the chemical potential. We see the emission of particles in the course of time evolution.}
    \label{fig:fig5}
\end{figure}
As the excitation energy is below chemical potential, such a perturbation does not lead to the emission of particles from the droplet. After some time, the nonlinear character of evolution brings additional frequencies to the dynamics. Nevertheless, for short times the evolution is governed by a single frequency, given by the excitation energy. This can be seen from the good agreement between revival time of the excited wave and $T_\text{exc}=2 \pi /\epsilon_3$, cf. black line on Fig. \ref{fig:fig5}(a). The situation is drastically different when we consider a mode with an excitation energy higher than the chemical potential. Droplet for our parameters supports eight bound modes, and perturbing it with higher cosine modes should lead to emission of particles. As we see in Fig.~\ref{fig:fig5}(b), where the dynamics of $m=14$ cosine mode is presented, this is indeed the case. After short time, particles are emitted and the dynamics in the bulk cannot be described as a single wave that propagates inside and is reflected from the edges.


In Figs.\ref{fig:fig5}(c,d) we consider evolution of the same modes, but in the system with nonzero range of dipolar interactions $\sigma=0.05$. We observe that this change makes no qualitative difference in the behavior of the system, suggesting that the general picture of excitations seen in the analytical regime may be valid for flat-top droplets with longer range of dipolar forces.
The dynamics presented in this section confirms and supplements the conclusions drawn from the analysis of excitations in the BdG framework. 

\section{Summary and outlook}\label{sec:summ}
We have studied flat-top dipolar Tonks-Girardeau droplets, namely quasi 1D droplets in the regime of strong contact interaction. We derived simple approximate formulas for the most important quantities characterizing such a droplet. In  between we have found the formulas for the energy of the bulk, and the energy of a droplet edge, i.e. its "surface". We presented the excitations spectrum and the shapes of the excitations. Our analysis was based mainly on two approaches.

In the first approach we focused on the {\it analytical regime}, where the energy becomes a local operator, and its ground state is analytically known \cite{baizakov_solitons_2009}. We linearize the equation of motion around the ground state to derive the local BdG equations and find numerically the exact excitation spectrum.  Our results can be applied in the limit of large $g$ and when the dipolar interaction coefficient $g_{\rm dd}$, the number of atoms  $N$ and the transverse oscillator length $\sigma$ fulfill the conditions  $N \gg 3 g_\text{dd} \sigma/2 \pi^2 \gg 1$.

The second approach was based on a trivial ansatz -- approximating droplet density with a rectangular shape and minimizing total energy with respect to its width. Using it, we derived formulas for the excitation spectrum and the number of "bound modes", i.e. internal droplet excitations, also beyond the validity regime of the first approach.

We benchmarked the two approaches in the common validity regime, obtaining excellent agreement. We complement this analysis by the simulations of the perturbed droplet dynamics, using the full LLGPE \eqref{eq:LLGPdyn}. We confirmed the stability of the droplets to perturbation with a bound mode and show that to perturbations with energy above chemical potential induce particle emission from the droplet. 

To summarize, the elementary excitations can be characterized as follows. There are two zero-energy modes related to broken symmetries in the system. Next, there is a finite number of bound modes with excitation energies below the chemical potential, whose shapes strongly resemble standing waves of a box with width given by the size of the droplet. These excitations have phononic dispersion relation in the limit of large particle number $N$, moreover the number of these modes is finite and scales proportionally to $N$. In addition to that, a droplet has a continuum spectrum of scattering modes characterized by a non-zero probability of finding a particle outside the droplet.


\section*{Acknowledgements}
We thank Jan Chwedeńczuk and Maciej Marciniak for valuable discussions.
J.K., M.Ł. and K.P. acknowledges support from the (Polish) National
Science Center Grant No. 2019/34/E/ST2/00289. 
W.G. was supported by the Foundation for Polish Science (FNP) via the START
scholarship.
Center for Theoretical Physics of the Polish Academy of Sciences is a member of the National Laboratory of Atomic, Molecular and Optical Physics (KL FAMO).

