# Tuning Kerr-Soliton Frequency Combs to Atomic Resonances

Su-Peng Yu<EMAIL_ADDRESS>Time and Frequency Division, NIST, Boulder,
Colorado, USA  Department of Physics, University of Colorado, Boulder,
Colorado, USA  Travis C. Briles  Gregory T. Moille  Xiyuan Lu  Center for
Nanoscale Science and Technology, NIST, Gaithersburg, Maryland, USA  Scott A.
Diddams  Time and Frequency Division, NIST, Boulder, Colorado, USA  Department
of Physics, University of Colorado, Boulder, Colorado, USA  Kartik Srinivasan
Center for Nanoscale Science and Technology, NIST, Gaithersburg, Maryland, USA
Scott B. Papp<EMAIL_ADDRESS>Time and Frequency Division, NIST, Boulder,
Colorado, USA  Department of Physics, University of Colorado, Boulder,
Colorado, USA

###### Abstract

Frequency combs based on nonlinear-optical phenomena in integrated photonics
are a versatile light source that can explore new applications, including
frequency metrology, optical communications, and sensing. We demonstrate
robust frequency-control strategies for near-infrared, octave-bandwidth
soliton frequency combs, created with nanofabricated silicon-nitride ring
resonators. Group-velocity-dispersion engineering allows operation with a 1064
nm pump laser and generation of dual-dispersive-wave frequency combs linking
wavelengths between approximately 767 nm and 1556 nm. To tune the mode
frequencies of the comb, which are spaced by 1 THz, we design a photonic chip
containing 75 ring resonators with systematically varying dimensions and we
use 50 ∘C of thermo-optic tuning. This single-chip frequency comb source
provides access to every wavelength including those critical for near-infrared
atomic spectroscopy of rubidium, potassium, and cesium. To make this possible,
solitons are generated consistently from device-to-device across a single
chip, using rapid pump frequency sweeps that are provided by an optical
modulator.

##  I Introduction

Recently, there has been significant effort to develop chip-scale and
integrated-photonics light sources and devices through nonlinear optics Moss
_et al._ (2013). These have ranged from the development of ultra-stable CW
lasers based on stimulated Brillouin scattering Li  _et al._ (2013); Loh  _et
al._ (2015) to metrology with frequency combs based on supercontinuum
generation in photonic waveguides Carlson  _et al._ (2017); Lamb  _et al._
(2018) to dissipative-Kerr-soliton (DKS) formation in microresonators
Kippenberg  _et al._ (2018). Such soliton microcombs exhibit interesting
nonlinear optical phenomena, including breathing oscillations Bao  _et al._
(2016); Lucas  _et al._ (2017); Yu  _et al._ (2017), dark-pulse formation Xue
_et al._ (2015) and soliton crystallization Cole  _et al._ (2017), and have
enabled diverse applications like communications Marin-Palomo  _et al._
(2017), and optical-frequency synthesis Spencer  _et al._ (2018). DKS combs
have been demonstrated with a variety of materials, including fused silica
fiber and microresonators Yi  _et al._ (2015); Wang  _et al._ (2017),
crystalline magnesium fluoride Matsko  _et al._ (2016) and silicon nitride
(Si3N4 hereafter SiN)Brasch  _et al._ (2016); Jaramillo-Villegas  _et al._
(2015); Briles  _et al._ (2018a). SiN is an appealing material for broad
bandwidth combs due to its moderately high Kerr nonlinearity and low optical
propagation loss. Additionally, it is compatible with standard lithographic
fabrication techniques, allowing for reliable fabrication of resonators and
integration with other photonic elements Moss  _et al._ (2013), providing a
path for broadband, low-power consumption frequency combs on a photonic chip
Stern  _et al._ (2018).

Realizing DKS frequency combs that interface with optical-atomic transitions
opens up timekeeping and other measurement applications. However alkali or
alkaline-earth atomic resonances are largely prevalent in the 700 to 1000 nm
spectral region, whereas many DKS experiments have operated in the 1550 nm
band due to the availability of pump lasers, optical amplifiers, filters,
detectors, and other components. Previous optical-clock demonstrations with
soliton microcombs were pumped at 1550 nm Papp  _et al._ (2014), and required
second-harmonic generation to access atomic transitions. Moreover, atomic
transitions are spectrally narrow and thus require finely tuning a frequency
comb to reach resonance.

A central facet of designing DKS combs for the near-infrared Karpov  _et al._
(2018) is that the very high microwave-to-THz repetition frequencies of
microcombs and their ultrabroad-bandwidth spectra are almost exclusively
determined by the geometry of the devices. A recipe has emerged to design and
realize soliton microcombs at 1550 nm for some applications like optical-
frequency synthesis Briles  _et al._ (2018a). The first step in designing the
soliton spectrum is to select device geometries with anomalous group-velocity-
dispersion (GVD) near the chosen pump wavelength Agha  _et al._ (2009). A
second step is to design higher than 2${}^{\textnormal{nd}}$\- order GVD
contributions for generation of dispersive waves (DWs) Okawachi  _et al._
(2014), which lead to a localized increase in optical power at the edges of
the soliton spectrum. A third step is tuning the carrier-envelope offset
frequency of the soliton microcomb for compatibility with electronic detection
on available highspeed photodetectors with bandwidth $\leq$ 20GHz.

Figure 1: Resonator dispersion design and experimental setup. (a) The
calculated GVD profile (top) and integrated modal dispersion (middle) for ring
resonators with material thickness 667nm, radius of 22.5 $\mu$m and RW of 810
nm, 850 nm, and 890 nm. The experimentally generated single soliton spectrum
from the RW=850 nm device (bottom) has its strongest SDW line aligned to the
Rb 2-h$\nu$ transition (red dashed line). (b) The experimental setup for
soliton generation and comb diagnostics. An ECDL is coupled to a Mach-Zehnder
intensity-modulator, amplified in a YDFA and coupled into the photonic chip
using lensed fibers and precision positioning stages. The generated comb light
is separated into spectral bands centered at 780, 1064 and 1550 nm for
subsequent diagnostics. Each branch can be monitored in real time or undergo
heterodyne with local oscillators. The notch filter used in the 1064 nm branch
to reject pump light when monitoring soliton generation, is not shown. The
long-wavelength DW branch is amplified using EDFA to provide optical power for
further measurement steps.

Here, we explore frequency control and tuning of octave-spanning, near-
infrared (NIR) soliton microcombs pumped at 1064 nm, towards use in atomic
spectroscopy. We demonstrate a combination of GVD control of the comb’s
spectrum and thermo-optic frequency tuning for a single-chip solution to
interface with rubidium (780 nm), cesium (852 nm), and potassium (767 nm)
atoms. The advantage of a 1064 nm pump laser is that it bisects the frequency
span between telecom wavelengths (1300 to 1600 nm) and atomic species that are
commonly used for atomic clocks like strontium (699 nm) and calcium ion (729
nm) Diddams  _et al._ (2004); Ludlow  _et al._ (2015); Targat  _et al._
(2013). Our soliton microcomb provides up to 0.3$\mu$W per mode over the 760
to 900 nm range of interest for atomic spectroscopy, and also few $\mu$W per
mode in the Telecom bands. The comb is pumped with standard 1064 nm laser
technology, and we demonstrate high coherence with the soliton microcomb
through detection of its 1 THz repetition frequency and $f$-$2f$ measurements.

##  II Soliton Microcomb Design

Figure 2: Soliton microcomb spectrum control through nanofabrication. (a)
Single-soliton spectra obtained from resonators with varying waveguide widths
on the same chip. Relevant atomic transitions are indicated on the plot. (b)
Comparison between measured dispersive wave frequencies and that predicted by
simulated waveguide dispersion show good agreement within fabrication
imprecision of resonator dimensions. (c) A zoom in plot near the rubidium
optical transition frequencies.

The photonic-chip resonators in this study are SiN microring resonators
embedded in a silicon-dioxide cladding. A ring radius of 22.5 $\mu$m (center
of waveguide) is chosen to give a free-spectral range (FSR) of $\approx$ 1
THz, which reduces the comb operating power Briles  _et al._ (2018a). The
devices are fabricated at wafer scale by Ligentec, using LPCVD stoichiometric
silicon-nitride deposition, DUV stepper lithography, dry chemical etching, and
separation to individual chips. Broadband, weakly anomalous GVD is created
near the 1064 nm pump laser wavelength by design of the resonator dimensions,
namely the SiN layer thickness and resonator waveguide width (RW). GVD
simulations for resonators with SiN layer thicknesses of 667 nm and RW = 810,
850 and 890 nm illustrate the soliton microcomb design process; see the top
panel of Fig. 1(a). We indicate the telecom bands (1500-1600 nm) and the NIR
range relevant for rubidium atomic spectroscopy (778-795 nm) with vertical
gray boxes.

Our resonators also support the generation of DWs that extend the soliton
spectrum outside the anomalous-GVD region. The emitted DWs correspond
approximately to the wavelengths of zero integrated dispersion, $D_{int}$
Brasch  _et al._ (2016). The integrated dispersion is equivalent to the
detuning of the cold-resonator frequencies $\nu_{\mu}$ from an equidistant
grid centered about the pump laser and spaced by the free spectral range (FSR)
about the pump laser,

$D_{int}(\mu)=\nu_{\mu}-\left(\nu_{p}-\mu\times\frac{D_{1}}{2\pi}\right)$ (1)

Here, $\nu_{p}$ is the pump laser frequency, $\frac{D_{1}}{2\pi}$ is the
resonator FSR near the pump and $\mu$ is an integer indexing the resonances
relative to the pumped mode $\mu=0$. This quantity, given in units of GHz, is
plotted in the middle panel of Fig. 1(a) for the three corresponding traces in
the top panel. We note that increasing the resonator waveguide-width (RW)
parameter leads to an overall redshift of the GVD profile as well as the
wavelengths of $D_{int}=0$. Since both the device layer thickness and RW
affect the GVD profile, we can quasi-independently design and control the peak
wavelength of the two DWs.

The linear approximation for the DW wavelengths given by $D_{int}(\mu)=0$
neglects higher-order effects arising from soliton self-phase modulation
Skryabin and Kartashov (2017). However, it is a reasonable criteria for
predicting the performance of resonator designs given the fabrication
tolerances of $\sim$5 nm in device layer thickness and $\sim$20 nm in pattern
dimensions, including effects of waveguide sidewall angle. These small
deviations are accounted for by fabricating resonators of various RW values,
stepped by 10’s of nanometers, on a single chip. We design chips composed of
$\sim$75 individual ring resonators to systematically cover the microcomb
design space.

The bottom panel of Fig. 1(a) shows the result of controllably taking
advantage of the microcomb design space in SiN. The soliton microcomb spectrum
that we present features a long-wavelength DW in the telecom band and short-
wavelength DW with its strongest line aligned to the Rb 2-photon clock
transition at 778 nm Hilico  _et al._ (1998); Zameroski  _et al._ (2014);
Bernard  _et al._ (2000).

##  III Optical Measurement

The setup we use for our experiments is shown in Fig. 1(b). In particular, the
apparatus is composed of elements for soliton microcomb generation and thermo-
optic tuning, and for coherence measurements of the soliton microcomb’s
carrier-envelope-offset frequency and repetition frequency. The pump laser we
use is an external cavity diode laser, providing coarse tuning from 1020 to
1070 nm to address a few of the TE1 ring resonator modes that are spaced by 1
THz.

In order to provide laser light to our resonators, coupling waveguides are
fabricated in the silicon nitride device layer. Inverse tapering of the
coupling waveguide at chip edge enables approximately 3 dB per facet insertion
loss from standard lensed fibers. A ‘pulley’ design Hosseini  _et al._ (2010);
Li  _et al._ (2015); Briles  _et al._ (2018a) for the coupling waveguide
achieves adequate coupling rates for low-power operation of the comb and also
at the edges of the spectrum for output-coupling of the dispersive waves. We
include a variation in the coupling gap between the pulley coupler and the
resonator to achieve critical coupling, namely, to set the coupling strength
between the resonator and the coupling waveguide to match the internal loss of
the resonator. This coupling condition approximately maximizes the intra-
resonator power for a given pump power, reducing the pump power required to
support soliton microcombs. We extract information on the intrinsic and
coupled quality factor by measuring the resonance linewidth and minimum
transmission as a function of gap. The spectra shown in Figs 1(a) and 2 were
obtained using resonators with intrinsic quality factors of $\approx 1\times
10^{6}$ and comb generation threshold powers as low as 9 mW in the coupling
waveguide. Ongoing experiments show that intrinsic Q’s up to $\approx 3\times
10^{6}$ and threshold powers of 1 mW are possible at 1064 nm.

We efficiently generate soliton microcombs by controlling the pump-resonator
detuning at a rate faster than the resonator thermal dynamics Briles  _et al._
(2018b); Stone  _et al._ (2018); Li  _et al._ (2017). This is implemented
using an ultrafast electro-optic frequency shifting technique similar to that
reported in Ref Briles  _et al._ (2018a). Here we utilize a standard lithium-
niobate waveguide Mach-Zehnder intensity modulator in place of the single-
sideband, suppressed-carrier modulator. The output of the intensity modulator
is passed to a polarization controller and amplified in an ytterbium-doped
fiber amplifier (YDFA) before being coupled into the chip. The blue-detuned
sideband frequency is tuned to resonance and then decreased by $\approx 5$ GHz
in 100 ns. This rapid pump frequency tuning technique minimizes the thermo-
optic shift of the resonance frequency, which can destabilize the soliton. The
red sideband and residual carrier are far detuned (by $\geq$20 resonator
linewidths) from any resonances and do not adversely affect soliton
generation. We filter out the pump light to monitor the generated comb power.
The emergence of soliton microcombs from chaotic states as we vary the
detuning is identified from the photodetected comb power, using an
oscilloscope; see the trace plotted in Fig. 1(b).

Optical spectra of single-soliton microcombs that we create are shown in Fig.
2(a). We present spectra for five RW values, spaced by 20 nm from RW = 810-890
nm. Each spectrum features two prominent dispersive waves on either side of
the pump laser. The measured DW frequencies in the telecom C/L band and the
NIR 780 nm band are plotted versus RW in Fig. 2(b) as hollow red and blue
diamonds, respectively. For comparison, we also present DW frequency
measurements obtained with blue-detuned chaotic combs (solid circles). This
thorough and unique comparison of microcomb spectra is enabled by precision
fabrication of numerous devices and reliable soltion microcomb generation. Our
measurements of the DW frequencies agree with our designs, which are based on
finite-element calculations of $D_{int}$ (open circles), using accurate
Sellmeier coefficients for SiN Luke  _et al._ (2015) and silica Malitson
(1965) and a fabrication precision that we assess to be $\sim$20 nm. The
predictability of dispersive wave frequencies enables placement of dispersive
waves on the targeted Rb resonances, as shown in the zoom-in plot in Fig.
3(c).

Figure 3: Soliton coherence measurements. (a) The filtered and amplified long-
wavelength DW before (black) and after (red) being sent to electro-optic phase
modulators. A heterodyne beat is made at the wavelength marked with an arrow
to evaluate fluctuations in THz comb mode spacing. This beat is shown in (b).
(c,d) Heterodyne beats of each DW with a low-noise auxiliary laser: (c) LDW at
1563 nm and (d) SDW at 777.3 nm. The traces (b) through (d) are taken with
resolution bandwidth 1 kHz, 100 kHz, and 300 kHz, respectively.

We measure the repetition frequency of the soliton microcomb by electro-optic
modulation of the LDW near 1550 nm. Such measurements are shown for the
RW=0.83 $\mu$m device in Fig. 3. This procedure creates a low frequency $<$100
MHz optical heterodyne beatnote $f_{beat}$ that represents the 0.999291 THz
repetition frequency. To facilitate these measurements, we use optical filters
to isolate the dispersive waves near 1550 nm and 780 nm. The long-wavelength
DW is amplified with a C-band erbium-doped fiber amplifier (EDFA), filtered,
and sent to two electro-optic phase modulators operating in series and driven
at $f_{eo}\approx 16.65$ GHz to create many coherent sidebands around each
comb tooth Del’Haye  _et al._ (2012). We detect $f_{beat}$ by photodetecting
the phase modulation sidebands nearly in between the two soliton microcomb
modes, allowing the repetition frequency to be determined via $f_{rep}=60\cdot
f_{eo}+f_{beat}$. A narrow optical bandpass filter selects the desired
interfering modes (marked with an arrow in Fig. 3(a)) before photodetecting
$f_{beat}$. The highly coherent beat obtained from this measurement (shown in
Fig. 3(b)) confirms low-noise soliton operation.

We evaluate the frequency noise on individual comb modes by recording optical
heterodyne beats between the soliton comb and an auxiliary low-noise CW laser.
We focus on large comb mode numbers at a large frequency offset from the 1064
nm pump laser, which are interesting to test since they are most sensitive to
multiplicative noise of the comb. A C-band external cavity diode laser is
heterodyned directly with the long-wavelength DW near 1563 nm, and it is also
frequency doubled in a bulk periodically poled lithium niobate crystal for the
short-wavelength DW near 777.3 nm. For both heterodyne beatnotes, which are
shown in Figs. 3 (c) and (d) for the long- and short-wavelength DW,
respectively, we observe linewidths of a few megahertz that we associate with
the linewidth of the soliton microcomb modes. This is reasonable considering
that the free-running repetition rate has a frequency linewidth of $\approx$5
kHz and the comb modes we characterize in Figs. 3 (c) and (d) are $\approx$100
modes from the pump.

We note that the DW wavelengths were chosen in this case to be near-harmonic
to facilitate $f$-$2f$ referencing. For this particular device, the offset
frequency $f_{0}$ is $\approx 80\pm 1$ GHz. While such a high frequency makes
a direct electronic measurement of $f_{0}$ challenging, in future designs we
would systematically vary the resonator radius in steps of $\approx 10$ nm to
create a specific device that satisfies the conditions $f_{0}\leq 20$ GHz
Briles  _et al._ (2018a).

##  IV Thermo-optic Fine Tuning

We finely shift the absolute frequencies of the soliton comb modes beyond the
level provided by nanofabrication process control, using thermo-optic tuning
of the resonator. This approach combines the high optical power per mode of
DKS combs with the high spectral resolution typical of lower repetition rate
combs Heinecke  _et al._ (2009). The temperature of the chip is set with a
PID-controlled thermoelectric heater; see Fig. 1 (b) for the setup. Fig 4 (a)
and (b) analyze the thermal tuning of solitons created in RW = 0.87 $\mu$m
resonators (second to last panel in Fig. 2 (a) and (c)). Comb modes near 778
nm have a measured tuning coefficient of -5.48 GHz/ ∘C which is sufficiently
large to access the Rb $2h\nu$ and D2 resonances at chip temperatures of
$\approx 37$ and $\approx 44$ ∘C respectively; see Fig. 4(a). The change in
the comb modes power between these two temperatures is associated with
hysteresis in switching behavior between distinct operating regimes of the
single soliton pulse, which we plan to explore in more detail in future work.
The measured linear thermal tuning of the soliton repetition rate of $15.86\pm
0.34$ MHz $/^{\circ}C$ shown in Fig. 4 (b) amounts to a fractional change of
$\approx 16$ ppm on the $1.0$ THz spacing. This is in reasonable agreement
with measured thermo-optic coefficients of bulk silicon nitride at 20 to 30
ppm Ikeda  _et al._ (2008)as well as experimental measurements on resonators
pumped at 1550 nm Xue  _et al._ (2016). In future work we plan to incorporate
integrated heaters Xue  _et al._ (2016); Joshi  _et al._ (2016), which will
reduce the overall electrical power consumption.

Figure 4: Thermo-optic tuning of soliton spectrum. (a) The thermal tuning of
two modes of the soliton spectrum across the Rb two-photon and Rb D2
transitions. (b) The repetition rate of the soliton can also be thermally
tuned. (c) Overlay of the spectra presented in Fig. 2 plotted over the range
from 333 to 337 THz. Note that any arbitrary frequency can be reached with the
current RW sweep and a 50 ∘C thermal tuning range.

Through both lithographic variation in resonator dimensions and thermo-optic
tuning, we exert comprehensive frequency control over the SiN comb. A typical
chip used in this study contains 70 devices, which is sufficient to include
combinations of fine resonator width and gap sweeps for desired dispersion and
coupling. A modest $\Delta$T = $50^{\circ}C$ thermal-optic tuning of the chip
is sufficient to place a comb line onto any desired wavelength within the
design bandwidth of the resonators. For example, to probe the cesium D1
transition, we can select a device with ring width 870 nm or 890 nm, and
thermally tune the modes toward lower frequency to align to the desired
transition, as shown in Fig. 4(c).

##  V Conclusion

We have demonstrated optical-frequency-comb generation in the near infrared,
using silicon nitride ring resonators. These devices demonstrate octave-span
comb generation that can address optical resonances of various atomic species.
Sufficient power output in the far short- and long-wavelength sides of the
comb spectra is achieved utilizing engineered dispersive waves and efficient
out-coupling waveguides. Single Kerr solitons are reliably generated, using a
fast frequency sweep technique that works with all devices presented in this
study, and the coherent properties of the soliton microcombs are verified by
heterodyne beatnote measurements.

From the perspective of building a capable frequency-comb system, our 1064 nm
pumped soliton microcombs could be combined with mature optical amplifier,
optical modulator, and integrated second-harmonic generation technology in the
1550 nm band. These tools present a route to $f$-$2f$ detection of our combs
that are frequency tunable continuously across the near infrared. Indeed, we
have demonstrated repetition rate detection of the soliton microcomb with
phase modulation. Further, we have demonstrated the flexibility of our system
by thermally tuning optical comb modes onto Rb and Cs resonances, which would
be enabling for future experiments that explore spectroscopy and metrology
with atomic resonances.

## Funding Information

Funding provided by DARPA ACES and DODOS, Air Force Office of Scientific
Research (FA9550-16-1-0016), NIST, NASA, and NRC.

## Acknowledgments

We thank Jordan R. Stone and Connor Fredrick for a careful reading of the
manuscript and Ligentec for fabrication of the optical resonators. This work
is a contribution of the U.S. Government and is not subject to copyright.
Mention of specific companies or trade names is for scientific communication
only, and does not constitute an endorsement by NIST.

## References

  * Moss  _et al._ (2013) D. J. Moss, R. Morandotti, A. L. Gaeta, and M. Lipson, Nature Photonics 7, 597 (2013).
  * Li  _et al._ (2013) J. Li, H. Lee, and K. J. Vahala, Nature Communications 4, 2097 (2013).
  * Loh  _et al._ (2015) W. Loh, A. A. Green, F. N. Baynes, D. C. Cole, F. J. Quinlan, H. Lee, K. J. Vahala, S. B. Papp, and S. A. Diddams, Optica 2, 225 (2015).
  * Carlson  _et al._ (2017) D. R. Carlson, D. D. Hickstein, A. Lind, S. Droste, D. Westly, N. Nader, I. Coddington, N. R. Newbury, K. Srinivasan, S. A. Diddams, _et al._ , Optics letters 42, 2314 (2017).
  * Lamb  _et al._ (2018) E. S. Lamb, D. R. Carlson, D. D. Hickstein, J. R. Stone, S. A. Diddams, and S. B. Papp, Physical Review Applied 9, 024030 (2018).
  * Kippenberg  _et al._ (2018) T. J. Kippenberg, A. L. Gaeta, M. Lipson, and M. L. Gorodetsky, Science 361, eaan8083 (2018).
  * Bao  _et al._ (2016) C. Bao, J. A. Jaramillo-Villegas, Y. Xuan, D. E. Leaird, M. Qi, and A. M. Weiner, Physical Review Letters 117, 163901 (2016).
  * Lucas  _et al._ (2017) E. Lucas, M. Karpov, H. Guo, M. Gorodetsky, and T. J. Kippenberg, Nature Communications 8, 736 (2017).
  * Yu  _et al._ (2017) M. Yu, J. K. Jang, Y. Okawachi, A. G. Griffith, K. Luke, S. A. Miller, X. Ji, M. Lipson, and A. L. Gaeta, Nature Communications 8, 14569 (2017).
  * Xue  _et al._ (2015) X. Xue, Y. Xuan, Y. Liu, P.-H. Wang, S. Chen, J. Wang, D. E. Leaird, M. Qi, and A. M. Weiner, Nature Photonics 9, 594 (2015).
  * Cole  _et al._ (2017) D. C. Cole, E. S. Lamb, P. Del’Haye, S. A. Diddams, and S. B. Papp, Nature Photonics 11, 671 (2017).
  * Marin-Palomo  _et al._ (2017) P. Marin-Palomo, J. N. Kemal, M. Karpov, A. Kordts, J. Pfeifle, M. H. Pfeiffer, P. Trocha, S. Wolf, V. Brasch, M. H. Anderson, R. Rosenberger, K. Vijayan, W. Freude, T. J. Kippenberg, and C. Koos, Nature 546, 274 (2017).
  * Spencer  _et al._ (2018) D. T. Spencer, T. Drake, T. C. Briles, J. Stone, L. C. Sinclair, C. Fredrick, Q. Li, D. Westly, B. R. Ilic, A. Bluestone, N. Volet, T. Komljenovic, L. Chang, S. H. Lee, D. Y. Oh, M.-G. Suh, K. Y. Yang, M. H. P. Pfeiffer, T. J. Kippenberg, E. Norberg, L. Theogarajan, K. Vahala, N. R. Newbury, K. Srinivasan, J. E. Bowers, S. A. Diddams, and S. B. Papp, Nature 557, 81 (2018).
  * Yi  _et al._ (2015) X. Yi, Q.-F. Yang, K. Y. Yang, M.-G. Suh, and K. Vahala, Optica 2, 1078 (2015).
  * Wang  _et al._ (2017) Y. Wang, F. Leo, J. Fatome, M. Erkintalo, S. G. Murdoch, and S. Coen, Optica 4, 855 (2017).
  * Matsko  _et al._ (2016) A. B. Matsko, W. Liang, A. A. Savchenkov, D. Eliyahu, and L. Maleki, Optics letters 41, 2907 (2016).
  * Brasch  _et al._ (2016) V. Brasch, M. Geiselmann, T. Herr, G. Lihachev, M. H. Pfeiffer, M. L. Gorodetsky, and T. J. Kippenberg, Science 351, 357 (2016).
  * Jaramillo-Villegas  _et al._ (2015) J. A. Jaramillo-Villegas, X. Xue, P.-H. Wang, D. E. Leaird, and A. M. Weiner, Optics Express 23, 9618 (2015).
  * Briles  _et al._ (2018a) T. C. Briles, J. R. Stone, T. E. Drake, D. T. Spencer, C. Fredrick, Q. Li, D. Westly, B. Ilic, K. Srinivasan, S. A. Diddams, and S. B. Papp, Optics Letters 43, 2933 (2018a).
  * Stern  _et al._ (2018) B. Stern, X. Ji, Y. Okawachi, A. L. Gaeta, and M. Lipson, Nature , 1 (2018).
  * Papp  _et al._ (2014) S. B. Papp, K. Beha, P. Del’Haye, F. Quinlan, H. Lee, K. J. Vahala, and S. A. Diddams, Optica 1, 10 (2014).
  * Karpov  _et al._ (2018) M. Karpov, M. H. P. Pfeiffer, J. Liu, A. Lukashchuk, and T. J. Kippenberg, Nature Communications 9 (2018), 10.1038/s41467-018-03471-x.
  * Agha  _et al._ (2009) I. H. Agha, Y. Okawachi, and A. L. Gaeta, Optics Express 17, 16209 (2009).
  * Okawachi  _et al._ (2014) Y. Okawachi, M. R. Lamont, K. Luke, D. O. Carvalho, M. Yu, M. Lipson, and A. L. Gaeta, Optics Letters 39, 3535 (2014).
  * Diddams  _et al._ (2004) S. A. Diddams, J. C. Bergquist, S. R. Jefferts, and C. W. Oates, Science 306, 1318 (2004).
  * Ludlow  _et al._ (2015) A. D. Ludlow, M. M. Boyd, J. Ye, E. Peik, and P. O. Schmidt, Reviews of Modern Physics 87, 637 (2015).
  * Targat  _et al._ (2013) R. L. Targat, L. Lorini, Y. L. Coq, M. Zawada, J. Guéna, M. Abgrall, M. Gurov, P. Rosenbusch, D. G. Rovera, B. Nagórny, R. Gartman, P. G. Westergaard, M. E. Tobar, M. Lours, G. Santarelli, A. Clairon, S. Bize, P. Laurent, P. Lemonde, and J. Lodewyck, Nature Communications 4 (2013), 10.1038/ncomms3109.
  * Skryabin and Kartashov (2017) D. Skryabin and Y. Kartashov, Optics Express 25, 27442 (2017).
  * Hilico  _et al._ (1998) L. Hilico, R. Felder, D. Touahri, O. Acef, A. Clairon, and F. Biraben, The European Physical Journal-Applied Physics 4, 219 (1998).
  * Zameroski  _et al._ (2014) N. D. Zameroski, G. D. Hager, C. J. Erickson, and J. H. Burke, Journal of Physics B: Atomic, Molecular and Optical Physics 47, 225205 (2014).
  * Bernard  _et al._ (2000) J. Bernard, A. Madej, K. Siemsen, L. Marmet, C. Latrasse, D. Touahri, M. Poulin, M. Allard, and M. Têtu, Optics Communications 173, 357 (2000).
  * Hosseini  _et al._ (2010) E. S. Hosseini, S. Yegnanarayanan, A. H. Atabaki, M. Soltani, and A. Adibi, Optics Express 18, 2127 (2010).
  * Li  _et al._ (2015) Q. Li, T. C. Briles, D. Westly, J. Stone, R. Ilic, S. Diddams, S. Papp, and K. Srinivasan, in  _Frontiers in Optics_ (Optical Society of America, 2015) pp. FW6C–5.
  * Briles  _et al._ (2018b) T. C. Briles, J. R. Stone, T. E. Drake, D. T. Spencer, S. A. Diddams, and S. B. Papp, In preparation (2018b).
  * Stone  _et al._ (2018) J. R. Stone, T. C. Briles, T. E. Drake, D. T. Spencer, D. R. Carlson, S. A. Diddams, and S. B. Papp, Physical Review Letters 121, 063902 (2018).
  * Li  _et al._ (2017) Q. Li, T. C. Briles, D. A. Westly, T. E. Drake, J. R. Stone, B. R. Ilic, S. A. Diddams, S. B. Papp, and K. Srinivasan, Optica 4, 193 (2017).
  * Luke  _et al._ (2015) K. Luke, Y. Okawachi, M. R. E. Lamont, A. L. Gaeta, and M. Lipson, Optics Letters 40 (2015), 10.1364/OL.40.004823.
  * Malitson (1965) I. H. Malitson, Journal of the Optical Society of America 55, 1205 (1965).
  * Del’Haye  _et al._ (2012) P. Del’Haye, S. B. Papp, and S. A. Diddams, Physical Review Letters 109, 263901 (2012).
  * Heinecke  _et al._ (2009) D. C. Heinecke, A. Bartels, T. M. Fortier, D. A. Braje, L. Hollberg, and S. A. Diddams, Physical Review A 80, 053806 (2009).
  * Ikeda  _et al._ (2008) K. Ikeda, R. E. Saperstein, N. Alic, and Y. Fainman, Optics Express 16, 12987 (2008).
  * Xue  _et al._ (2016) X. Xue, Y. Xuan, C. Wang, P.-H. Wang, Y. Liu, B. Niu, D. E. Leaird, M. Qi, and A. M. Weiner, Optics Express 24, 687 (2016).
  * Joshi  _et al._ (2016) C. Joshi, J. K. Jang, K. Luke, X. Ji, S. A. Miller, A. Klenner, Y. Okawachi, M. Lipson, and A. L. Gaeta, Optics Letters 41, 2565 (2016).


11footnotetext: Work partially supported by U.S. National Science Foundation
grant DMS 2055282.

# A monotonicity version of a concavity theorem of Lieb

Eric A. Carlen  $\\\
\vspace{5pt}\small{DepartmentofMathematics,HillCenter,}\\\
\small{RutgersUniversity,110FrelinghuysenRoadPiscatawayNJ08854-8019USA}\\\ $

###### Abstract

We give a simple proof of a strengthened version of a theorem of Lieb that
played a key role in the proof of strong subadditivity of the quantum entropy.

##  1 Introduction

We write $M_{n}({\mathord{\mathbb{C}}})$ to denote the $n\times n$ complex
matrices, $M_{n}^{+}({\mathord{\mathbb{C}}})$ the subset consisting of
positive semidefinite matrices, and $M_{n}^{++}({\mathord{\mathbb{C}}})$ the
subset consisting of positive definite matrices. The following theorem was
proved by Lieb in 1973 [7, Theorem 6].

######  1.1 THEOREM (Lieb).

For all self-adjoint $H\in M_{n}({\mathord{\mathbb{C}}})$ the function

$Y\mapsto{\rm Tr}\left[\exp(H+\log Y)\right]$ (1.1)

is concave on $M_{n}^{++}({\mathord{\mathbb{C}}})$.

As a simple consequence of this, Lieb deduced his triple matrix inequality, a
generalization of the Golden-Thompson inequality to three self adjoint
matrices. This played a fundamental role in the proof of strong subadditivity
of the quantum entropy [8]. For more recent applications of Theorem 1.1, see
the influential paper of Tropp [14].

We now prove a stronger version of this theorem in terms of monotonicity
instead of concavity. Recall that a linear map
$\Phi:M_{n}({\mathord{\mathbb{C}}})\to M_{m}({\mathord{\mathbb{C}}})$ is
positive if $\Phi(A)\in M_{m}^{+}({\mathord{\mathbb{C}}})$ whenever $A\in
M_{n}^{+}({\mathord{\mathbb{C}}})$, and is unital if $\Phi(I)=I$, and trace
preserving if ${\rm Tr}[\Phi(X)]={\rm Tr}[X]$ for all $X\in
M_{n}({\mathord{\mathbb{C}}})$. We equip $M_{n}({\mathord{\mathbb{C}}})$ with
the Hilbert-Schmidt inner product, and we use $\Phi^{\dagger}$ to denote the
corresponding adjoint of a linear map on $M_{n}({\mathord{\mathbb{C}}})$. Note
that $\Phi$ is unital if and only if $\Phi^{\dagger}$ is trace preserving. Our
main result is:

######  1.2 THEOREM.

Let $\Phi:M_{n}({\mathord{\mathbb{C}}})\to M_{m}({\mathord{\mathbb{C}}})$ be
unital and positive. Then for all self-adjoint $H\in
M_{n}({\mathord{\mathbb{C}}})$ and all $Y\in
M_{m}^{++}({\mathord{\mathbb{C}}})$,

${\rm Tr}\left[\exp(H+\log\Phi^{\dagger}(Y))\right]\geq{\rm
Tr}\left[\exp(\Phi(H)+\log Y)\right]\ .$ (1.2)

Before giving the very simple proof, we explain how Theorem 1.2 implies
Theorem 1.1. Let $\Phi:M_{n}({\mathord{\mathbb{C}}})\to
M_{2n}({\mathord{\mathbb{C}}})$ be defined by

$\Phi(H)=\left[\begin{array}[]{cc}H&0\\\ 0&H\end{array}\right]\quad{\rm and\
hence}\quad\Phi^{\dagger}\left(\left[\begin{array}[]{cc}Y_{1}&0\\\
0&Y_{2}\end{array}\right]\right)=Y_{1}+Y_{2}\ .$ (1.3)

Thus,

${\rm Tr}[e^{H+\log(Y_{1}+Y_{2})}]\geq{\rm Tr}[e^{H+\log Y_{1}}]+{\rm
Tr}[e^{H+\log Y_{2}}]\ .$

Since $Y\mapsto{\rm Tr}[e^{H+\log Y}]$ is homogeneous of degree one, this is
the same as concavity. Note that this particular map $\Phi$ is not only
positive; it is completely positive.

Our proof is based on the well-known and elementary Gibbs variational
principle for the free energy in terms of the entropy $S(X)=-{\rm Tr}[X\log
X]$ of a density matrix $X$. This states that for all self-adjoint $K\in
M_{n}({\mathord{\mathbb{C}}})$

$\log({\rm Tr}[e^{K}])=\sup\\{{\rm Tr}[XK]-{\rm Tr}[X\log X]\ :\ X\in
M_{n}^{++}({\mathord{\mathbb{C}}})\ ,{\rm Tr}[X]=1\ \\}\ .$ (1.4)

A short proof from scratch can be found in Appendix A of [5]. There is a
simple variant involving the relative entropy

$D(X||Y)={\rm Tr}[X(\log X-\log Y)]$

of two density matrices: For $W\in M_{n}^{++}({\mathord{\mathbb{C}}})$ replace
$K$ with $K+\log W$ in (1.4) to conclude that for all self-adjoint $K\in
M_{n}({\mathord{\mathbb{C}}})$ and all $W\in
M_{m}^{++}({\mathord{\mathbb{C}}})$,

$\log({\rm Tr}[e^{K+\log W}])=\sup\\{{\rm Tr}[XK]-D(X||W)\ :\ X\in
M_{n}^{++}({\mathord{\mathbb{C}}})\ ,{\rm Tr}[X]=1\ \\}\ .$ (1.5)

We shall also use the result due to Müller-Hermes and Reeb [10] that the
relative entropy is monotone under positive trace-preserving maps
$\Phi^{\dagger}$. That is, for all such maps $\Phi^{\dagger}$, and all density
matrices $X,Y$

$D(\Phi^{\dagger}(X)||\Phi^{\dagger}(Y))\leq D(X||Y)\ .$ (1.6)

A somewhat weaker result is due to Uhlmann, who proved in 1977 that (1.6) is
true whenever $\Phi^{\dagger}$ is the dual of a unital Schwarz map; i.e., a
unital map $\Phi$ such that $\Phi(A^{*}A)\geq\Phi(A)^{*}\Phi(A)$. Earlier
still in 1973, Lindblad had shown that (1.6) was valid unital completely
positive $\Phi$. Lindblad’s proof relied on the Lieb Concavity Theorem [7]. It
is well known that every completely positive unital map is a Schwarz map, and
it is evident that Schwarz maps are positive. By now, quite simple proofs of
the monotonicity of the relative entropy under the dual of a unital Schwarz
map are known; see e.g., [11]. The deep result in [10] turns on ideas
introduced by Beigi [2].

###### Proof of Theorem 1.2.

By (1.5),

$\displaystyle\log({\rm Tr}[e^{H+\log\Phi^{\dagger}(Y)}])$ $\displaystyle=$
$\displaystyle\sup\\{{\rm Tr}[XH]-D(X||\Phi^{\dagger}(Y))\ :\ X\in
M_{n}^{++}({\mathord{\mathbb{C}}})\ ,{\rm Tr}[X]=1\ \\}$ $\displaystyle\geq$
$\displaystyle\sup\\{{\rm
Tr}[\Phi^{\dagger}(W)H]-D(\Phi^{\dagger}(W)||\Phi^{\dagger}(Y))\ :\ W\in
M_{m}^{++}({\mathord{\mathbb{C}}})\ ,{\rm Tr}[W]=1\ \\}$ $\displaystyle\geq$
$\displaystyle\sup\\{{\rm Tr}[W\Phi(H)]-D(W||Y)\ :\ W\in
M_{m}^{++}({\mathord{\mathbb{C}}})\ ,{\rm Tr}[W]=1\ \\}$ $\displaystyle=$
$\displaystyle\log{\rm Tr}\left[e^{\Phi(H)+\log Y}\right]\ .$

where we used $\\{\ \Phi^{\dagger}(W)\ :\ W\in
M_{m}^{++}({\mathord{\mathbb{C}}})\ ,{\rm Tr}[W]=1\ \\}\subset\\{\ X\in
M_{n}^{++}({\mathord{\mathbb{C}}})\ ,{\rm Tr}[X]=1\ \\}$ and (1.6).
Exponentiating both sides of the inequality yields (1.2). ∎

While Theorem 1.2 is strictly stronger than Theorem 1.1, the greatest interest
in it may lie in the very simple proof that its proof provides of Theorem 1.1.
This is an interesting example of how it may be easiest to prove a concavity
result by first proving a monotonicity result, and then applying that to the
particular map $\Phi$ that is defined in (1.3).

It is interesting to observe another advantage of the monotonicity approach to
convexity or concavity inequalities. A duality method for proving convexity
and concavity inequalities was introduced by myself and Lieb in [4] which uses
the lemma:

######  1.3 LEMMA.

If $f(x,y)$ is jointly convex in $x,y$, then $g(x):=\inf_{y}f(x,y)$ is convex.
If $f(x,y)$ is jointly concave in $x,y$, then $g(x):=\sup_{y}f(x,y)$ is
concave.

This may be found in [12, Theorem 1], and the simple proof is also given in
[4]. Since $(X,W)\mapsto D(X||W)$, is jointly convex, for fixed $K$,
$(X,W)\mapsto{\rm Tr}[XK]-D(X||W)$, is jointly concave. Then by Lemma 1.3 and
(1.5), $W\mapsto\log{\rm Tr}[e^{K+\log W}]$ is concave. However, this is a
weaker statement than Theorem 1.1 since if $g(x)=\log(f(x))$ with $f$ positive
and twice continuously differentiable.

$g^{\prime\prime}(x)=-\left(\frac{f^{\prime}(x)}{f(x)}\right)^{2}+\frac{1}{f(x)}f^{\prime\prime}(x)\
.$

Thus, concavity of $f$ implies concavity of $\log f$, but not the other way
around. However, monotonicity of $f$ is equivalent to monotonicity of $\log
f$. For this reason, we could simply use the Gibbs variational principle to
prove our Theorem 1.2. Tropp [13] found an ingenious variational
representation of $W\mapsto{\rm Tr}[e^{K+\log W}]$ which allowed him to give a
proof of Theorem 1.1 using the joint convexity of the relative entropy and
Lemma 1.3, along the lines of [4]. However, when using duality to to prove the
monotonicity theorem, the logarithm is not an issue because log monotonicity
is the same as monotonicity, and we do not need Lemma 1.3. Already in 1973
Epstein [6] gave a second proof of Theorem 1.1 using the theory of Herglotz
functions, but this is considerably more involved than the present proof. For
more infromation, see [3].

Aclknowledgement I thank Aleksander Müller-Hermes for helpful correspondence
on his work with Reeb.

## References

  * [1]
  * [2] S. Beigi. Sandwiched Rényi Divergence Satisfies Data Processing Inequality. J. Math. Phys., 54 (2013), 122202.
  * [3] E. A. Carlen, On some convexity and monotonicity inequalities of Elliott Lieb, arXiv:2202.03591
  * [4] E. A. Carlen and E. H. Lieb. A Minkowski-type trace inequality and strong subadditivity of quantum entropy II: Convexity and concavity. Lett. Math. Phys., 83 (2008)107-126,
  * [5] E. A. Carlen and E. H. Lieb, Some trace inequalities for exponential and logarithmic functions, Bull. Math. Sciences, 1, No. 2 (2019) 1950008
  * [6] H. Epstein On a concavity theorem of Lieb, Comm.. Math. Phys. 31 (1973), 317–327
  * [7] E. H. Lieb, Convex trace functions and the Wigner-Yanase-Dyson conjecture. Advances in Math. 11 (1973), 267–288.
  * [8] E. H. Lieb and M. B. Ruskai: Proof of the strong subadditivity of quantum-mechanical entropy, J. Math. Phys. 14, 1973, pp. 1938-1941.
  * [9] G. Lindblad, Expectations and entropy inequalities for finite quantum systems, Comm. Math. Phys. 39 (1974), 111–119.
  * [10] A. Müller-Hermes and D. Reeb. Monotonicity of the quantum relative entropy under positive maps. Annales Henri Poincaré, 18 (2107), 1777–1788.
  * [11] D. Petz, Monotonicity of the qunatum relative entropy revisited, Rev. Math. Phys. 15 (2003), 79-91.
  * [12] R. T. Rockafellar, Conjugate duality and optimization, Vol. 16, Regional conference series in applied mathematics, SIAM, Philadelphia, 1974
  * [13] J. Tropp, From joint convexity of quantum relative entropy to a concavity theorem of Lieb, Proceedings of the American Mathematical Society 140 (2012), 1757–1760
  * [14] J. Tropp, User-Friendly Tail Bounds for Sums of Random Matrices, Foundations of Computational Mathematics 12 (2012), 389–434


11institutetext:  Yunnan Observatories, Chinese Academy of Sciences, Kunming
650011, China<EMAIL_ADDRESS>22institutetext: Key Laboratory for the
Structure and Evolution of Celestial Objects, Chinese Academy of Sciences,
Kunming 650011 33institutetext: Key Laboratory of Dark Matter and Space
Astronomy, Purple Mountain Observatory, Chinese Academy of Science, Nanjing
210008, China<EMAIL_ADDRESS>44institutetext: Department of Astronomy,
Yunnan University, Kunming 650091, China<EMAIL_ADDRESS>55institutetext:
Los Alamos National Laboratory, Los Alamos, New Mexico, 87545, USA

# The structure of TeV-bright shell-type supernova remnants

Chuyuan Yang 1122 Siming Liu 33 Jun Fang 44 Hui Li 55

(Received / Accepted)

###### Abstract

Aims. Two-dimensional magnetohydrodynamic (MHD) simulations are used to model
the emission properties of TeV-bright shell-type supernova remnants (SNRs) and
to explore their nature.

Methods.  In the leptonic scenario for the TeV emission, the $\gamma$-ray
emission is produced via Inverse Compton scattering of background soft photons
by high-energy electrons accelerated by the shocks of the SNRs. Since the
energy density of the Cosmic Microwave Background (CMB) radiation and that of
the IR/Optical background photons are much higher than that of the photons
produced by the same high-energy electrons via the synchrotron process, in a
previous paper, we showed that the observed correlation between X-ray and TeV
brightness of SNR RX J1713.7-3946 can be readily explained with the assumption
that the energy density of relativistic electrons is proportional to that of
the magnetic field. The TeV emissivity is therefore proportional to the
magnetic field energy density and MHD simulations can be used to model the TeV
structure of such remnants directly. 2D MHD simulations for SNRs are then
performed under the assumption that the ambient interstellar medium is
turbulent with the magnetic field and density fluctuations following a
Kolmogorov-like power-law spectrum.

Results. (1) As expected, these simulations confirm early 1D and 2D modelings
of these sources, namely the hydrodynamical evolution of the shock waves and
amplification of magnetic field by Rayleigh-Taylor convective flows and by
shocks propagating in a turbulent medium; (2) We reproduce rather complex
morphological structure for $\gamma$-rays, for example, bright thin rim and
significant asymmetry, suggesting intrinsic variations of the source
morphology not related to the structure of the progenitor and environment; (3)
The observed radial profile of several remnants are well reproduced with an
ambient medium density of $0.1-1$ cm-3. An even lower ambient density leads to
a sharper drop of the TeV brightness with radius than what is observed near
the outer edge of these remnants.

Conclusions. In a turbulent background medium, we can reproduce the observed
characteristics of several shell-type TeV SNRs with reasonable parameters
except for a higher ambient density than that inferred from X-ray
observations.

###### Key Words.:

ISM: supernova remnants – MHD – plasmas – radiation mechanisms: non-thermal –
shock waves – turbulence

##  1 Introduction

It is generally believed that supernova remnants (SNRs) are the main
acceleration sites of Galactic cosmic rays. Since charged cosmic ray particles
escaping from SNRs are mixed in the Galactic magnetic fields and can not be
used to image their sources of origin directly, detection of high-energy
$\gamma$-rays from SNRs is very useful for testing whether SNRs are
responsible for the origin of the bulk of the cosmic rays observed at the
Earth. However, an open question is whether the observed high-energy
$\gamma$-ray emission from an SNR has a hadronic or a leptonic origin (e.g.,
Katz & Waxman, 2008; Yuan et al., 2012).

Nonthermal photons from SNRs in the radio to X-ray bands are produced by the
synchrotron radiation of relativistic electrons. The detection of synchrotron
X-rays in some SNRs shows that electrons can be accelerated efficiently to TeV
energies by shocks of SNRs (e.g. Koyama et al., 1995). Although several young
SNRs have been detected to emit GeV to TeV $\gamma$-rays (e.g. Abdo et al.,
2011; Abramowski et al., 2011), it is not clear whether these SNRs produce
high-energy protons or not. High-energy $\gamma$-ray photons can be produced
by either inverse-Compton(IC) process or $\pi^{0}$-decay corresponding to the
leptonic or hadronic origin, respectively. For example, $\gamma$-ray emission
from SNR RX J1713.7-3946 is likely produced by very energetic electrons
radiating through IC scattering against some photon background (Abdo et al.,
2011), and objects like the Tycho SNR likely emit GeV $\gamma$-rays via the
decay of neutral pions produced in nuclear collisions between shock
accelerated relativistic nuclei and those in the background plasma (Morlino &
Caprioli, 2012; Ackermann et al., 2013).

Recent X-ray and $\gamma$-TeV observations have revealed a group of TeV-bright
shell-type SNRs with SNR RX J1713.7-3946 as a prototype. Similar to SNR RX
J1713.7-3946, no thermal emission has been detected from SNRs RX J0852.0-4622
and HESS J1731-347; the $\gamma$-ray luminosity is only a factor of a few
lower than the X-ray luminosity; and they have a very hard spectrum in the
$\sim$ GeV energy range, which favor the leptonic origin for the GeV$-$TeV
emission with a low mean magnetic field of $\sim 10\ \mu$Gauss (Li et al.,
2011; Ellison et al., 2012; Yuan et al., 2012; Yang et al., 2014). Rapid
variations of X-ray filaments with a width of about 0.1 light year on a
timescale of $\sim 1$ year have been interpreted as due to radiative energy
loss of TeV electrons in strong ($\sim 1$ mG) fields (Uchiyama,et al., 2007).
The strong field may only exist in these filaments with a very small volume
filling factor, and alternative explanations for this variability with a weak
mean field have been proposed as well (Bykov et al., 2008; Liu et al., 2008).
Although thermal emission has been detected from SNR RCW 86 (Williams et al.,
2011), the characteristics of nonthermal emission from this remnant is very
similar to other TeV-bright shell-type SNRs (Lemoine-Goumard et al., 2012;
Yuan et al., 2014). These results imply a universal mechanism of particle
acceleration by shocks of SNRs (Yuan et al., 2012).

The nonthermal emission of these remnants has a well-defined shell structure
and there is evidence that the relative mean thickness of the emitting shell
varies from remnant to remnant (Aharonian et al., 2006, 2007). Previously 1D
modeling of particle acceleration in SNRs has been used to study the radial
brightness profile in the hadronic scenario (Berezhko & Völk, 2010). Detailed
study of the source structure with MHD simulations hasn’t been carried out.
Although these remnants have complicated appearance in X-rays and
$\gamma$-rays with significantly fluctuations in both the radial and azimuthal
directions, there is a good correlation between the X-ray and $\gamma$-ray
brightness (Acero et al., 2009). Magnetic fields play an important role in
collisionless astrophysical shocks and in the acceleration of charged
particles. They can be amplified by cosmic ray induced streaming instability
in the presence of efficient particle acceleration (e.g. Bell, 2004; Riquelme
& Spitkovsky, 2009). In a turbulent medium, they can also be amplified by
turbulent motion in the shock downstream(Giacalone & Jokipii, 2007). The
latter has been explored recently by Guo et al. (2012) who performed extensive
2D numerical simulations for SNR blast wave interacting with a turbulent
plasma background and found that the magnetic field can also be amplified by
Rayleigh-Taylor convective flows induced at the contact discontinuity of the
shock flows. In our previous paper (Yang & Liu, 2013), we have shown that the
observed correlation between X-ray and $\gamma$-ray brightness for SNR RXJ
1713.7-3946 suggests that the energy density of the accelerated electrons is
proportional to that of the magnetic field in the leptonic scenario for the
$\gamma$-ray emission. In this paper, we study the structure of these remnants
and the effect of magnetic inhomogeneity on the TeV profile and morphology by
using 2D MHD simulations.

Our 2D MHD simulations are based on the PLUTO code (Mignone et al., 2007). In
Section 2, we describe the MHD model and show the magnetic field structure. In
section 3, an algorithm is proposed to convert the 2D magnetic field structure
into a $\gamma$-ray image of the remnant and several simulations are done to
fit the azimuthal averaged radial profile of the TeV brightness of the
remnants mentioned above. Discussion and our conclusions are given in Section
4.

##  2  Model Description

The numerical model has been described in Yang & Liu (2013). For the sake of
completeness, we briefly summarize the key issues here. The dynamical
evolution of an SNR shock propagating into a turbulent ambient medium is
simulated with the time-dependent ideal MHD equations of mass, momentum, and
energy conservation:

$\partial_{t}\rho+\nabla\cdot(\rho\bf u)=0\,,$ (1) $\partial_{t}\rho{\bf
u}+\nabla\cdot\left(\rho{\bf uu-{BB\over 4\pi}}\right)+\nabla P^{\prime}=0\,,$
(2) $\partial E+\nabla\cdot\left[(E+P^{\prime}){\bf u-{B(u\cdot B)\over
4\pi}}\right]=0\,,$ (3)

and the induction equation:

$\partial_{t}{\bf B}+\nabla\times({\bf u\times B})=0\,,$ (4)

where $P^{\prime}=P+B^{2}/8\pi$ is the total pressure with thermal pressure
$P$ and magnetic pressure $B^{2}/8\pi$, and cgs units have been adopted in
this paper. $E$ is the total energy density:

$E=\frac{P}{\gamma-1}+\frac{1}{2}\rho u^{2}+\frac{B^{2}}{8\pi},$ (5)

$\rho$, $\bf u$, $\bf B$, and $\gamma=5/3$ are the plasma mass density, fluid
velocity, magnetic field, and the adiabatic index, respectively.

The numerical scheme to solve the ideal MHD equations was presented by Mignone
et al. (2007) in detail. We model the simulation in a two-dimensional
Cartensian coordinate ($x,y$) with uniform grids. The size of simulation
domain is chosen to be 40 pc $\times$ 40 pc to cover the extension of young
SNRs . The supernova blast wave is driven by the injection of internal energy
and ram pressure in a small circular region at the center of simulation box,
in which the density is assumed to be constant corresponding to a plateau
volume ($4\pi r_{\rm inj}^{3}/3$) with $r_{\rm inj}=0.5$ pc. The initial
magnetic field and density in the background plasma include an average
component and a turbulent component. We assume a constant average magnetic
field $B_{0}$ along the $x$ direction and a constant average gas density
$n_{0}$ (Guo et al., 2012).

Both density and magnetic fluctuations are generated by the assumption of a
Kolmogorov-like power-law spectrum of the form

$P(k)\propto\frac{1}{1+(kL)^{\Gamma}},$ (6)

where the spectral index $\Gamma$ depends on the dimensionality and equals
$8/3$ for 2D system. $k$ is the magnitude of the wavevector and $L=3\
\mathrm{pc}$ is the turbulence coherence length. The turbulence is generated
by summing a large number of discrete wave modes with random phases (Giacalone
& Jokipii, 1999). The random component of magnetic field is given by

$\displaystyle\delta\textbf{B}(x,y)=\sum^{N_{m}}_{n=1}$
$\displaystyle\sqrt{C_{B}2\pi k_{n}\Delta
k_{n}P_{B}(k_{n})}(\sin\theta_{n}\hat{x}-\cos\theta_{n}\hat{y})$ (7)
$\displaystyle\times\exp(i\cos\theta_{n}k_{n}x+i\sin\theta_{n}k_{n}y+i\phi_{n})$

where $P_{B}(k_{n})$ represents the power of the wave mode $n$ with wavenumber
$k_{n}$. The turbulence distributes randomly in propagation direction
$-1<\cos\theta_{n}<1$ and with a random phase $0<\phi_{n}<2\pi$. $C_{B}$ is a
normalization constant determined by $\langle\delta B^{2}\rangle=B_{0}^{2}$,
where $\langle\rangle$ represents averaging in space.

The density fluctuations satisfy the following probability distribution
(Giacalone & Jokipii, 2007; Burlaga & Lazarus, 2000):

$n(x,y)=n_{0}\exp(f_{0}+\delta f)$ (8)

where $f_{0}$ is a constant chosen to give the average density $n_{0}$ and the
description of $\delta f$ is similar to the turbulent part of the magnetic
field:

$\displaystyle\delta f(x,y)=$ $\displaystyle\sum^{N_{m}}_{n=1}\sqrt{C_{f}2\pi
k_{n}\Delta k_{n}P_{f}(k_{n})}$ (9)
$\displaystyle\times\exp(i\cos\theta_{n}k_{n}x+i\sin\theta_{n}k_{n}y+i\phi_{n})\,,$

with the normalization constant $C_{f}$ determined by $\langle(\delta
n)^{2}\rangle=0.4n_{0}^{2}$.

Guo et al. (2012) carried out extensive study of 2D MHD simulations of SNRs in
a turbulent medium. We first have a test run with $n_{0}=0.04$ cm-3,
$B_{0}=0.8\ \mu$G, the total injected internal energy $E_{SN}=2\times 10^{51}$
ergs, and the total mass of the ejector $M_{ej}=2M_{\odot}$. The left panel of
Figure 1 shows the magnetic field structure at 1600 years. The units of the
coordinates are pc and the color scale is the logarithmic of the magnetic
field in units of $\mu$G with the peak magnetic field reaching about 100
$\mu$G. The magnetic field near the shock front is mostly amplified by
turbulent motion in the downstream (Giacalone & Jokipii, 2007). The strongest
magnetic field however is produced by Rayleigh-Taylor convective flows near
the contact discontinuity. Since the ejector is uniform, the magnetic field is
very weak in the downstream of the reverse shock.

In order to construct the profile of $\sim$TeV emission, we adopt the leptonic
scenario where the $\gamma$-ray is produced via the IC of background soft
photons off high-energy electrons accelerated by the shocks. Besides the CMB
radiation, we also include an IR background photon component with a
temperature of $T=30$ K and an energy density of 1.2 eV cm-3 (Yang & Liu,
2013). The energy density of radio to X-ray band synchrotron component of the
SNRs is much lower than the energy density of the CMB radiation, the
synchrotron self-Comptonization process therefore can be ignored. Since the
background photon field is approximately uniform, the $\gamma$-ray emission
characteristics are then determined by the properties of energetic electrons.
Yang & Liu (2013) showed that the observed correlation between the X-ray and
$\gamma$-ray brightness of these SNRs implies that the electron energy density
is proportional to the magnetic field energy density, we then can describe the
electron spatial distribution simply as the following:

$\frac{dN(x,y)}{d\gamma_{e}}=CB(x,y)^{2}F_{\gamma_{e}},$ (10)

where $F_{\gamma_{e}}\sim\gamma_{e}^{\delta}\exp(-\gamma_{e}/\gamma_{\rm
max})$ is a power law energy distribution with an exponential cutoff and
$\gamma_{e}$ is the electron Lorentz factor. From the overall $\gamma$-ray
spectrum of SNR RX J1713.7-3946, we find $\delta=-2$ and $\gamma_{\rm
max}=1.3\times 10^{7}$ (Li et al., 2011; Yang & Liu, 2013). The
proportionality of the magnetic field and high-energy electron energy
densities then implies that the electron distribution only depends on the
magnetic field energy density and the structure of TeV emission is determined
by $B^{2}$. A recent study of TeV bright shell-type SNRs by Yang et al. (2014)
showed that the energy partition between energetic electrons above 1 GeV and
the magnetic field may be different for different SNRs with a tendency that
the total magnetic field energy within the volume enclosed by the front shock
of SNRs increases faster than that of energetic electrons as the SNRs evolve.
But the proportionality of these energy densities is inferred from the
correlation between the X-ray and $\gamma$-ray brightness in the leptonic
scenario, we therefore expect our treatment be valid for sources with a good
correlation between the X-ray and $\gamma$-ray brightness as far as the
spatial structure of the $\gamma$-ray emission is concerned.

##  3 Results

In the following, we use SNR RX J1713.7-3946 as a prototype to describe the
modeling. The remnant was first discovered in X-rays with the ROSAT all-sky
survey in 1996 (Pfeffermann & Aschenbach, 1996). It is close to the Galactic
plane and its distance and ages are about 1 kpc and 1600 yrs respectively
(e.g. Fukui et al., 2003; Cassam-Chenaï et al., 2004), and the radius of the
remnant is about 10 pc (Acero et al., 2009). The gamma-ray image of SNR RX
J1713.7-3946 is consistent with a thick emitting shell with a thickness close
to half of the radius and with the emissivity of southeast part is a factor of
two lower than the rest (Aharonian et al., 2006).

###  3.1 Structure and Radial Profile of TeV Emission

Figure 1: Left: The structure of the magnetic field $B$ in units of $\mu$G at
1600 years of a test run with $n_{0}=0.04$ cm-3. The rest of the model
parameters are shown in the text. Middle: TeV image derived from the 2D MHD
simulation. The image has been normalized to its peak value. Right: Comparison
of the observed radial profile averaged along the azimuthal direction with
results (solid line) derived from the 2D MHD simulation. The dotted line is
obtained by averaging the radial brightness profiles obtained in the position
angle range of $0.7\pi$ to $0.95\pi$ in the middle panel.

To model the observed morphology of the TeV emission, we also need to
construct a pseudo 3D structure with 2D simulations. One can integrate
$B^{2}$, which is proportional to the TeV emissivity, along a direction
(presumably corresponding to the line-of-sight) in the 2D simulation domain to
obtain a brightness profile in the radial direction of the SNR. Such a radial
profile can be obtained for all angular directions in the 2D domain. To mimic
the observed TeV image, we stack the radial profile obtained above
continuously in a polar coordinate and smooth the resultant 2D structure with
the point spread function of the HESS. The middle panel of Figure 1 shows such
an image derived from the structure of the magnetic field shown in the left
panel for an assumed distance of 1 kpc. The arcs connecting the center and the
outer rim are artifacts of the algorithm proposed above. Although this image
does not recover all details of the TeV image obtained with the HESS
(Aharonian et al., 2006), the strong variation of the brightness in the
azimuthal directions and the shell structure is well reproduced. The one-sided
structure of the simulated TeV image is closely related to the asymmetry of
the magnetic field structure. Due to the presence of a large scale magnetic
field in the horizontal direction, the magnetic field at the bottom and top
segments of the SNR shown in the left panel is preferentially amplified with
the bottom segment amplified the most. Since the TeV emissivity is
proportional to $B^{2}$, the asymmetry in the magnetic field structure is
further amplified in the TeV image giving rise to the bright southern rim in
the TeV image shown in the middle panel.

The right panel of Figure 1 compares the model predicted azimuthal averaged
radial profile with the observed profile. The radial profiles have been
normalized to their peak value. It can be seen that the solid line of the
model result is much narrower than the observed profile. This is not
surprising since the left panel of Figure 1 shows that the emission comes from
a shell with a thickness less than one-third of the radius while the observed
profile is consistent with an emitting shell with a thickness about half of
the radius (Aharonian et al., 2006). In the 1D model proposed by Berezhko &
Völk (2010), it is claimed that a broad radial profile can be produced by
considering the angular resolution of the observations. This 1D model predicts
a very bright rim at the shock front. The smoothing effect caused by the
angular resolution of the observations can reduce the brightness of this rim
dramatically and therefore suppress the brightness contrast between the center
and the shock front regions. In our 2D model, we don’t have a very bright rim
and the azimuthal variance also reduces the intrinsic brightness contrast of
the averaged radial profile. The angular resolution of the observations does
not change the brightness contrast between the center and the shock front
regions drastically.

Considering the complex structure of the magnetic field in the left panel
Figure 1 and its dependence on the large scale magnetic field, projection
effect may play a role in the brightness profile. Instead of averaging of the
whole $2\pi$ angle, the dotted line shows the profile averaged from $0.7\pi$
to $0.95\pi$. The model prediction is more consistent with observations than
the solid line. However, the brightness profile of the model still falls
faster than the observed profile near the edge of the remnant. 3D simulations
are needed to clarify this issue.

###  3.2 Application to TeV Images of SNRs

Table 1: Model Parameters.

Cases | $E_{\rm SN}$ | $M_{\rm ej}$ | $R_{s}$ | $n_{0}$ | Age | $V_{s}$ | $D$ | Ref. |
---|---|---|---|---|---|---|---|---|---
| (erg) | ($M_{\sun}$) | (pc) | ($\rm cm^{-3}$) | (year) | (km/s) | (kpc) |  |
RX J0852.0-4622 | $1\times 10^{51}$ | 1.4 | 8 | 0.1 | 990 | 5200 | 0.5 | 1-3 |
RX J1713.7-3946 | $2\times 10^{51}$ | 2 | 9 | 0.8 | 1600 | 3700 | 1.0 | 4-7 |
RCW 86 | $1\times 10^{51}$ | 1.4 | 11 | 0.4 | 1800 | 3900 | 1.9 | 8-11 |
HESS J1731-347 | $1\times 10^{51}$ | 2 | 13 | 0.1 | 2000 | 4300 | 3.0 | 12-15 |

  * Reference - (1) Aharonian et al. (2005); (2) Aharonian et al. (2007); (3) Bamba et al. (2005); (4) Fukui et al. (2003); (5) Cassam-Chenaï et al. (2004); (6) Acero et al. (2009); (7) Aharonian et al. (2006); (8) Aharonian et al. (2009); (9) Rosado, et al. (1996); (10) Bocchino et al. (2000); (11) Vink et al. (2006); (12) Tian, et al. (2008); (13) Aharonian et al. (2008); (14) Abramowski et al. (2011); (15) Bamba et al. (2012)

Figure 2: The structure of magnetic energy density $B^{2}$ for four SNRs
investigated in this paper. The parameters of model are described in Table 1.

Figure 3: Same as the middle panel of Figure 1 but for the four SNRs in Figure
2.  Figure 4: Same as the right panel of Figure 1 bur for the four SNRs in
Figure 2. The data of RX J0852.0-4622, RX J1713.7-3946, RCW86, and HESS
J1731-347 are from HESS observations (Aharonian et al., 2006, 2007, 2009;
Abramowski et al., 2011).

Besides SNR RX J1713.7-3946, SNR RX J0852.0-4622 is another very important
Galactic SNR. Its X-ray and gamma-ray images are consistent with a very thin
shell of $1^{\circ}$ in radius and $\sim 0.2^{\circ}$ in thickness. Its radius
is $5-10$ pc(Aharonian et al., 2005, 2007) for a distance of $0.26-0.50$ kpc,
and its age is about $420-1400$ yrs (Bamba et al., 2005). The young shell-type
SNR RCW 86 (also known as G315.4-2.3 and MSH 14-63) has been observed in radio
(Kesteven & Caswell, 1987), optical (Smith, 1997), X-ray (Pisarskietal et al.,
1984), and gamma-ray (Yuan et al., 2014; Aharonian et al., 2009) with a nearly
circular shape of $40^{\prime}$ in diameter. Rosado, et al. (1996) deduced a
distance of 2.8 kpc and an age of $\sim 10,000$ years by using optical
observations, whereas recent observations of the northeast part of the remnant
with the Chandra and XMM-Newton satellites strengthen the case that the event
recorded by the Chinese in 185 AD was a supernova and RCW 86 is likely its
remnant (Vink et al., 2006). Tian, et al. (2008) reported the discovery of a
faint radio and X-ray shell with a distance of $3.2\pm 0.8$ kpc closely
matching the extended TeV source HESS J1731-347 (Aharonian et al., 2008).
Deeper $\gamma$-ray observations of the source have revealed a large shell-
type structure with similar position and extension ($r\sim 0.25^{\circ}$) as
the radio SNR, confirming their association (Abramowski et al., 2011). These
authors suggested that its distance is not much larger than the derived low
limit of $3.2$ kpc, and a physical size of $\leq 15$ pc is also reasonable.

To improve the fitting to the azimuthal averaged radial profile of the TeV
brightness, one may adjust the mean density $n_{0}$ of the ambient medium. In
the following, we will adjust $n_{0}$ as the major free parameter to fit the
radial brightness profile of these SNRs. Considering the ambiguity of the
progenitor, we picked the mass and total energy of the ejector in a reasonable
region to reproduce the physical size of these SNRs. A fit to the observed
radial TeV brightness profile is then used to constrain the density $n_{0}$.
The model parameters are shown in Tab.1 including the output parameters, such
as the shock velocity $V_{s}$ and age.

Figure 2 shows the magnetic energy density. It is evident that the shell
thickness increases with the increase of the mean density of the background.
The results for RX J1713.7-3946 and RCW 86 have thicker shells than the other
two sources. Figures 3 and 4 show that the TeV morphology and radial profile
of these sources can be approximately reproduced with our model. To improve
this model, one may consider more realistic 3D simulations with the structure
of the progenitor wind and environment of these sources taken into account
properly. We note that the background density derived from our model is higher
than the upper limit inferred from the X-ray observations, which poses the
most serious challenge to our model. The thermal X-ray emission must be
suppressed in the shock downstream to be consistent with the absence or
weakness of thermal X-ray emission from these sources. Alternatively, in a 3D
simulation, a lower background density may be derived from the TeV radial
brightness profiles.

##  4 Discussion and Conclusions

Based on 2D MHD simulations of shell-type supernova remnants evolving in a
turbulent background, the evolution of the magnetic field for these remnants
can be obtained. In the leptonic scenario for the $\gamma$-ray emission, the
observed correlation between the X-ray and $\gamma$-ray brightnesses implies
that the energy density of high-energy electrons is proportional to that of
the magnetic field. One can then model the multi-band nonthermal emission with
the MHD simulations with the energy partition between energetic electrons and
the magnetic field as a free parameter. It is shown that the radial profile of
the TeV emission is very sensitive to the density of the background plasma.
With reasonable parameters for four TeV bright shell-type SNRs, we reproduce
the azimuthal averaged radial brightness profile of the TeV emission obtained
with the HESS observations.

Using multi-dimensional hydrodynamic simulations, Obergaulinger et al. (2014)
interpreted the inhomogeneities of the emission from the SNR RX J0852.0-4622
by assuming the blast wave evolved in a cloudy environment. To reproduce the
morphology of two opposite lobes for the SNR Kepler, Tsebrenko & Soker (2013)
proposed that jets existed before the explosion, which can produce the
axisymmetric structure with lobes. Our simulations are intrinsically 2D, and a
3D model including the details of both the ejecta and the background medium
around the TeV remnants is needed to reproduce the detailed structure of
nonthermal emission of these four remnants studied here.

###### Acknowledgements.

We thank Dr. Fan Guo for providing a code to generate turbulent ISM. This work
is partially support by the Strategic Priority Research Program C The
Emergence of Cosmological Structures of the Chinese Academy of Sciences, Grant
No. XDB09040200, NSFC grants: 11163006, 11173064, 11233001, and 11233008. Fang
is supported by the NSFC grant (11103016), and by the Key Project of Chinese
Ministry of Education (212160). HL is supported by the LDRD program at LANL.

## References

  * Abdo et al. (2011) Abdo, A. A. et al. 2011, ApJ, 734,28
  * Acero et al. (2009) Acero, F., Ballet, J., Decourchelle, A., et al. 2009, A&A, 505, 157
  * Ackermann et al. (2013) Ackermann, M., et al. 2013, Science, 339, 807-811
  * Aharonian et al. (2005) Aharonian, F. et al. 2005, A&A, 437,L7
  * Aharonian et al. (2006) Aharonian, F. et al. 2006, A&A, 449, 223
  * Aharonian et al. (2007) Aharonian, F. et al. 2007, ApJ,661,236
  * Aharonian et al. (2008) Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2008, A&A, 477, 353
  * Aharonian et al. (2009) Aharonian, F. et al. 2009, ApJ,692,1500
  * Abramowski et al. (2011) Abramowski A. et al. 2011, A&A, 531,1
  * Bamba et al. (2005) Bamba,A., Yamazaki,R. & Junko S. Hiraga 2005, ApJ,632,294
  * Bamba et al. (2012) Bamba,A. et al. 2012, ApJ, 756,149
  * Bell (2004) Bell, A. R. 2004, MNRAS, 353, 550
  * Berezhko & Völk (2010) Berezhko,E. & Völk 2010, A&A, 511,34
  * Bocchino et al. (2000) Bocchino, F., Vink, J., Favata, F., Maggio, A., & Sciortino, S. 2000, A&A, 360, 671
  * Burlaga & Lazarus (2000) Burlaga, L. F., & Lazarus, A. J. 2000, J. Geophys. Res., 105, 2357
  * Bykov et al. (2008) Bykov, A. M., Uvarov, Y. A., & Ellison, D. C. 2008, ApJ, 689, L133
  * Cassam-Chenaï et al. (2004) Cassam-Chenaï, G., Decourchelle, A., Ballet, J., et al. 2004, A&A, 427, 199
  * Ellison et al. (2010) Ellison, D. C., Patnaude, D. J., Slane, P., & Raymond, J. 2010, ApJ, 712, 287
  * Ellison et al. (2012) Donald C. Ellison, Patrick Slane, Daniel J. Patnaude, & Andrei M. Bykov 2012, ApJ, 744, 39
  * Fan et al. (2010) Fan, Z. H., Liu, S. M., & Fryer, C. L. 2010, MNRAS, 406, 1337-1349
  * Fukui et al. (2003) Fukui, Y., Moriguchi, Y., Tamura, K., et al. 2003, PASJ, 55, L61
  * Giacalone & Jokipii (1999) Giacalone, J., & Jokipii, J. R. 1999, ApJ, 520, 204
  * Giacalone & Jokipii (2007) Giacalone, J., & Jokipii, J. R. 2007, ApJ, 663, L41
  * Guo et al. (2012) Fan Guo, Shengtai Li, Hui Li, Joe Giacalone, J. R. Jokipii, & David Li 2012,ApJ, 747,98
  * Inoue et al. (2009) Inoue, T., Yamazaki, R., & Inutsuka, S. 2009, ApJ, 695, 825
  * Katz & Waxman (2008) Katz, B., & Waxman, E. 2008, JCAP, 01, 018
  * Kesteven & Caswell (1987) Kesteven, M. J., & Caswell, J. L. 1987, A&A, 183, 118
  * Koyama et al. (1995) Koyama, K., et al. 1995, Nature, 378, 225
  * Lemoine-Goumard et al. (2012) Lemoine-Goumard, M., Renaud, M., Vink, J., Allen, G. E., Bamba, A., Giordano, F., & Uchiyama, Y. 2012, A&A, 545, A28
  * Li et al. (2011) Li, H., Liu, S., & Chen, Y. 2011, ApJ, 742,L10
  * Liu et al. (2008) Liu, S., Fan, Z., Fryer, C. L., Wang, J. M., & Li, H. 2008, ApJ, 683, L163
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S. et al. 2007, ApJS, 170, 228
  * Moriguchi et al. (2005) Moriguchi, Y., Tamura, K.,& Tawara, Y., et al. 2005, ApJ, 631, 947
  * Morlino & Caprioli (2012) Morlino, G. & Caprioli, D. 2012, A&A, 538,81
  * Morlino et al. (2009) Morlino, G., Amato, E.& Blasi, P. 2009, MNRAS, 392,240
  * Obergaulinger et al. (2014) Obergaulinger, M., Iyudin, A. F., M ller, E., Smoot, G. F. 2014, MNRAS, 437, 976
  * Pacholczyk (1970) Pacholczyk, A. G. 1970, Radio Astrophysics (San Francisco: Freeman)
  * Pfeffermann & Aschenbach (1996) Pfeffermann, E. & Aschenbach, B. 1996, in Roentgenstrahlung from the Universe, ed. H. U. Zimmermann, J. Tr mper, & H. Yorke, 267
  * Pisarskietal et al. (1984) Pisarski, P. L., Helfand, D. J., & Kahn, S. M. 1984, ApJ, 277, 710
  * Porter et al. (2006) Porter, T. A., Moskalenko, I. V., & Strong, A. W. 2006, ApJ, 648, L29-L32
  * Riquelme & Spitkovsky (2009) Riquelme, M. A., & Spitkovsky, A. 2009, ApJ, 694, 626
  * Rosado, et al. (1996) Rosado, M., Ambrocio-Cruz, P., Le Coarer, E., & Marcelin, M. 1996, A&A, 315, 243
  * Sano et al. (2010) Sano, H., Sato, J., Horachi, H., et al. 2010, ApJ, 724, 59
  * Smith (1997) Smith, R. C. 1997, AJ, 114, 2664
  * Tanaka et al (2008) Tanaka, T., Uchiyama, Y., Aharonian, F. A., et al. 2008, ApJ, 685, 988
  * Tian, et al. (2008) Tian, W. W., Leahy, D. A., Haverkorn, M., & Jiang, B. 2008, ApJ, 679, L85
  * Tsebrenko & Soker (2013) Tsebrenko, D., & Soker, N. 2013, MNRAS, 435, 320
  * Uchiyama,et al. (2007) Uchiyama, Y., Aharonian, F. A., Tanaka, T., Takahashi, T., & Maeda, Y. 2007, Nature, 449, 576
  * Williams et al. (2011) Williams, B. J. et al. 2011, ApJ, 741, 96(15pp)
  * Vink et al. (2006) Vink, J., Bleeker, J., Van Der Heyden, K., Bykov, A., Bamba, A., & Yamazaki, R. 2006, ApJ, 648, L33
  * Yamazaki et al. (2009) Yamazaki, R., Kohri, K.& Katagiri, H. 2009, A&A, 495,9
  * Yang & Liu (2013) Yang,Chuyuan & Liu, S. M. 2013, ApJ, 773,138
  * Yang et al. (2014) Yang,R. Z., Zhang, X., Yuan, Q., & Liu, S. M. 2014, A&A, 567, A23 (7pp)
  * Yuan et al. (2011) Yuan, Q., Liu, S. M., Fan, Z. H., Bi, X. J., & Fryer, C. L. 2011, ApJ, 735, 120(9pp)
  * Yuan et al. (2012) Yuan, Q., Liu, S. M., & Bi, X. J. 2012, ApJ, 761, 133
  * Yuan et al. (2014) Yuan, Q., Huang, X. Y., Liu, S. M., & Zhang, B. 2014, ApJL, 785, L22(5pp)


\ X_{\nu,\mu}-X_{\mu,\nu} = A_{\nu,\mu}-A_{\mu,\nu}\,,\\
\ \mu_0 J^\mu =g^{\nu\sigma}g^{\mu\tau}\left(X_{\tau\sigma,\nu}-X_{\sigma\tau,\nu}+\gamma_\nu X_{\tau,\sigma}-\gamma_\nu X_{\sigma,\tau}\right) \,.
\label{ContactEOM}
\end{dcases}
\end{equation}

We recover only a part of the holonomy, as it is expected in premulticontact Lagrangian systems.
Imposing that ${\bf X}_\L$ is holonomic and denoting $X_{\mu\nu,\rho}-X_{\nu\mu,\rho}=F_{\nu\mu,\rho}$, the last equation of (<ref>) is

$$\mu_0 J^\mu =g^{\nu\sigma}g^{\mu\tau}\left(F_{\sigma\tau,\rho}+\gamma_\nu F_{\sigma\tau}\right)\,.$$

In order to recover an expression involving the electric and magnetic fields, we consider the case $\displaystyle g^{\mu\nu} = \frac{1}{\sqrt{1+\chi_m}}\text{diag}\Big((1+\chi_e)(1+\chi_m),-1,-1,-1 \Big)$, where $\chi_e$ and $\chi_m$ are the electric and magnetic susceptibilities respectively. Then, letting $\displaystyle\gamma^\mu = \left(\frac{\gamma}{c},\pmb{\gamma}\right)$ and $\displaystyle J^\mu = \left(c\rho,\textbf{j}\right)$, the last equation of (<ref>) reads in vector notation as

\begin{align*}
\frac{\rho}{\epsilon_0}&=(1+\chi_e)\big(\nabla\cdot\textbf{E}+\pmb{\gamma}\cdot\textbf{E}\big)\,,
\\
\textbf{j}&=\-(1+\chi_e)\epsilon_0\left(\frac{\partial\textbf{E}}{\partial t}+\gamma\textbf{E}\right)+\frac{1}{1+\chi_m}\left(\nabla\times\frac{\textbf{B}}{\mu_0}+\pmb{\gamma}\times\frac{\textbf{B}}{\mu_0}\right) \,,
\end{align*}

which are the contact Maxwell's equations.
These equations are the so-called “second pair of Maxwell's equations” and correspond to the Gauss law (for electric fields) and the Ampère–Maxwell law for linear materials (or also for vacuum) when $\gamma_\nu = 0$.
The “first pair of Maxwell equations”; i.e., the  Gauss law (for magnetic fields) and the Faraday–Henry–Lenz law
are the same as in the vacuum case,
since these laws are just stating that the curvature of the connection, whose coefficients are the electromagnetic tensor field $F$, is closed.

§ CONCLUSIONS AND OUTLOOK

In this paper we have introduced and studied (pre)multicontact structures
as a generalization of (pre)cocontact and (pre)multisymplectic structures.

We have stated the properties defining (pre)multicontact manifolds in general.
The bundle structure of these kinds of manifolds has been analyzed, also proving a Darboux-type theorem to find charts of adapted coordinates.
This is useful, in first place, to study these manifolds from an abstract point of view and, at the same time, they are the setting for describing non-conservative first-order classical field theories.
To do the latter, we need to restrict the type of (pre)multicontact structures involved.
This is achieved by introducing an additional condition that ensures that the pdes associated with the structure
(which are the field equations in the non-conservative field theories) are of variational type.
We have characterized the field equations using different geometrical tools (sections, multivector fields and Ehresmann connections) that are useful in different contexts.
Then, the Lagrangian and the Hamiltonian formalisms of non-conservative field theories have been developed for the regular and the almost-regular situations.
Finally, the extension of the Herglotz variational Principle has been studied in the Lagrangian, Hamiltonian and general cases, in order to derive the (pre)multicontact field equations variationally.

As illustrative examples, from this framework we have recovered the cocontact formulation describing non-conservative time-dependent mechanics,
we have analyzed the $1$-dimensional wave equation with time-dependent damping, and Maxwell's electromagnetism with damping terms.

The present work opens many new lines of research.
Some of the possible topics are:

* To do an accurate analysis of the case of singular systems and the corresponding Lagrangian and Hamiltonian constraint algorithms,
characterizing geometrically the constraint submanifolds
and studying the equivalence between the solutions to the Lagrangian and the Hamiltonian field equations.

* To develop a Hamilton–Jacobi equation for non-conservative field theories similar to that for contact Hamiltonian systems [26, 27] and multisymplectic field theories [33].

* To study the reduction of multicontact systems in the presence of a group of symmetries, introducing a good definition of momentum map.
For multisymplectic reduction see, for instance, [7, 40] and references therein.

* To study the different types of submanifolds of a multicontact manifold and to find an analogue to the coisotropic reduction.

* To study the formalism in Cauchy data space in the multicontact setting.

* To establish a bracket similar to that of multisymplectic theories, which in this case will not satisfy Leibniz's rule-like properties, as in the case of contact Hamiltonian systems.
Related to this, it would be interesting to study the relationship
between (pre)multicontact structures and Jacobi bundles.

* To provide a geometric framework for new gravitational theories to describe astrophysical observations [51, 82].

* It is expected that the multicontact framework can be used to formulate continuum mechanics equations [76].

§ ACKNOWLEDGMENTS

We acknowledge financial support of the
Ministerio de Ciencia, Innovación y Universidades (Spain), projects PGC2018-098265-B-C33,
PID2019-106715GB-C21, and PID2021-125515NB-C21,
the ICMAT Severo Ochoa project CEX2019-000904-S, and grant EIN2020-112197, funded by AEI/10.13039/501100011033 and European Union NextGenerationEU/PRTR,
and the financial support for research groups AGRUPS-2022 of the Universitat Politècnica de Catalunya (UPC).
X. Rivas acknowledges financial support of the Novee Idee 2B-POB II, project PSP: 501-D111-20-2004310 funded by the “Inicjatywa Doskonałości-Uczelnia Badawcza” (IDUB) program.
We thank the referees for their helpful comments.

[1]
R. Abraham, J. Marsden,
Foundations of mechanics ($2nd$ ed.),
Benjamin–Cummings, Redwood City CA, 1987.

[2]
V. Aldaya, J.A. de Azcárraga,
“Geometric formulation of classical mechanics and field theory”,
Riv. Nuovo Cim. 3(1)
(1980) 1–66. (<https://doi.org/10.1007/BF02906204>).

[3]
V. Apostolov, D. Calderbank, P. Gauduchon, E. Legendre,
“Toric contact geometry in arbitrary codimension”,
Int. Math. Res. Not. IMRN 2020(8) (2018) 2436–2467.

[4]
V.I. Arnold,
Mathematical methods of classical mechanics,
Graduate Texts in Mathematics 60. Springer-Verlag, New York, 1989.

[5]
G. Bande, A. Hadjar,
“Contact pairs”,
Tohoku Math. J. 57(2) (2005) 247–260.

[6]
A. Banyaga, D. F. Houenou,
A brief introduction to symplectic and contact manifolds,
World Scientific, Singapore, 2016.

[7]
C. Blacker, “Reduction of multisymplectic manifolds",
Lett. Math. Phys. 111 (2021) 64. (<https://doi.org/10.1007/s11005-021-01408-y>).

[8]
P. Bolle,
“Une condition de contact pour les sous-variétés coısotropes d’une variété symplectique”,
C. R. Math. Acad. Sci. Sér. 1 1 (1996) 83–86.

[9]
A. Bravetti,
“Contact Hamiltonian dynamics: the concept and its use”,
Entropy 19(10) (2017) 535.

[10]
A. Bravetti, H. Cruz, D. Tapias,
“Contact Hamiltonian mechanics”,
Ann. Phys. (N.Y.) 376 (2017) 17–39.

[11]
A. Bravetti, M. de León, J.C. Marrero, E. Padrón,
“Invariant measures for contact Hamiltonian systems: symplectic sandwiches with contact bread”,
J. Phys. A: Math. Gen. 53(45)
(2020) 455205.

[12]
J.F. Cariñena, M. Crampin,  L.A. Ibort,
“On the multisymplectic formalism for first order field theories”,
Diff. Geom. Appl. 1(4) (1991) 345–374.

[13]
J. Cariñena, P. Guha,
“Nonstandard Hamiltonian structures of the Liénard equation
and contact geometry”,
Int. J. Geom. Meth. Mod. Phys. 16(supp 01) (2019) 1940001.

[14]
F. Cantrijn, L.A. Ibort, M. de León,
“Hamiltonian structures on multisymplectic manifolds”,
Rend. Sem. Mat. Univ. Pol. Torino 54(3) (1996) 225–236.

[15]
F. Cantrijn, L.A. Ibort, M. de León,
“On the geometry of multisymplectic manifolds”,
J. Austral. Math. Soc. (Series A) 66(3) (1999) 303–330.

[16]
E. Cartan,
Leçons sur les invariants integraux, Hermann, Paris, 1922.

[17]
M. Castrillón López, J. Muñoz Masqué, “The geometry of the bundle of connections”, Math. Z. 236(4) (2001) 797–811.

[18]
F. M. Ciaglia, H. Cruz, G. Marmo,
“Contact manifolds and dissipation, classical and quantum”,
Ann. Phys. (N.Y.) 398 (2018) 159–179.

[19]
U.N. de Almeida,
Contact Anosov actions with smooth invariant bundles 59 (2018), Ph D. thesis,
Universidade de São Paulo.

[20]
T. de Donder,
Théorie Invariantive du Calcul des Variations.
Gauthier-Villars, Paris, 1935.

[21]
M. de León, J. Gaset, X. Grácia, M.C. Muñoz-Lecanda, X. Rivas,
“Time-dependent contact mechanics”,
Monatsh. Math. (2022).

[22]
M. de León, J. Gaset, M Laínz-Valcázar,
“Inverse problem and equivalent contact systems”,
J. Geom. Phys. 176 (2022) 104500.

[23]
M. de León, J. Gaset, M. Laínz-Valcázar, M.C. Muñoz-Lecanda, N. Román-Roy,
“Higher-order contact mechanics”,
Ann. Phys. 425 (2021) 168396.

[24]
M. de León, J. Gaset, M. Laínz-Valcázar, X. Rivas, N. Román-Roy,
“Unified Lagrangian–Hamiltonian formalism for contact systems”,
Fortsch. Phys. 68(8) (2020) 2000045.

[25]
M. de León, M. Laínz-Valcázar,
“Singular Lagrangians and precontact Hamiltonian Systems”,
Int. J. Geom. Meth. Mod. Phys.,  16(10) (2019) 1950158.

[26]
M. de León, M. Laínz, A. López-Gordón, X. Rivas,
“Hamilton–Jacobi theory and integrability for autonomous and non-autonomous contact systems”, arXiv: 2208.07436 [math.DG] (2022).

[27]
M. de León, M. Laínz, A. Muñiz-Brea,
“The Hamilton–Jacobi Theory for Contact Hamiltonian Systems”,
Mathematics 9(16) (2021) 1993. (<https://doi.org/10.3390/math9161993>).

[28]
M. de León, B. López, J. C. Marrero, E. Padrón,
“Lichnerowicz–Jacobi cohomology and homology of Jacobi manifolds: modular class and duality”, arXiv: 9910079 [math.DG] (1999).

[29]
M. de León, J. Marín-Solano, J.C. Marrero,
“A Geometrical approach to Classical Field Theories:  A constraint algorithm for singular theories”,
Proc.  New Develops. Dif. Geom.,
L. Tamassi, J. Szenthe eds., Kluwer Acad. Press, (1996) 291–312.

[30]
M. de León, D. Martín de Diego, A. Santamaría-Merino,
“Tulczyjew Triples and Lagrangian Submanifolds in Classical Field Theories”,
in Applied Differential Geometry and Mechanics, 189,
Academia Press: Ghent, Belgium, 2003.

[31]
M. de León, D. Martín de Diego, A. Santamaría-Merino,
“Symmetries in classical field theory”,
Int. J. Geom. Meth. Mod. Phys. 1(5) (2004) 651–710.

[32]
M. de León, P.R. Rodrigues,
Generalized classical mechanics and field theory. A geometrical approach of Lagrangian and Hamiltonian formalisms involving higher order derivatives,
North-Holland Math. Studies 112
Notes Pure Math. 102.
North-Holland Publishing Co., Amsterdam, 1985.

[33]
M. de León, J.C. Marrero, D. Martín de Diego,
“A geometric Hamilton–Jacobi theory for classical field theories”,
in Variations, Geometry and Physics 129–140. Nova Sci.Publ., New York, 2009.

[34]
M. de León, M. Salgado, S. Vilariño,
Methods of Differential Geometry in Classical Field Theories: $k$-Symplectic and $k$-Cosymplectic Approaches,
World Scientific, Hackensack, 2016.

[35]
J. de Lucas, X. Rivas,
“Contact Lie systems: theory and applications”, arXiv: 2207.04038 [math-ph] (2022).

[36]
A. Echeverría-Enríquez, M.C. Muñoz-Lecanda, N. Román-Roy,
“Geometry of Lagrangian first-order classical field theories”,
Fortsch. Phys. 44 (1996) 235–280.

[37]
A. Echeverría-Enríquez, M.C. Muñoz-Lecanda, N. Román-Roy,
“Multivector fields and connections: Setting Lagrangian equations in field theories”,
J. Math. Phys. 39(9) (1998) 4578–4603.

[38]
A. Echeverría-Enríquez, M.C. Muñoz-Lecanda, N.
“Multivector Field Formulation of Hamiltonian
Field Theories: Equations and Symmetries”,
J. Phys. A: Math. Gen. 32(48) (1999) 8461–8484.

[39]
A. Echeverría-Enríquez, M.C. Muñoz-Lecanda, N. Román-Roy,
“Geometry of multisymplectic Hamiltonian first-order field theories”,
J. Math. Phys. 41(11)  (2000) 7402–7444.

[40]
A. Echeverría-Enríquez, M.C. Muñoz-Lecanda, N. Román-Roy,
“Remarks on multisymplectic reduction”,
Rep. Math. Phys. 81(3)  (2018) 415–424.

[41]
C. Ehresmann,
“Espaces fibrés associéés”, C. R. Acad. Sci. Paris 213 (1941) 762–764.

[42]
C. Ehresmann,
“Sur les espaces fibrés associés à une variété difféérentiable”, C. R. Acad. Sci. Paris 216 (1943) 628–630.

[43]
C. Ehresmann, J. Feldbau,
“Sur les propriétés d'homotopie des espaces fibrés”,
C. R. Acad. Sci. Paris 212 (1941) 945–948.

[44]
D. Finamore,
“Contact foliations and generalised Weinstein conjectures”,
arXiv:2202.07622 [math.SG] (2022).

[45]
P.L. García,
“The Poincaré–Cartan invariant in the calculus of variations”,
Symp. Math. 14 (1973) 219–246.

[46]
J. Gaset, X. Gràcia, M. Muñoz-Lecanda, X. Rivas,
N. Román-Roy,
“A contact geometry framework for field theories with dissipation”,
Ann. Phys. 414 (2020) 168092.

[47]
J. Gaset, X. Gràcia, M. Muñoz-Lecanda, X. Rivas,  N. Román-Roy,
“New contributions to the Hamiltonian and Lagrangian contact formalisms for dissipative mechanical systems and their symmetries”,
Int. J. Geom. Meth. Mod. Phys. 17(6) (2020) 2050090.

[48]
J. Gaset, X. Gràcia, M. Muñoz-Lecanda, X. Rivas, N. Román-Roy,
“A $k$-contact Lagrangian formulation for nonconservative field theories”,
Rep. Math. Phys. 87(3) (2021) 347–368.

[49]
J. Gaset, M. Laínz, A. Mas, X. Rivas,
“The Herglotz principle for field theories”,
arXiv:2211.17058 [math-ph] (2022).

[50]
J. Gaset, A. Marín-Salvador,
“Application of Herglotz's variational principle to electromagnetic systems with dissipation”,
Int. J. Geom. Meth. Mod. Phys. 19(supp 10) (2022) 2250156. (<https://doi.org/10.1142/S0219887822501560>).

[51]
J. Gaset, A. Mas,
“A variational derivation of the field equations of an action-dependent Einstein-Hilbert Lagrangian”, arXiv:2206.13227 [gr-qc] (2022).

[52]
F. Gay-Balmaz, J.C. Marrero, N. Martínez,
“A new canonnical affine bracket formulationn of Hamiltonian classical field theories of first-order”,
arXiv:2209.08736 [math-ph] (2022).

[53]
H. Geiges,
An Introduction to Contact Topology,
Cambridge University Press, Cambridge, 2008.

[54]
B. Georgieva, R. Gunther, T. Bodurov,
“Generalized variational principle of Herglotz for several independent variables. First Noether-type theorem”,
J. Math. Phys. 44, 3911 (2003).
(<https://doi: 10.1063/1.1597419>).

[55]
G. Giachetta, L. Mangiarotti, G. Sardanashvily,
New Lagrangian and Hamiltonian methods in field theory,
World Scientific Publishing Co., Inc., River Edge, NJ, 1997.

[56]
C. Godbillon,
Géométrie différentielle et mécanique analytique,
Hermann, Paris, 1969.

[57]
H. Goldschmidt, S. Sternberg,
“The Hamilton–Cartan formalism in the calculus of variations”,
Ann. Inst. Fourier Grenoble 23(1) (1973) 203–267.

[58]
M.J. Gotay,
“A multisymplectic framework for classical field theory and the calculus of variations. I. Covariant Hamiltonian formalism”,
in Mechanics, analysis and geometry: 200 years after Lagrange, 203–235,
North-Holland Delta Ser., North-Holland, Amsterdam, 1991.

[59]
K. Grabowska, J. Grabowski,
“A novel approach to contact Hamiltonians and contact Hamilton–Jacobi theory”,
arXiv:2207.04484 [math.SG] (2022).

[60]
K. Grabowska, J. Grabowski,
“Contact geometric mechanics: the Tulczyjew triples”,
arXiv:2209.03154 [math.SG] (2022).

[61]
X. Gràcia, X. Rivas, N. Román-Roy.
“Skinner–Rusk formalism for $k$-contact systems”,
J. Geom. Phys. 172 (2022) 104429.

[62]
F. Hélein, J. Kouneiher,
“Covariant Hamiltonian formalism for the calculus of variations with several variables: Lepage–Dedecker versus De Donder–Weyl”,
Adv. Theor. Math. Phys. 8 (2004) 565–601.

[63]
G. Herglotz, “Berührungstransformationen”, Lectures at the University of Gottingen, 1930.

[64]
G. Herglotz,
Vorlesungen über die Mechanik der Kontinua.
Teubner-Archiv zur Mathematik 3;
Teubner, Leipzig, 1985.

[65]
L.A. Ibort, A. Echeverría-Enríquez, M.C Muñoz-Lecanda, N. Román-Roy,
“Invariant forms and automorphisms of locally homogeneous multisymplectic manifolds”,
J. Geom. Mech. 4(4) (2012) 397–419.

[66]
I.V. Kanatchikov,
“Canonical Structure of Classical Field Theory in the
Polymomentum Phase Space”,
Rep. Math. Phys. 4(1) (1998) 49–90.

[67]
J. Kijowski,
“A finite-dimensional canonical formalism in the classical field theory”,
Comm. Math. Phys. 30 (1973) 99–128.

[68]
J. Kijowski, W.M. Tulczyjew,
Symplectic Framework for Field Theories,
Lect. Notes Phys. 170, Springer-Verlag, Berlin-New York, 1979.

[69]
A. L. Kholodenko,
Applications of Contact Geometry and Topology in Physics,
World Scientific, Singapore, 2013.

[70]
O. Krupkova,
“Hamiltonian field theory”,
J. Geom. Phys. 43(2–3) (2002) 93–132

[71]
M. Laínz-Valcázar, M. de León,
“Contact Hamiltonian systems”,
J. Math. Phys. 60(10) (2019) 102902.

[72]
M.J. Lazo, J. Paiva, J.T.S. Amaral, G.S.F. Frederico,
“An action principle for action-dependent Lagrangians: Toward an action principle to non-conservative systems”,
J. Math. Phys. 59(3) (2018) 032902.

[73]
T.H.J. Lepage,
“Champs stationnaires, champs géodésiques et formes intégrables”,
Acad. Roy. Belgique. Bull. Cl. Sci.  (5) 28 (1942) 73–92, 247–265.

[74]
P. Libermann, C. Marle.
Symplectic Geometry and Analytical Mechanics.
D. Reidel, Dordrecht, 1987.

[75]
Q. Liu, P.J. Torres, C. Wang,
“Contact Hamiltonian dynamics: variational principles, invariants,
completeness and periodic behavior”,
Ann. Phys. 395 (2018) 26–44.

[76]
J.E. Marsden, T.J.R. Hughes,
Mathematical foundations of elasticity,
Dover Pub. Inc., New York, 1994.

[77]
J.E. Marsden, S. Shkoller,
“Multisymplectic geometry, covariant Hamiltonians and nonlinear pdes”,
Comm. Math. Phys. 199(2) (1998) 351–395.

[78]
G. Martin,
“A Darboux theorem for multisymplectic manifolds”,
Lett. Math. Phys. 16(2) (1988) 133–138.

[79]
B. Montano,
“Integral submanifolds of $r$-contact manifolds”,
Demonstr. Math. 41(1) (2008) 189–202.

[80]
R. Narasimhan,
Analysis on Real and Complex Manifolds ($3rd$ ed.), Elsevier, North Holland, 1985.

[81]
C. Paufler, H. Römer,
“Geometry of Hamiltonian $n$-vector fields in multisymplectic field theory”,
J. Geom. Phys 44(1) (2002)

[82]
J. Paiva, M.J. Lazo, V.T. Zanchin,
“Generalized nonconservative gravitational field equations from Herglotz action principle”,
Phys. Rev. D 105(12) (2022)

[83]
X. Rivas,
“Nonautonomous $k$-contact field theories”, arXiv: 2210.09166 [math-ph] (2022).

[84]
X. Rivas, D. Torres,
“Lagrangian–Hamiltonian formalism for cocontact systems”, J. Geom. Mech. 15(1) (2022) 1–26. (<https://doi.org/10.3934/jgm.2023001>).

[85]
N. Román-Roy,
“Multisymplectic Lagrangian and Hamiltonian formalisms of classical field theories”,
Symm. Integ. Geom. Meth. Appl. (SIGMA)
5 (2009) 100.

[86]
D.J. Saunders, The geometry of jet bundles, London Math.
Soc., Lect. Notes Ser. 142, Cambridge Univ. Press,
Cambridge, New York, 1989.

[87]
A. Tomassini, L. Vezzoni,
“Contact Calabi–Yau manifolds and special Legendrian submanifolds”,
Osaka J. Math. 45(1) (2008) 127–147.

[88]
L. Vitagliano,
“$L_\infty$-algebras from multicontact geometry”,
Diff. Geom. Appl. 59 (2015) 147–165.

[89]
H. Weyl,
“Geodesic fields in the calculus of variation for multiple integrals”,
Ann. Math. 36(3) (1935) 607–629.