# SOUP-GAN: Super-Resolution MRI Using Generative Adversarial Networks

Kuan Zhang  Haoji Hu  Kenneth Philbrick  Gian Marco Conte  Joseph D. Sobek
Pouria Rouzrokh  Bradley J. Erickson Corresponding author: Bradley J. Erickson
(e-mail: bje@mayo.edu).K. Zhang, G. M. Conte, J. D. Sobek, P. Rouzrokh and B.
J. Erickson are with the Department of Radiology, Mayo Clinic, Rochester, MN.
H. Hu is with the Department of Computer Science & Engineering, the University
of Minnesota, Minneapolis, MN.K. Philbrick is with Google LLC, Mountain View,
CA. Work was conducted while employed at Mayo Clinic.

###### Abstract

There is a growing demand for high-resolution (HR) medical images in both the
clinical and research applications. Image quality is inevitably traded off
with the acquisition time for better patient comfort, lower examination costs,
dose, and fewer motion-induced artifacts. For many image-based tasks,
increasing the apparent resolution in the perpendicular plane to produce
multi-planar reformats or 3D images is commonly used. Single image super-
resolution (SR) is a promising technique to provide HR images based on
unsupervised learning to increase resolution of a 2D image, but there are few
reports on 3D SR. Further, perceptual loss is proposed in the literature to
better capture the textual details and edges than using pixel-wise loss
functions, by comparing the semantic distances in the high-dimensional feature
space of a pre-trained 2D network (e.g., VGG). However, it is not clear how
one should generalize it to 3D medical images, and the attendant implications
are still unclear. In this paper, we propose a framework called SOUP-GAN:
$\mathbb{S}$uper-resolution $\mathbb{O}$ptimized $\mathbb{U}$sing
$\mathbb{P}$erceptual-tuned $\mathbb{G}$enerative $\mathbb{A}$dversarial
$\mathbb{N}$etwork (GAN), in order to produce thinner slice (e.g., high
resolution in the ‘Z’ plane) medical images with anti-aliasing and deblurring.
The proposed method outperforms other conventional resolution-enhancement
methods and previous SR work on medical images upon both qualitative and
quantitative comparisons. Specifically, we examine the model in terms of its
generalization for various SR ratios and imaging modalities. By addressing
those limitations, our model shows promise as a novel 3D SR interpolation
technique, providing potential applications in both clinical and research
settings.

super-resolution, magnetic resonance imaging (MRI), generative adversarial
networks (GAN), 3D perceptual loss, medical imaging interpolation, de-
aliasing, deep learning

##  1 Introduction

Medical imaging modalities (e.g., MRI, CT and ultrasound imaging) are widely
used for clinical diagnosis and research [1]. They reveal internal anatomy and
can provide quantitative measurements of the human body. Acquisition of a
medical image depends on both the technical quality of the instrument and the
conditions of the scan (e.g. patient properties, nature of the image being
acquired), and in most cases, involves compromise, i.e., choosing a shorter
acquisition time in MRI, reducing the radiation dose in CT, and using lower
ultrasound power levels in ultrasound imaging. Thus, the acquired images may
have limited spatial resolution, most commonly lower in the ‘z-plane’. For
instance, MRI using 2D spin-echo technique can provide good in-plane
resolution, but typically have a high thickness for a clinically reasonable
acquisition time, constraining the possibility of capturing important signals
of lesions between the slices. Moreover, the thick-slice images produce poor
quality oblique or perpendicular reconstructions. There is active demand in
practice to generate HR images in order to relieve the burden of patients, to
reduce examination costs and to decrease motion-induced artifacts.

One of the main approaches to accelerate MRI acquisition is through under-
sampling the data in k-space, which inevitably results in aliasing artifacts
by violating the Nyquist-Shannon sampling criteria [2], and then performing
reconstruction of the under-sampled MRI, by solving an ill-posed inverse
problem from an optimization point of view [3]. Currently, parallel imaging
(PI) and compressed sensing (CS) [4] are widely used in clinical practice.
While the former method adopts multiple independent channels to simultaneously
receive each view of the tissue area and utilizes software tools such as SENSE
[5] or GRAPPA [6] to combine the multiple signals, the latter provides a
framework for imaging reconstruction by adding a regularization term subject
to prior sparsity constraints into the optimization problem. Beyond the
conventional techniques, deep learning methods for MRI reconstruction have
been extensively studied to learn the relationship between high-quality ground
truths as target and the corresponding under-sampled images as input [7, 8, 9,
10]. More advanced networks and algorithms are proposed following the idea of
incorporating more “physics” from conventional methods (CS and PI), e.g., the
development of variational networks [11, 12, 13].

Single-image SR focuses on the 2D spatial domain, and refers to the process of
generating HR images from low-resolution (LR) images. Deep learning based SR
models have been actively explored and achieve state-of-the-art performance on
natural image datasets [14, 15, 16, 17, 18]. One important breakthrough is the
perceptual loss function [19, 20, 21], which is proposed to tackle the issue
of blurring generated by optimization of pixel-wise loss functions such as
MSE. Together with GAN loss, perceptual loss has the potential to create
images that are visually indistinguishable from the ground truth images.

Basic SR techniques have been introduced to MRI. Chaudhari et al. [22]
implements convolutional neural networks (CNNs) to learn residual information
between HR thin-slice images and LR thick-slice images using MSE loss. A self-
supervised anti-aliasing algorithm is developed by [23, 24] to resolve the LR
out-of-plane slices by learning from the HR in-plane slices, thus requiring no
external training dataset. High-quality results are reported by the SR MRI
approaches in literature [25, 26, 27, 28, 29]. However, there are certain
limitations that preclude their application in real-world scenarios. Given the
blurred textures trained by MSE loss, perceptual loss with GAN seems necessary
in the single-image SR model, but it is unclear how to generalize the pre-
trained 2D model based perceptual loss to 3D medical images. In modern medical
image deep learning tasks, data augmentation is often applied during data
preprocessing when the acquired data is insufficient. The adaptation of the SR
MRI model to address a range of SR scenarios such as various sampling ratios,
acquisition processes, or medical modalities is still limited. As a result,
interpolation is still used instead of SR for routine data augmentation.

In this study, we proposed a novel algorithm termed SOUP-GAN:
$\mathbb{S}$uper-resolution $\mathbb{O}$ptimized $\mathbb{U}$sing a
$\mathbb{P}$erceptual-tuned $\mathbb{GAN}$ network for generating HR thin-
slice images. Our main contributions are:

  * •

We generalize the application of perceptual loss previously defined on 2D pre-
trained VGG onto 3D medical images. Together with GAN, the new scheme with 3D
perceptual loss significantly improves the perceived image quality over the
MSE trained results.

  * •

We develop a multi-scale SR architecture that works for various sampling
factors with an overall satisfactory performance.

  * •

Two criteria are proposed accounting for different acquisition protocols.
Without examining the acquisition process carefully, the extracted datasets
will not provide accurate mapping from LR slices into HR slices, and therefore
may not work well for all the SR tasks.

  * •

We evaluate the feasibility of designing one model that works for all the
medical images. The proposed model is applied to datasets from other imaging
modalities, e.g., T2-weighted MRI and CT. Combining the above aspects, our
model shows potential to be applied for the novel SR interpolations on medical
images.

##  2 Method

###  2.1 Residual CNNs for 3D SR-MRI

Generally, image reconstruction is the transformation of signals acquired to
images that are meaningful to humans and can be formulated as the multi-
dimensional linear system:

$Ex=y+\varepsilon,$ (1)

which tries to generate the desired images $x={x_{i},i\leq N}$ based on the
observed signals $y={y_{i},i\leq n}$ with $N$ and $n$ as the numbers of
desired and observed slices (images in the Z-dimension), respectively. $E$ is
the encoder operator, representing various image generation or transformation
tasks, e.g., an identity operator with and without convolutions respectively
for image deblurring and denoising, image-space uniform under-sampling for SR,
local masking operators for inpainting and k-space under-samples for MRI
reconstruction. $\varepsilon$ is the coherent noise by the measurement.

SR tasks on 3D MRI specifically can be viewed as an under-sampling encoder
trained to generate HR slices from LR slices. While Eq. 1 itself is ill-posed,
it can be turned into an optimization problem:

$\hat{x}=\arg\min\left(\frac{1}{2}||Ex-y||^{2}_{2}\right),$ (2)

and solved by supervised learning. Given a set of HR ground-truth examples x,
the corresponding LR images y can be obtained through the encoder operator
$E=(x\ast\kappa)\downarrow_{s}$ that typically applies the anti-aliasing
filtering with a subsequent down-sampling by the sampling factor $s$ on the
ground truth slices. The training dataset is composed of the ground-truths $x$
and their associated under-sampled $y$, on which a generator network
$G_{\theta}$ is proposed and trained to approximate the function $F$, the
inverse form of the encoder $E$, such that $x=F(y)$.

The designed generator $G_{\theta}$ consists of a pre-upsampling layer by
cubic interpolation together with a residual network of 3D convolutional
layers, or 3D U-net, on which the data patches can be efficiently trained. For
the task involving a large training dataset, the so-called RRDB (Residual in
Residual Dense Block) introduced in ESRGAN [18] is adopted to better resolve
the residual spatial details, and to alleviate the burden of occupied memory
through such a post-upsampling scheme.

As an outline of the proposed framework SOUP-GAN, Algorithm 1 is provided with
pseudocode to summarize the training process including multiple processing
steps, which are depicted by Fig. 1 and explained in the later subsections.

Algorithm 1 Training the proposed framework SOUP-GAN.

0: HR 3D MRI volume.

0: 1) Data preprocessing:

Create LR volume by data-preprocessing following either the thin-to-thick or
thin-to-thin criteria.

2) Prepare training dataset:

Partition the LR data as input and HR data as ground truth into pairs of
$32\times 32\times 32$ patches.

3) Multi-scale SR network:

Input LR patches to the multi-scale SR network with an appropriate module
entrance and exit according to the associated sampling factor $s$.

4) 3D perceptual loss with GAN:

Based on the pre-trained MSE results, further tune the model by employing the
3D perceptual loss with GAN.

Figure 1: Architecture of our multi-scale SR network with GAN. 3D perceptual
loss is calculated by inputting all the axial, sagittal and coronal surfaces
of the 3D data into the pre-trained VGG, and comparing the differences between
predictions and ground-truths in the high-dimensional feature space.

###  2.2 Data Preprocessing

Data acquired for this study are 3D MRI datasets. Nine (9) T1-weighted brain
MRI volumes with 1 mm isotropic pixel spacing in both the in-plane and out-of-
plane directions were cropped and stacked into one $256\times 256\times 1170$
volume. Those thin slices were HR ground-truth images. The LR volume
($256\times 256\times 234$) was created by the data-preprocessing methods
explained later. Both the input and ground-truth datasets were then
partitioned into pairs of patches. While larger volumetric patches potentially
capture better semantic features [30], we selected a patch size of $32\times
32\times 32$ to reduce computational demands. All patch pairs were shuffled,
and distributed into training, validation and test datasets at a ratio of
$8:1:1$.

The specific way to generate the input datasets really depends on the
acquisition protocols and the up-sampling factor $s$ selected for the SR task.
We define two different criteria accounting for the situations broadly
categorized in the acquisition process, because without taking care of these
variations, the extracted datasets will not provide accurate mapping from LR
slices into HR slices, and therefore will fail to work well for various
resolution-enhanced tasks in practice. In detail, the acquisition process is
categorized into two types: 1) thin-to-thick; and 2) thin-to-thin.

The first type assumes the same FOV of the 3D data (i.e., same Z coverage),
and is subject to the selection of the slice thickness, which determines how
much tissue is to be covered within each slice. Thicker slices can achieve
more efficient scans without slice gaps or spatial aliasing, but may cause
image blurring as partial volume effects. Correspondingly, the thick-slice
representation of the HR ground-truth data is obtained by the finite impulse
response (FIR), low-pass filtering with the filter coefficients as $b=1/s$ to
mimic the averaging and anti-aliasing characteristic with the high thickness,
followed by a subsequent down-sampling with the factor $s$.

The second type is to expedite the scans by discretizing FOV in Z (i.e., same
2D X-Y resolution). This is useful for the acquisition protocols focused on
the in-plane resolution without considering the perpendicular (Z) planes. Thin
slices with large separation are acquired efficiently. As a result, they are
insufficient to cover the whole FOV in Z direction, therefore inevitably
introducing spatial aliasing in the perpendicular planes [31]. One clinical
example is for diagnosing interstitial lung disease (ILD) in which HR CT data
($\sim$ 1mm thin slice) is acquired at 10 mm intervals, since the disease is
fairly diffuse. Correspondingly to this type of clinical acquisition as well
as imaging interpolation in deep learning tasks, the LR representation of the
HR ground-truth data is obtained directly by the down-sampling with the factor
s. For other cases between the first and second types of data preprocessing
discussed here, a Gaussian filter with a tuned value of the standard deviation
can be adopted. Finally, the input dataset to the network is obtained by a
cubic spline interpolation with s on the LR representations that are obtained
from the algorithms based on the thin-to-thick or thin-to-thin criteria.

###  2.3 Multi-scale Model for SR Interpolation

Figure 2: Quantitative comparisons of measures: (a) RMSE, (b) PSNR and (c)
SSIM in terms of the sampling factor s between the single-scale model, the
multi-scale model, perceptual-tuned GAN (SOUP) and cubic interpolation.
Statistical significance is indicated on the x-axes with * for $p<0.05$ and **
for $p<0.001$.

Supervised medical deep-learning tasks (e.g., classification, detection and
segmentation) are based on the acquired images and associated labels. The
accuracy of the model inference is not only decided by the design of the
networks, but more importantly relies on the data size and quality. Data
augmentation is often performed when the training data is insufficient,
imbalanced or incompatible in shape. Conventional interpolations (such as
linear, square, and cubic) are widely used and can provide approximate
intermediate slices. However, when the sampling factor $s$ increases, the
interpolated new images become more blurred and further away from the optimal
resolution. Whether or not healthier and more realistic images could be
generated seems to become the bottle-neck to further improve the model
performance through data augmentation. From the observation on SR tasks, we
design a multi-scale model to provide image interpolation with superior
resolution, upon a user-selected sampling factor $s$.

In our multi-scale architecture, the single-scale model described earlier
provides the starting point for learning the residual between LR and HR
slices. Pre- and post- processing modules are added respectively at the top
and the bottom layers, to deal with the differential information between
scales. As a result, we construct one model that is able to provide
predictions for multiple cases with different sampling factors
$s\in\\{2,3,4,5,6\\}$, see Fig. 2, while most of the model parameters are
shared by the backbone across the scales.

For fractional values of the sampling factor $s$ between the selected
integers, we incorporate them by interpolating the corresponding parameters of
the two neighboring networks to derive an interpolated model:

$\theta_{G}^{s}=(1-\alpha)\theta_{G}^{m}+\alpha\theta_{G}^{m+1},$ (3)

where $s$ is between the integer $m$ and $m+1$ with the ratio $\alpha=s-m$. In
conclusion, the proposed multi-scale SR model together with the parameter
interpolation scheme provides an efficient tool to perform the image
interpolation with higher resolutions than conventional interpolation methods.

###  2.4 3D Perceptual Loss

In previous single image SR tasks studied on natural images, it was observed
that training based on a pixel-wise loss, e.g., L1, MSE, or PSNR can achieve
relatively high metric scores, but often results in overly-smooth SR images
with blurring, since the perceptual quality of textures, edges and high-
frequency details are ignored by the pixel-wise losses. Based on the idea of
being closer to semantic similarity, Gatys et al. [19] proposes perceptual
loss, which is defined by evaluating the Euclidean distances between the
feature maps extracted from a pre-trained VGG network instead of low-level
pixel-wise error. For instance, perceptual loss based on the high-level
representations at the $l$-th layer of the VGG19 network [32] is formulated
as:

$L_{per}=\frac{1}{hwc}\sum_{i,j,k}\left(\phi^{l}_{i,j,k}(x)-\phi^{l}_{i,j,k}(F(y))\right)^{2},$
(4)

where $h$, $w$ and $c$ are the height, width and channel number of the
representation on the layer $l$, respectively.

It is worth noting that SR MRI tasks deal with 3D medical images, while the
VGG network takes 2D images as input. To generalize the single image SR with
perceptual loss into 3D, we propose to adopt the original slices together with
adjacent planes of the ground truth volume $x$ and the predicted volume
$F(y)$. That is, all the axial, sagittal and coronal images of the 3D data
will be imported into the pre-trained VGG and contribute to the perceptual
loss calculation (see Fig. 1). In this way, the 3D MRI volume is perceptually
resolved in terms of projected slices along the three principal directions. We
also note that we used the deep intermediate layer VGG19-54 before activation
rather than after activation.

Figure 3: Test examples on brain and knee MRI showing perpendicular
reconstructions from original thick slices, using tricubic interpolation
(TCI), MSE and GAN using perceptual loss (SOUP). We note [22] proposed a
single-scale residual-based 3D convolutional neural network with MSE, similar
to the third column here. Our SOUP approach better resolves the perceptual
details and is more generally applicable to various sampling factors and
different imaging modalities.  Figure 4: Comparison between single-scale model
and multi-scale model on T1-weighted brain MRI. The figure also includes
original slices, cubic interpolation and ground-truth slices as a reference.
Residual images are added to better visualize the difference of the generated
images to the ground-truths.

###  2.5 GAN

In addition to the perceptual loss, we define a discriminator network
$D_{\theta}$, together with the generator model $G_{\theta}$, based on the
Generative Adversarial Networks structure. The generator $G_{\theta}$ and
discriminator $D_{\theta}$ are optimized in an alternating manner to solve the
adversarial min-max problem:

$\min_{G}\max_{D}\mathbb{E}_{x\sim p_{\rm{gt}}}[\log D(x)]+\mathbb{E}_{y\sim
p_{\rm{train}}}[\log(1-D(G(y)))],$ (5)

which aims to train a generative model $G_{\theta}$ trying to fool the
discriminator $D_{\theta}$, while the discriminator is trained simultaneously
to distinguish the generated SR images from ground truth images. With this
approach of adversarial training between $G_{\theta}$ and $D_{\theta}$, the
generated images are encouraged residing in the featured manifold of real
images, potentially with perceptually superior quality that is visually
indistinguishable from the ground truth images.

##  3 Experiments

###  3.1 Training details

The training process is divided into two stages. The generator model
$G_{\theta}$ is first trained by the MSE loss. The Adam optimization is used.
The learning rate is initialized as $3\times 10^{-4}$, and decayed by a factor
of 3 upon three cycles, if the validation loss is no longer updated for the
current cycle. Then, the model is further tuned by employing the perceptual
loss and GAN. The total loss for $G_{\theta}$ upon the refined tuning is:

$L_{G}=L_{per}+\lambda L_{GAN}+\mu L_{MSE},$ (6)

where the weights $\lambda$ and $\mu$ are selected as 0.01 and 0.001. The pre-
trained pixel-wise model is used to initialize the second stage of training.
Transfer learning helps the perceptual-tuned GAN model converge more
efficiently and obtain more visually pleasing results. The generator
$G_{\theta}$ and the discriminator $D_{\theta}$ are alternately updated until
the model converges to produce visually satisfactory images and no further
improvement is visually apparent.

###  3.2 3D perceptual loss performance

As discussed in section 2, pixel-wise loss usually achieves a relatively high
metric score but fails in resolving the textures, edges, and high-frequency
details, resulting in blurring of the images that are produced. The scheme
that we propose generalizes the single-image perceptual loss based on pre-
trained 2D VGG network into 3D, which is applicable for sliced medical
datasets.

As the first example of our qualitative results, Fig. 3 shows the difference
between using MSE loss and 3D perceptual loss with GAN, to train the network
respectively. The two test datasets included in Fig. 3 are the T1-weighted
brain MRI (reformatted coronal images in the first row, and reformatted
sagittal images in the second), and knee MRI in the third row. The images
appear better compared to the original thick slices, the cubic-spline
interpolated images, and the MSE results. They are further improved by 3D
perceptual loss with GAN, achieving the highest perceptual quality.

###  3.3 Single-scale and multi-scale model comparison

Fig. 3 reveals the possibility of using SR technique as a better
“interpolation”, given the active demand for data augmentation in medical-
imaging deep-learning tasks. The feasibility of this idea requires an adaptive
SR model, which is designed to incorporate random sampling factors. Fig. 4
compares the single-scale SR model with our multi-scale SR model at the
sampling factor equal to 5. We can imagine that, to some extent, the multi-
scale model has to sacrifice partial resolutions, as a trade-off for allowing
multiple scale entries. But still, the images obtained by the multi-scale
model possess a satisfactory quality compared to the TCI and single-scale
model.

Figure 5: Evaluating on T1-weighted brain MRI test dataset with quantitative
comparisons of measures: (a) RMSE, (b) PSNR and (c) SSIM in terms of the
sampling factor s. Statistical significance is indicated on the x-axes with *
for $p<0.05$ and ** for $p<0.001$.  Figure 6: Comparison between multi-scale
model and cubic interpolation at different scales.  Figure 7: Comparison
between single-scale model and multi-scale model performance on T2-FLAIR brain
MRI.  Figure 8: Evaluation using T2-weighted brain MRI images with
quantitative comparisons of measures: (a) RMSE, (b) PSNR and (c) SSIM in terms
of the sampling factor s. Statistical significance is indicated on the x-axes
with * for $p<0.05$ and ** for $p<0.001$.  Figure 9: Comparison between multi-
scale model and cubic interpolation on T2-FLAIR brain MRI at different scales.
Figure 10: Example of SR interpolation on CT. Edges and textual details are
better resolved by SR. Blurring is reduced.

Beyond the single-scale limit, our multi-scale model is demonstrated to work
for user-selected scale values, which is the key to turning the SR model into
a useful interpolation tool. Fig. 5 reports the quantitative data performed on
T1-weighted brain MRI test dataset with the metrics at different scales,
including the root of mean squared error (RMSE), peak signal to noise ratio
(PSNR) and the structural similarity (SSIM). Fig. 6 shows the generated images
at different scales, obtained by SOUP (the multi-scale model) and compared to
the conventional interpolation results. As a sum-up of the comparison, while
the multi-scale SR model and the conventional interpolation method converge to
a similar quality at the scale $s\leq 2$, the SR model starts to outperform
the interpolation and achieves better and better image qualities, when the
scale $s$ increases.

###  3.4 Application to other contrast types of MRI images

In the literature, single image SR models are trained and tested on datasets
of natural images, collected from different categories with various
distributions. Similarly, an open question for medical images is whether it is
possible to develop a generalized SR model that works for various medical
image data types (e.g., modalities or contrast types) with an arbitrary
sampling factor $s$ selected by the user. To address this question, we apply
our multi-scale SR model to the test datasets of other medical imaging types,
e.g., T2-FLAIR MRI and CT. Fig. 7 shows the comparison between single-scale
model and multi-scale model on FLAIR brain MRI, Fig. 8 quantitatively
evaluates the performances of the models on the FLAIR MRI dataset with
different metrics, and Fig. 9 compares SOUP with conventionally interpolated
images at different sampling factor $s$.

Figure 11: Another example of SR interpolation on CT data.

###  3.5 Generalization to other medical imaging modalities, e.g., CT

We also apply our multi-scale SR model on CT data. Fig. 10 and 11 show the SR-
interpolated (SOUP) and conventionally interpolated images on abdominal and
pelvic CTs. It is clear from both examples that the SR-interpolated results
achieve significantly higher quality than the conventionally interpolated ones
in regard to the edges, textual details and blurring. Those comparisons
demonstrate our multi-scale SR model (SOUP) as a widely generalizable tool for
different sampling factors and various applied imaging modalities, which paves
the way for more advanced medical image interpolation through deep-learning
SR.

##  4 Conclusions

To our knowledge, this is the first published work to combine the deep-
learning SR method with a 3D perceptual-tuned GAN network, and to examine its
application and generalization in terms of different sampling scales and
various medical imaging modalities, as a novel 3D image interpolation
technique. The superior-resolution thin slices obtained by the model
outperform other conventional resolution-enhancement methods and previous SR
work on medical images upon both qualitative and quantitative comparisons.
Based on the promising simulation data, the proposed work shows potential in
both clinical tasks, such as reducing the acquisition time and further
resolving the MRI scans in retrospective studies, and research tasks, e.g.,
serving as a novel SR interpolation method for data augmentations with further
applications in lesion measurements, radiomics and automatic segmentations.

## References

  * [1] M. Mahesh, “The essential physics of medical imaging, third edition.” _Medical Physics_ , vol. 40, no. 7, p. 077301, 2013.
  * [2] H. Nyquist, “Certain topics in telegraph transmission theory,” _Transactions of the American Institute of Electrical Engineers_ , vol. 47, no. 2, pp. 617–644, 1928.
  * [3] M. Benning and M. Burger, “Modern regularization methods for inverse problems,” _Acta Numerica_ , vol. 27, pp. 1–111, 2018.
  * [4] D. L. Donoho, “Compressed sensing,” _IEEE Transactions on Information Theory_ , vol. 52, no. 4, pp. 1289–1306, 2006.
  * [5] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, and P. Boesiger, “Sense: sensitivity encoding for fast mri.” _Magn Reson Med_ , vol. 42, no. 5, pp. 952–962, Nov 1999.
  * [6] M. A. Griswold, P. M. Jakob, R. M. Heidemann, M. Nittka, V. Jellus, J. Wang, B. Kiefer, and A. Haase, “Generalized autocalibrating partially parallel acquisitions (grappa),” _Magnetic Resonance in Medicine_ , vol. 47, no. 6, pp. 1202–1210, 2002.
  * [7] G. Yang, S. Yu, H. Dong, G. Slabaugh, P. L. Dragotti, X. Ye, F. Liu, S. Arridge, J. Keegan, Y. Guo, and D. Firmin, “Dagan: Deep de-aliasing generative adversarial networks for fast compressed sensing mri reconstruction,” _IEEE Transactions on Medical Imaging_ , vol. 37, no. 6, pp. 1310–1321, 2018.
  * [8] F. Liu, A. Samsonov, L. Chen, R. Kijowski, and L. Feng, “Santis: Sampling-augmented neural network with incoherent structure for mr image reconstruction,” _Magnetic Resonance in Medicine_ , vol. 82, no. 5, pp. 1890–1904, 2019.
  * [9] H. K. Aggarwal, M. P. Mani, and M. Jacob, “Modl: Model-based deep learning architecture for inverse problems,” _IEEE Transactions on Medical Imaging_ , vol. 38, no. 2, pp. 394–405, 2019.
  * [10] T. M. Quan, T. Nguyen-Duc, and W. K. Jeong, “Compressed sensing mri reconstruction using a generative adversarial network with a cyclic loss,” _IEEE Transactions on Medical Imaging_ , vol. 37, no. 6, pp. 1488–1497, 2018\.
  * [11] K. Hammernik, T. Klatzer, E. Kobler, M. P. Recht, D. K. Sodickson, T. Pock, and F. Knoll, “Learning a variational network for reconstruction of accelerated mri data,” _Magnetic Resonance in Medicine_ , vol. 79, no. 6, pp. 3055–3071, 2018.
  * [12] A. Sriram, J. Zbontar, T. Murrell, A. Defazio, C. L. Zitnick, N. Yakubova, F. Knoll, and P. Johnson, “End-to-end variational networks for accelerated mri reconstruction,” 2020.
  * [13] K. Zhang, K. A. Philbrick, G. M. Conte, Q. Huang, J. C. Cai, P. Rouzrokh, H. Hu, W. A. Jons, J. M. M. Jagtap, Y. Singh, E. G. Fuemmeler, J. D. Sobek, J. G. Pipe, and B. J. Erickson, “A dual-space variational network for mr imaging reconstruction,” _SIIM, CMIMI20_ , 2020.
  * [14] C. Dong, C. C. Loy, K. He, and X. Tang, “Image super-resolution using deep convolutional networks,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 38, no. 2, pp. 295–307, 2016.
  * [15] W. Shi, J. Caballero, F. Huszár, J. Totz, A. P. Aitken, R. Bishop, D. Rueckert, and Z. Wang, “Real-time single image and video super-resolution using an efficient sub-pixel convolutional neural network,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016, pp. 1874–1883.
  * [16] C. Ledig, L. Theis, F. Huszár, J. Caballero, A. Cunningham, A. Acosta, A. Aitken, A. Tejani, J. Totz, Z. Wang, and W. Shi, “Photo-realistic single image super-resolution using a generative adversarial network,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2017, pp. 105–114.
  * [17] B. Lim, S. Son, H. Kim, S. Nah, and K. M. Lee, “Enhanced deep residual networks for single image super-resolution,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition Workshops (CVPRW)_ , 2017, pp. 1132–1140.
  * [18] X. Wang, K. Yu, S. Wu, J. Gu, Y. Liu, C. Dong, Y. Qiao, and C. C. Loy, “Esrgan: Enhanced super-resolution generative adversarial networks,” in _The European Conference on Computer Vision Workshops (ECCVW)_ , September 2018.
  * [19] L. A. Gatys, A. S. Ecker, and M. Bethge, “Image style transfer using convolutional neural networks,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016, pp. 2414–2423.
  * [20] J. Johnson, A. Alahi, and L. Fei-Fei, “Perceptual losses for real-time style transfer and super-resolution,” in _European Conference on Computer Vision_ , 2016.
  * [21] Y. Blau, R. Mechrez, R. Timofte, T. Michaeli, and L. Zelnik-Manor, “The 2018 pirm challenge on perceptual image super-resolution,” 2019.
  * [22] A. S. Chaudhari, Z. Fang, F. Kogan, J. Wood, K. J. Stevens, E. K. Gibbons, J. H. Lee, G. E. Gold, and B. A. Hargreaves, “Super-resolution musculoskeletal mri using deep learning,” _Magnetic Resonance in Medicine_ , vol. 80, no. 5, pp. 2139–2154, 2018.
  * [23] C. Zhao, B. E. Dewey, D. L. Pham, P. A. Calabresi, D. S. Reich, and J. L. Prince, “Smore: A self-supervised anti-aliasing and super-resolution algorithm for mri using deep learning,” _IEEE Transactions on Medical Imaging_ , pp. 1–1, 2020.
  * [24] C. Zhao, A. Carass, B. E. Dewey, and J. L. Prince, “Self super-resolution for magnetic resonance images using deep networks,” in _2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018)_ , 2018, pp. 365–368.
  * [25] Y. Chen, Y. Xie, Z. Zhou, F. Shi, A. G. Christodoulou, and D. Li, “Brain mri super resolution using 3d deep densely connected neural networks,” in _2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018)_ , 2018, pp. 739–742.
  * [26] I. Sanchez and V. Vilaplana, “Brain mri super-resolution using 3d generative adversarial networks,” 2018.
  * [27] J. Wang, Y. Chen, Y. Wu, J. Shi, and J. Gee, “Enhanced generative adversarial network for 3d brain mri super-resolution,” in _2020 IEEE Winter Conference on Applications of Computer Vision (WACV)_ , 2020, pp. 3616–3625.
  * [28] J. Shi, Z. Li, S. Ying, C. Wang, Q. Liu, Q. Zhang, and P. Yan, “Mr image super-resolution via wide residual networks with fixed skip connection,” _IEEE Journal of Biomedical and Health Informatics_ , vol. 23, no. 3, pp. 1129–1140, 2019.
  * [29] J. Du, Z. He, L. Wang, A. Gholipour, Z. Zhou, D. Chen, and Y. Jia, “Super-resolution reconstruction of single anisotropic 3d mr images using residual convolutional neural network,” _Neurocomputing_ , vol. 392, pp. 209–220, 2020.
  * [30] W. Luo, Y. Li, R. Urtasun, and R. Zemel, “Understanding the effective receptive field in deep convolutional neural networks,” in _Advances in Neural Information Processing Systems_ , D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, Eds., vol. 29. Curran Associates, Inc., 2016, pp. 4898–4906.
  * [31] H. Greenspan, S. Peled, G. Oz, and N. Kiryati, “Mri inter-slice reconstruction using super-resolution,” in _Medical Image Computing and Computer-Assisted Intervention – MICCAI 2001_ , W. J. Niessen and M. A. Viergever, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2001, pp. 1204–1206.
  * [32] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” in _International Conference on Learning Representations_ , 2015.


# A Schmidt decomposition approach to quantum thermodynamics

André Malavazi<EMAIL_ADDRESS>Frederico Brito<EMAIL_ADDRESS>Instituto de Física de São Carlos, Universidade de São Paulo, C.P. 369,
13560-970 São Carlos, SP, Brasil

###### Abstract

The development of a self-consistent thermodynamic theory of quantum systems
is of fundamental importance for modern physics. Still, despite its essential
role in quantum science and technology, there is no unifying formalism for
characterizing the thermodynamics within general autonomous quantum systems,
and many fundamental open questions remain unanswered. Along these lines, most
current efforts and approaches restrict the analysis to particular scenarios
of approximative descriptions and semi-classical regimes. Here we propose a
novel approach to describe the thermodynamics of arbitrary bipartite
autonomous quantum systems based on the well-known Schmidt decomposition. This
formalism provides a simple, exact and symmetrical framework for expressing
the energetics between interacting systems, including scenarios beyond the
standard description regimes, such as strong coupling. We show that this
procedure allows a straightforward identification of local effective operators
suitable for characterizing the physical local internal energies. We also
demonstrate that these quantities naturally satisfy the usual thermodynamic
notion of energy additivity.

††preprint: APS/123-QED

During the last decades, there have been many efforts into the extension of
thermodynamics both to encompass microscopic systems out of equilibrium [1, 2]
and to include genuine quantum phenomena. Along with the state-of-the-art
capability of controlling fragile quantum systems in a wide variety of
physical platforms, this context has paved the way for the current strategic
endeavours to develop a self-consistent thermodynamic theory of quantum
systems [3, 4, 5, 6]. Interestingly, despite all recent extensive research and
advances in quantum thermodynamics [7, 8, 9, 10], some open questions
concerning central aspects of the theory are exceptionally notorious. In
particular, there are no universally accepted general definitions for the
quantum counterparts of the most basic thermodynamic quantities and,
therefore, no formalism consolidated for a quantum thermodynamic description
suitable for arbitrary quantum systems. This context poses a fundamental
restriction to many cases of interest and highlights a core issue of the
field.

While there is no ambiguity in defining proper thermodynamic quantities for
macroscopic systems in equilibrium, fundamental subtleties arise once one must
account for non-negligible interactions and intrinsically quantum properties
of systems out-of-equilibrium. In this sense, the identification of physical
local internal energy and quantum thermodynamic entropy is not as
straightforward as in the classical scenario since it is unclear how to
proceed in the presence of these notions. Although different definition
approaches have been reported in the literature [11, 12, 13, 14, 15, 16, 17,
18], this task remains elusive and far from consensus. Naturally, any attempt
to describe energy exchanges, such as work and heat, inherits this fundamental
uncertainty and brings its own conceptual issues and challenges [19, 20, 21,
22, 23, 24].

Additionally, most current approaches in quantum thermodynamics commonly rely
on a particular class of assumptions and approximations, which inevitably
limits their application to a small set of regimes. In particular, it is often
assumed weakly-coupled systems, Markovian dynamics, and semi-classical
external control [11, 25, 26, 27]. Despite being well-suitable for specific
scenarios, such assumptions clearly exclude several relevant situations.

Along these lines, it is clear that the current status of the theory is
restrictive and does not allow the description and characterization of
thermodynamic processes for general contexts, especially for autonomous
quantum systems and in the presence of strong coupling. In this respect, many
efforts have been attempting to expand quantum thermodynamics to these broader
conditions [12, 28, 13, 29, 30, 31, 32, 33, 34, 15]. Nevertheless, there is no
general framework suitable for arbitrary quantum systems.

This work is inserted in this context. Here, we introduce a novel approach for
analyzing the energetic exchange occurring between parts of a general
autonomous quantum system. This proposal is exact and based on the well-known
mathematical procedure of the Schmidt decomposition for bipartite systems.
Interestingly, despite being simple and providing a powerful statement, it is
still not explicitly explored in the context of quantum thermodynamics. This
framework will allow us to describe the dynamics and energetics within generic
interacting subsystems in a symmetrical fashion, i.e., regardless of their
individual properties, details, and dimension, they will be treated on equal
footing. In addition, it will not require any complementary hypotheses and
approximations, such as the commonly used ones concerning the Hamiltonian
structure, interaction regimes, and type of dynamics, i.e., strict energy
conservation, weak-coupling, and Markovianity. Formally, we will introduce
time-dependent local effective Hamiltonians that naturally embrace both their
respective bare ones and the contributions of the interaction term. These
elements will be identified as the representative operators for characterizing
the subsystem’s physical internal energies and will allow us to extend the
usual classical thermodynamic notion of energy additivity to general
interacting quantum systems. Such identification represents one of our main
results.

##  I Results

### System description

As a quantum universe, we consider a finite and isolated quantum system
composed of two interacting subsystems with dimensions $d^{(1)}$ and
$d^{(2)}$. In addition, let us assume that - without any loss of generality -
$d^{(1)}\leq d^{(2)}$. The whole system Hamiltonian is given by

$\hat{H}^{(0)}:=\hat{H}^{(1)}\otimes\hat{1}^{(2)}+\hat{1}^{(1)}\otimes\hat{H}^{(2)}+\hat{H}_{int},$
(1)

where $\hat{H}^{(k)}$ is the local bare Hamiltonian of subsystem $(k)$ with
$k=1,2$ and $\hat{H}_{int}$ is the term encompassing all interactions between
them.

At any time $t$, the universe’s quantum state is depicted by a pure state
represented by the state vector $|\Psi(t)\rangle$, whose dynamics is unitary
and governed by the usual Schrödinger equation
$i\hbar\frac{d}{dt}|\Psi(t)\rangle=\hat{H}^{(0)}|\Psi(t)\rangle$. The local
state description for each subsystem is expressed by the reduced density
matrix of the universe, such that $\hat{\rho}^{(1,2)}(t)\equiv
tr_{2,1}\\{|\Psi(t)\rangle\langle\Psi(t)|\\}$. No additional assumptions
regarding the universe, its constituents or the Hamiltonians will be
considered.

### Schmidt decomposition approach

Naturally, $|\Psi(t)\rangle$ can be written according to any conceivable
basis. Nevertheless, the well known Schmidt decomposition [35] provides the
following particular and convenient form

$|\Psi(t)\rangle=\sum_{j=1}^{d^{(1)}}\lambda_{j}(t)|\varphi_{j}^{(1)}(t)\rangle\otimes|\varphi_{j}^{(2)}(t)\rangle,$
(2)

for every instant $t$, where $\\{\lambda_{j}(t)\geq 0;j=1,...,d^{(1)}\\}$ and
$\\{|\varphi_{j}^{(k)}(t)\rangle;j=1,...,d^{(1)}\\}$ are the time-local
Schmidt coefficients and local Schmidt basis of subsystem $(k)$, respectively.
The representation above is compelling and useful for a number of reasons:
notice that despite a potential huge discrepancy between $d^{(1)}$ and
$d^{(2)}$ there is a single sum bounded by the smallest dimension in question,
by hypothesis $d^{(1)}$; apart from that, there is a one-to-one relation
between the Schmidt basis elements, i.e., for each Schmidt ket
$|\varphi_{j}^{(1)}(t)\rangle$ from subsystem $(1)$ there is a related ket
$|\varphi_{j}^{(2)}(t)\rangle$ from $(2)$; also, it is guaranteed that the
Schmidt coefficients are unambiguously defined, while the Schmidt basis are
unique up to eventual degenerate coefficients and a phase degree of freedom,
in the sense that equation (2) is invariant over simultaneous local phases
changes; additionally, it gives all necessary information for representing the
subsystem’s local density matrices $\hat{\rho}^{(k)}(t)$, i.e., their
simultaneous spectral decomposition is given by the Schmidt coefficients
squared $\\{\lambda_{j}^{2}(t)\\}_{j}$ and the Schmidt basis
$\\{|\varphi_{j}^{(k)}(t)\rangle\\}_{j}$, such that
$\hat{\rho}^{(k)}(t)=\sum_{j=1}^{d^{(1)}}\lambda_{j}^{2}(t)|\varphi_{j}^{(k)}(t)\rangle\langle\varphi_{j}^{(k)}(t)|$.
In particular, notice that both local state representations are - in general -
mixed state and have the same spectrum whenever the universe is pure.

The autonomous time evolution executed by the whole system at any time
interval $[t_{0},t_{1}]$ also induces the local dynamic changes, which, in
turn, are mainly non-unitary and strongly influenced by the presence of the
interaction term $\hat{H}_{int}$. Let us now focus on the changes in the
Schmidt decomposition elements and define the local dynamical maps
$\tilde{\mathcal{U}}^{(k)}:\mathcal{H}^{(k)}\rightarrow\mathcal{H}^{(k)}$ as
the operators characterizing the Schmidt basis dynamics, i.e.,
$|\varphi_{j}^{(k)}(t)\rangle=\tilde{\mathcal{U}}^{(k)}(t,t_{0})|\varphi_{j}^{(k)}(t_{0})\rangle$,
for any $t\geq t_{0}$, with $\underset{t\rightarrow
t_{0}}{\lim}\,\tilde{\mathcal{U}}^{(k)}(t,t_{0})=\hat{1}^{(k)}$. Since the
Schmidt basis at distinct times corresponds to different orthonormal basis for
the same Hilbert space, the relationship above is unitary, and one can define
the local time-translation generators
$\tilde{H}^{(k)}(t)=\tilde{H}^{(k)\dagger}(t)$ of the Schmidt basis
$\\{|\varphi_{j}^{(k)}(t)\rangle\\}_{j}$, such that
$i\hbar\frac{d}{dt}\tilde{\mathcal{U}}^{(k)}(t,t_{0})=\tilde{H}^{(k)}(t)\tilde{\mathcal{U}}^{(k)}(t,t_{0})$,
and, therefore,

$i\hbar\frac{d}{dt}|\varphi_{j}^{(k)}(t)\rangle=\tilde{H}^{(k)}(t)|\varphi_{j}^{(k)}(t)\rangle,$
(3)

for all $j$. Alternatively, $\tilde{H}^{(k)}(t)$ can be simply cast as

$\tilde{H}^{(k)}(t)\equiv
i\hbar\sum_{j=1}^{d^{(k)}}\left(\frac{d}{dt}|\varphi_{j}^{(k)}(t)\rangle\right)\langle\varphi_{j}^{(k)}(t)|.$
(4)

Still, we can go one step further and show that equation (4) above can be
directly related to the local bare Hamiltonian $\hat{H}^{(k)}$ as follows (see
Methods section)

$\tilde{H}^{(k)}(t)=\hat{H}^{(k)}+\hat{H}_{LS}^{(k)}(t)+\hat{H}_{X}^{(k)}(t),$
(5)

i.e., $\tilde{H}^{(k)}(t)$ can be split into the sum of three distinct
elements, including the bare Hamiltonian. The additional operators are
responsible for the time-dependency of $\tilde{H}^{(k)}(t)$, where
$\hat{H}_{LS}^{(k)}(t)$ is a general Lamb-shift like term, in the sense that
it is diagonal in the bare eigenbasis, i.e.,
$[\hat{H}_{LS}^{(k)}(t),\hat{H}^{(k)}]=0$ for all $t$, and
$\hat{H}_{X}^{(k)}(t)$ contains only non-diagonal elements. Note that in the
absence of interactions between the partitions, both subsystems would
individually evolve in time in a unitary fashion. Under these circumstances,
the role of the time-translation generator of the Schmidt basis would be
naturally assigned to the respective local bare Hamiltonians. That is, the
free evolution of the Schmidt basis would be simply
$|\bar{\varphi}_{j}^{(k)}(t)\rangle=e^{-\frac{i}{\hbar}\hat{H}^{(k)}(t-t_{0})}|\varphi_{j}^{(k)}(t_{0})\rangle$.
In this sense, it is clear that (i) the elements $\hat{H}_{LS}^{(k)}(t)$ and
$\hat{H}_{X}^{(k)}(t)$ are direct local byproducts of the interaction term and
vanish for non-interacting systems; (ii) equation (5) acts as an effective
Hamiltonian. Thus, from now on, $\tilde{H}^{(k)}(t)$ will be referred to as
the local effective Hamiltonian of subsystem $(k)$.

Also, one finds that the local effective Hamiltonians are responsible for
guiding the unitary component of the dynamics of $\hat{\rho}^{(k)}(t)$, i.e.,

$i\hbar\frac{d}{dt}\hat{\rho}^{(k)}(t)=[\tilde{H}^{(k)}(t),\hat{\rho}^{(k)}(t)]+\mathfrak{\hat{D}}_{t}^{(k)}\hat{\rho}^{(k)}(t),$
(6)

where
$\mathfrak{\hat{D}}_{t}^{(k)}\hat{\rho}^{(k)}(t)=i\hbar\sum_{j=1}^{d^{(1)}}\left(\frac{d}{dt}\lambda_{j}^{2}(t)\right)|\varphi_{j}^{(k)}(t)\rangle\langle\varphi_{j}^{(k)}(t)|$
is the non-unitary part that explicitly depends on the Schmidt coefficients
time dependency.

### Internal energies and additivity

By construction, in the Schrödinger picture, the time-translation generator of
isolated and closed quantum systems is the functional of energy. Thus, the
total internal energy, $U^{(0)}$, contained within the universe under
consideration is completely specified by the bare Hamiltonian $\hat{H}^{(0)}$,
i.e., one can readily identify

$U^{(0)}\equiv\langle\hat{H}^{(0)}\rangle$

where $\langle.\rangle\equiv\langle\Psi(t)|(.)|\Psi(t)\rangle$. Additionally,
since the bipartition is isolated, no energy flows inward or outward of the
system. Along these lines, it is easy to see that the unitary evolution of
$|\Psi(t)\rangle$ naturally guarantees internal energy conservation, such that

$\frac{d}{dt}U^{(0)}=0.$ (7)

Interestingly, it is unclear how to obtain a consistent and meaningful
analogous for the interacting subsystems. In general, the notion of local
internal energy is blurred by non-negligible interactions within the whole
system, i.e., the total internal energy can be written as the sum of the
individual local contributions of the expectation values of the bare
Hamiltonians and the interaction term,

$U^{(0)}=\langle\hat{H}^{(1)}\rangle(t)+\langle\hat{H}^{(2)}\rangle(t)+\langle\hat{H}_{int}\rangle(t).$

Thus, unless considering particular scenarios, it is clear that the
interaction is indispensable for the computation of the universe’s energy.
Yet, its role is uncertain when one attempts to establish local observable for
characterizing the subsystem’s internal energy. Along these lines, given that
the interaction term actively influences their local dynamics, it is
reasonable to expect that its contribution should be shared between them, in
one way or another.

In this context, it is desirable two relevant properties for any suitable
definition of local internal energies: (i) be obtained by local measurements,
i.e., associated with the expectation value of a local operator; (ii) be an
additive quantity. While the first condition guarantees a local description
and accessibility, the second also allows the intuitive picture of energy
flowing from one system to another without including energetic sinks or
sources, i.e., the sum of the local internal energies is a conserved quantity.
These features, of course, are not trivial, especially because the interaction
term acts on the whole Hilbert space, which means that it is a global property
per se. This fact, however, suggests that an effective approach for describing
local internal energy provides the most promising route. Otherwise, a global
picture would be necessary for characterizing the energetic flux, which is
impractical for most realistic scenarios.

We first remark that the local effective Hamiltonians presented earlier
satisfy both desired features for being interpreted as the representative
operators for characterizing the subsystem’s physical internal energies. Thus,
by identifying

$U^{(k)}(t):=\langle\tilde{H}^{(k)}(t)\rangle,$ (8)

the local internal energies become locally accessible properties, and one can
show that their sum is equal to the constant universe’s internal energy, such
that

$U^{(0)}=U^{(1)}(t)+U^{(2)}(t)$ (9)

and all energy flowing from subsystem $(1)$ is fully transferred to subsystem
$(2)$, i.e., no additional energetic source is required (see Methods section).

Note that such description is general, symmetrical and does not rely on
approximation, particular regimes and additional restrictive hypotheses
concerning the subsystem dynamics, the Hamiltonian structure and the coupling
strength. This identification represents one of our main results.

### Phase gauge

As mentioned earlier, the Schmidt basis is unique up to eventual degenerate
Schmidt coefficients and a phase choice freedom. While the former is
unimportant to our current purposes, the latter brings non-trivial
consequences to our formalism. From a global perspective, it is easy to see
that the simultaneous addition of local phases
$\\{\theta_{j}(t)\\}_{j}\in\mathbb{R}$, such that

$|\varphi_{j}^{\prime(k)}(t)\rangle=e^{(-1)^{k-1}i\theta_{j}(t)}|\varphi_{j}^{(k)}(t)\rangle,$
(10)

maintains the universe state form unchanged, i.e., the phases cancel out and
$|\Psi^{\prime}(t)\rangle=|\Psi(t)\rangle$. Such phase invariance expresses an
internal gauge freedom choice inbuilt in the Schmidt decomposition that
directly affects the identification of the local effective Hamiltonians, i.e.,
instead of describing the state vector with the basis
$\\{|\varphi_{j}^{(k)}(t)\rangle\\}_{j}$, one can use
$\\{|\varphi_{j}^{\prime(k)}(t)\rangle\\}_{j}$ and find its time-translation
generator $\tilde{H}^{\prime(k)}(t)$. Given equation (4), the gauge
transformation defined in equation (10) above implies the following expression
for $\tilde{H}^{\prime(k)}(t)$ in terms of $\tilde{H}^{(k)}(t)$ and
$\\{\theta_{j}(t)\\}_{j}$

$\tilde{H}^{\prime(k)}(t)=\tilde{H}^{(k)}(t)+(-1)^{k}\hbar\sum_{j=1}^{d^{(k)}}\left(\frac{d\theta_{j}(t)}{dt}\right)|\varphi_{j}^{(k)}(t)\rangle\langle\varphi_{j}^{(k)}(t)|.$
(11)

Naturally, its expectation value also transforms accordingly, such that

$\langle\tilde{H}^{\prime(k)}(t)\rangle=\langle\tilde{H}^{(k)}(t)\rangle+(-1)^{k}\hbar\sum_{j=1}^{d^{(k)}}\lambda_{j}^{2}(t)\left(\frac{d\theta_{j}(t)}{dt}\right).$
(12)

Thus, it is clear that from a local point of view these phases are relevant
and should be carefully scrutinized.

Note that a gauge change adds an extra diagonal term on the Schmidt basis that
only depends on the time derivative of the phases. In general, these
additional quantities change the local effective Hamiltonian structure in a
way that both their eigenbasis and eigenvalues are affected, which also
implies that the spectral gaps are not necessarily invariant. Still, as
expected, such transformations maintain the description of the local density
matrices and their dynamical equation (6) unchanged, given that
$[|\varphi_{j}^{(k)}(t)\rangle\langle\varphi_{j}^{(k)}(t)|,\hat{\rho}^{(k)}(t)]=0$
for all $j$ and $t$. More importantly, these terms shift the expectation
values $\langle\tilde{H}^{(k)}(t)\rangle$ in such a way the change obtained by
subsystem $(1)$ is compensated by the one acquired by subsystem $(2)$. Hence,
both the universe’s internal energy $U^{(0)}$ and the additivity property are
guaranteed to be invariant over these modifications, i.e.,
$U^{\prime(0)}=U^{(0)}$ and
$\langle\tilde{H}^{\prime(1)}(t)\rangle+\langle\tilde{H}^{\prime(2)}(t)\rangle=\langle\tilde{H}^{(1)}(t)\rangle+\langle\tilde{H}^{(2)}(t)\rangle$.

Hence, given the freedom of choosing the Schmidt basis for writing the state
vector, it is clear that for the same bipartite quantum system, several local
effective Hamiltonians are suitable for consistently describing the
subsystem’s dynamics, any observable physical quantity and the universe
energetics. However, despite $\tilde{H}^{(k)}(t)$ possessing the desired
properties for successfully fulfilling the role of physical local internal
energies, such freedom also brings a significant ambiguity in identifying a
single pair of preferential operators for quantifying them.

Notably, such mathematical freedom of adding local phases is not necessarily
physically consistent for characterizing energy in general when seeking its
association with the time-translation generator of the dynamics 111In a
private communication, Prof. Shang-Yung Wang presented an approach based on
the introduction of gauge potentials such that the local effective Hamiltonian
is gauge invariant. Even though this feature is quite neat, it does not
preserve the association of the Hamiltonian as the time-translation generator,
differing from the approach adopted here.. As mentioned earlier, if there is
no interaction between the bipartitions, they would behave independently as
isolated quantum systems. In this scenario, the Schmidt basis dynamics would
be described by
$i\hbar\frac{d}{dt}|\bar{\varphi}_{j}^{(k)}(t)\rangle=\hat{H}^{(k)}|\bar{\varphi}_{j}^{(k)}(t)\rangle$,
and the individual local internal energies would be readily identified as
$U^{(k)}\equiv\langle\hat{H}^{(k)}\rangle$, where $\frac{d}{dt}U^{(k)}=0$.
Along these lines, once considering the interaction term negligible, similar
conclusions should be obtained regardless of the chosen gauge. In this
respect, if $\\{|\varphi_{j}^{(k)}(t)\rangle\\}_{j}$ is indeed the Schmidt
basis such that
$|\varphi_{j}^{(k)}(t)\rangle\rightarrow|\bar{\varphi}_{j}^{(k)}(t)\rangle$
when $\hat{H}_{int}\rightarrow 0$, then it is guaranteed that
$\tilde{H}^{(k)}(t)\rightarrow\hat{H}^{(k)}$ and, since the phases are
arbitrary and independent of the interaction term,
$\tilde{H}^{\prime(k)}(t)\rightarrow\hat{H}^{(k)}+(-1)^{k}\hbar\sum_{j=1}^{d^{(k)}}\left(\frac{d\theta_{j}(t)}{dt}\right)|\bar{\varphi}_{j}^{(k)}(t)\rangle\langle\bar{\varphi}_{j}^{(k)}(t)|$.
Notice that, while the former is the time-translation generator of the basis
$\\{|\bar{\varphi}_{j}^{(k)}(t)\rangle\\}_{j}$, the latter is the one relative
to the basis
$\\{e^{(-1)^{k-1}i\theta_{j}(t)}|\bar{\varphi}_{j}^{(k)}(t)\rangle\\}_{j}$.
More importantly, it is clear that, in general, not all local effective
Hamiltonians (phases) are suitable for correctly describing local energy
measurements in the limit behaviour of non-interacting systems. Nevertheless,
this is achieved iff $\frac{d}{dt}\theta_{j}(t)=\alpha$ for all $j$, since
additive constants only shift the bare Hamiltonian spectrum. Thus, as long all
the phases are equal linear functions of time, we have a physically consistent
set of gauges satisfying
$\tilde{H}^{\prime(k)}(t)=\tilde{H}^{(k)}(t)+(-1)^{k}\hbar\alpha\hat{1}^{(k)}$
and
$\langle\tilde{H}^{\prime(k)}(t)\rangle=\langle\tilde{H}^{(k)}(t)\rangle+(-1)^{k}\hbar\alpha$
in a way that the remaining freedom is also guaranteed to maintain the
spectral gaps of the local effective Hamiltonians invariant. Despite being
local operators, observe that such a procedure requires knowledge of the
“whole” in order to determine the “parts” energetics, i.e., obtaining the
physical local effective Hamiltonians demands knowledge about the interaction
of the parts (a non-local property) to fix the correct physical phase gauge.
Hence, in order to construct a local and consistent energy description for the
parts, one cannot rely solely on local features (see Reference [37] for a
recent discussion).

In summary, by restricting ourselves to gauges such that
$\left\\{\frac{d}{dt}\theta_{j}(t)=\alpha\right\\}_{j}$, where
$\alpha\in\mathbb{R}$, one can identify the expectation values
$\langle\tilde{H}^{(k)}(t)\rangle$ as the local physical internal energies up
to an additive constant, $U^{(k)}(t):=\langle\tilde{H}^{(k)}(t)\rangle$ and
$U^{\prime(k)}(t)=U^{(k)}(t)+(-1)^{k}\hbar\alpha$. In this sense, even though
different gauges provide distinct absolute internal energy values, energy
measurements differences remain identical.

##  II Conclusion

In this work, we introduced a novel formalism suitable for describing the
energetics within arbitrary autonomous quantum systems. Despite being
presented for pure quantum states, the generalization for mixed ones is
straightforward. The formal procedure is based on the well-known Schmidt
decomposition given by equation (2), which provides the basis for the
identification of the local effective Hamiltonians defined by equation (4). We
highlighted the fact that these operators possess the desired properties for
being considered suitable candidates for characterizing the subsystem’s
internal energies, i.e., their expectation values are local quantities and
naturally satisfy the classical thermodynamic notion of energy additivity, as
shown in equation (9). In contrast with currently used approaches, such a
framework is exact, treats both partitions on equal footing and does not rely
on approximations and additional restrictive hypotheses, i.e., particular
coupling regimes, convenient Hamiltonian structures, semi-classical
description and specific types of dynamics.

Along with such an identification, we verified the existence of an internal
phase gauge freedom corresponding to the phases choices
$\\{\theta_{j}(t)\\}_{j}$ within the Schmidt decomposition procedure. As shown
in equations (11, 12), such arbitrary choice affects the structure of the
local effective Hamiltonians and, in general, manifests as a time-dependent
contribution to their expectation values. Nevertheless, despite the
mathematical freedom, physical consistency during gauges transformations is
only achieved for phases such that $\frac{d}{dt}\theta_{j}(t)=\alpha$ for all
$j$, where $\alpha\in\mathbb{R}$. In this case, we have
$\tilde{H}^{\prime(k)}(t)=\tilde{H}^{(k)}(t)+(-1)^{k}\hbar\alpha\hat{1}^{(k)}$
and
$\langle\tilde{H}^{\prime(k)}(t)\rangle=\langle\tilde{H}^{(k)}(t)\rangle+(-1)^{k}\hbar\alpha$,
and the limit behaviour
$\tilde{H}^{\prime(k)}(t)\rightarrow\hat{H}^{(k)}+(-1)^{k}\hbar\alpha\hat{1}^{(k)}$
is satisfied for $\hat{H}_{int}\rightarrow 0$. Thus, by identifying the local
internal energies as $U^{(k)}(t):=\langle\tilde{H}^{(k)}(t)\rangle$, different
gauges would provide distinct absolute values for these quantities but
equivalent spectral gaps. In this sense, the energy shift
$(-1)^{k}\hbar\alpha$ is analogous to the classical thermodynamic freedom in
the definition of internal energy [38].

In the current status of quantum thermodynamics, where there is no general
methodology suitable for describing the energetic exchanges of arbitrary
interacting quantum systems on equal footing, the procedure devised here
represents a viable alternative to fill this important gap and a promising
route to characterize general quantum thermodynamic processes. In particular,
providing means to describe scenarios that do not fall under the standard
description regimes, especially for strongly coupled systems, and in the
absence of external classical agents.

##  III Methods

### Local effective Hamiltonian

Given the following spectral decomposition for the local bare Hamiltonian

$\hat{H}^{(k)}\equiv\sum_{j=1}^{d^{(k)}}b_{j}^{(k)}|b_{j}^{(k)}\rangle\langle
b_{j}^{(k)}|,$ (13)

where $\\{b_{j}^{(k)}\\}_{j}$ and $\\{|b_{j}^{(k)}\rangle\\}_{j}$ are its
respective bare eigenenergies and eigenbasis. One can define the projections
$\langle
b_{j}^{(k)}|\varphi_{l}^{(k)}(t)\rangle:=r_{jl}^{(k)}(t)e^{-\frac{i}{\hbar}b_{j}^{(k)}t}$,
such that

$\langle
b_{j}^{(k)}|\frac{d}{dt}|\varphi_{l}^{(k)}(t)\rangle=\left(\frac{d}{dt}r_{jl}^{(k)}(t)\right)e^{-\frac{i}{\hbar}b_{j}^{(k)}t}\\\
-\frac{i}{\hbar}b_{j}^{(k)}r_{jl}^{(k)}(t)e^{-\frac{i}{\hbar}b_{j}^{(k)}t}.$

Also, given the orthonormality $\langle
b_{\alpha}^{(k)}|b_{\beta}^{(k)}\rangle=\delta_{\alpha\beta}$, it is easy to
see that

$\sum_{l=1}^{d^{(k)}}r_{\alpha l}^{(k)}(t)\left(r_{\beta
l}^{(k)}(t)\right)^{*}e^{\frac{i}{\hbar}\left(b_{\beta}^{(k)}-b_{\alpha}^{(k)}\right)t}=\delta_{\alpha\beta}.$

Finally, by casting $\tilde{H}^{(k)}(t)$ in the bare eigenbasis representation
and using the previous relations, one can rewrite the local effective
Hamiltonian according to equation (5), i.e.,

$\tilde{H}^{(k)}(t)=\hat{H}^{(k)}+\hat{H}_{LS}^{(k)}(t)+\hat{H}_{X}^{(k)}(t),$

where

$\hat{H}_{LS}^{(k)}(t)\equiv
i\hbar\sum_{j=1}^{d^{(k)}}\left(\sum_{l=1}^{d^{(k)}}\left(\frac{d}{dt}r_{jl}^{(k)}(t)\right)r_{jl}^{(k)*}(t)\right)|b_{j}^{(k)}\rangle\langle
b_{j}^{(k)}|,$

and

$\hat{H}_{X}^{(k)}(t)\equiv i\hbar\sum_{j=1}^{d^{(k)}}\sum_{m\neq
j}^{d^{(k)}}\left(\sum_{l=1}^{d^{(k)}}\frac{d}{dt}r_{jl}^{(k)}(t)r_{ml}^{(k)*}(t)\right)\times\\\
e^{\frac{i}{\hbar}\left(b_{m}^{(k)}-b_{j}^{(k)}\right)t}|b_{j}^{(k)}\rangle\langle
b_{m}^{(k)}|.$

### Local internal energies additivity

Given equation (2) and (3) for the Schmidt decomposition and the Schmidt basis
dynamics, we have the following equation

$i\hbar\frac{d}{dt}|\Psi(t)\rangle=\left(\tilde{H}^{(1)}(t)+\tilde{H}^{(2)}(t)\right)|\Psi(t)\rangle\\\
+\sum_{j=1}^{d^{(1)}}\left(i\hbar\frac{d}{dt}\lambda_{j}(t)\right)|\varphi_{j}^{(1)}(t)\rangle\otimes|\varphi_{j}^{(2)}(t)\rangle.$

Thus, since $i\hbar\frac{d}{dt}|\Psi(t)\rangle=\hat{H}^{(0)}|\Psi(t)\rangle$,
we obtain

$\langle\Psi(t)|\hat{H}^{(0)}|\Psi(t)\rangle=\langle\Psi(t)|\left(\tilde{H}^{(1)}(t)+\tilde{H}^{(2)}(t)\right)|\Psi(t)\rangle\\\
+\langle\Psi(t)|\sum_{j=1}^{d^{(1)}}\left(i\hbar\frac{d}{dt}\lambda_{j}(t)\right)|\varphi_{j}^{(1)}(t)\rangle\otimes|\varphi_{j}^{(2)}(t)\rangle.$

However, notice that due to normalization of $|\Psi(t)\rangle$, the second
contribution is necessarily null. Therefore,

$\langle\hat{H}^{(0)}\rangle=\langle\tilde{H}^{(1)}(t)\rangle+\langle\tilde{H}^{(2)}(t)\rangle=U^{(0)}.$
(14)

##  IV Acknowledgements

AM and FB acknowledge financial support in part by the Coordenação de
Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code
001. FB is also supported by the Brazilian National Institute for Science and
Technology of Quantum Information (INCT- IQ). The authors would like to thank
Prof. Shang-Yung Wang for the fruitful discussions.

## References

  * Jarzynski [1997] C. Jarzynski, Nonequilibrium equality for free energy differences, Physical Review Letters 78, 2690 (1997), 10.1103/PhysRevLett.78.2690.
  * Seifert [2012] U. Seifert, Stochastic thermodynamics, fluctuation theorems and molecular machines, Reports on Progress in Physics 75, 126001 (2012), 10.1088/0034-4885/75/12/126001.
  * Gemmer  _et al._ [2009] J. Gemmer, M. Michel, and G. Mahler,  _Quantum thermodynamics_ , Lecture notes in physics, v. 784 (Springer, Cham, 2009) 10.1007/978-3-540-70510-9.
  * Vinjanampathy and Anders [2016] S. Vinjanampathy and J. Anders, Quantum thermodynamics, Contemporary Physics 57, 545 (2016), 10.1080/00107514.2016.1201896.
  * Binder  _et al._ [2018] F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso,  _Thermodynamics in the quantum regime_ , Fundamental theories of physics, v. 195 (Springer, Cham, 2018) 10.1007/978-3-319-99046-0.
  * Deffner and Campbell [2019] S. Deffner and S. Campbell,  _Quantum thermodynamics_ (Morgan $\&$ Claypool Publishers, San Rafael, CA, 2019) 10.1088/2053-2571/ab21c6.
  * Batalhão  _et al._ [2014] T. B. Batalhão, A. M. Souza, L. Mazzola, R. Auccaise, R. S. Sarthour, I. S. Oliveira, J. Goold, G. De Chiara, M. Paternostro, and R. M. Serra, Experimental reconstruction of work distribution and study of fluctuation relations in a closed quantum system, Physical Review Letters 113, 140601 (2014), 10.1103/PhysRevLett.113.140601.
  * An  _et al._ [2015] S. An, J.-N. Zhang, M. Um, D. Lv, Y. Lu, J. Zhang, Z.-Q. Yin, H. Quan, and K. Kim, Experimental test of the quantum jarzynski equality with a trapped-ion system, Nature Physics 11, 193 (2015), 10.1038/nphys3197.
  * Roßnagel  _et al._ [2016] J. Roßnagel, S. T. Dawkins, K. N. Tolazzi, O. Abah, E. Lutz, F. Schmidt-Kaler, and K. Singer, A single-atom heat engine, Science 352, 325 (2016), 10.1126/science.aad6320.
  * Klatzow  _et al._ [2019] J. Klatzow, J. N. Becker, P. M. Ledingham, C. Weinzetl, K. T. Kaczmarek, D. J. Saunders, J. Nunn, I. A. Walmsley, R. Uzdin, and E. Poem, Experimental demonstration of quantum effects in the operation of microscopic heat engines, Physical Review Letters 122, 110601 (2019), 10.1103/PhysRevLett.122.110601.
  * Alicki [1979] R. Alicki, The quantum open system as a model of the heat engine, Journal of Physics A: mathematical and general 12, L103 (1979), 10.1088/0305-4470/12/5/007.
  * Weimer  _et al._ [2008] H. Weimer, M. J. Henrich, F. Rempp, H. Schröder, and G. Mahler, Local effective dynamics of quantum systems: a generalized approach to work and heat, EPL (Europhysics Letters) 83, 30008 (2008), 10.1209/0295-5075/83/30008.
  * Alipour  _et al._ [2016] S. Alipour, F. Benatti, F. Bakhshinezhad, M. Afsary, S. Marcantoni, and A. T. Rezakhani, Correlations in quantum thermodynamics: heat, work, and entropy production, Scientific Reports 6, 1 (2016), 10.1038/srep35568.
  * Valente  _et al._ [2018] D. Valente, F. Brito, R. Ferreira, and T. Werlang, Work on a quantum dipole by a single-photon pulse, Optics Letters 43, 2644 (2018), 10.1364/OL.43.002644.
  * Colla and Breuer [2022] A. Colla and H.-P. Breuer, Open-system approach to nonequilibrium quantum thermodynamics at arbitrary coupling, Physical Review A 105, 052216 (2022), 10.1103/PhysRevA.105.052216.
  * Esposito  _et al._ [2010] M. Esposito, K. Lindenberg, and C. Van den Broeck, Entropy production as correlation between system and reservoir, New Journal of Physics 12, 013013 (2010), 10.1088/1367-2630/12/1/013013.
  * Polkovnikov [2011] A. Polkovnikov, Microscopic diagonal entropy and its connection to basic thermodynamic relations, Annals of Physics 326, 486 (2011), 10.1016/j.aop.2010.08.004.
  * Šafránek  _et al._ [2019] D. Šafránek, J. M. Deutsch, and A. Aguirre, Quantum coarse-grained entropy and thermodynamics, Physical Review A 99, 010101 (2019), 10.1103/PhysRevA.99.010101.
  * Talkner  _et al._ [2007] P. Talkner, E. Lutz, and P. Hänggi, Fluctuation theorems: work is not an observable, Physical Review E 75, 050102 (2007), 10.1103/PhysRevE.75.050102.
  * Elouard  _et al._ [2017] C. Elouard, D. A. Herrera-Martí, M. Clusel, and A. Auffeves, The role of quantum measurement in stochastic thermodynamics, npj Quantum Information 3, 1 (2017), 10.1038/s41534-017-0008-4.
  * Perarnau-Llobet  _et al._ [2017] M. Perarnau-Llobet, E. Bäumer, K. V. Hovhannisyan, M. Huber, and A. Acin, No-go theorem for the characterization of work fluctuations in coherent quantum systems, Physical Review Letters 118, 070601 (2017), 10.1103/PhysRevLett.118.070601.
  * Ahmadi  _et al._ [2021] B. Ahmadi, S. Salimi, and A. S. Khorashad, Refined definitions of heat and work in quantum thermodynamics, arXiv preprint arXiv:1912.01983 (2021), 10.48550/arXiv.1912.01983.
  * Alipour  _et al._ [2022] S. Alipour, A. T. Rezakhani, A. Chenu, A. del Campo, and T. Ala-Nissila, Entropy-based formulation of thermodynamics in arbitrary quantum evolution, Physical Review A 105, L040201 (2022), 10.1103/PhysRevA.105.L040201.
  * Lima Bernardo [2020] B. Lima Bernardo, Unraveling the role of coherence in the first law of quantum thermodynamics, Physical Review E 102, 062152 (2020), 10.1103/PhysRevE.102.062152.
  * Kosloff [2013] R. Kosloff, Quantum thermodynamics: a dynamical viewpoint, Entropy 15, 2100 (2013), 10.3390/e15062100.
  * Kosloff [2019] R. Kosloff, Quantum thermodynamics and open-systems modeling, Journal of chemical physics 150, 204105 (2019), 10.1063/1.5096173.
  * Dann and Kosloff [2021a] R. Dann and R. Kosloff, Open system dynamics from thermodynamic compatibility, Physical Review Research 3, 023006 (2021a), 10.1103/PhysRevResearch.3.023006.
  * Hossein-Nejad  _et al._ [2015] H. Hossein-Nejad, E. J. O’Reilly, and A. Olaya-Castro, Work, heat and entropy production in bipartite quantum systems, New Journal of Physics 17, 075014 (2015), 10.1088/1367-2630/17/7/075014.
  * Dann and Kosloff [2021b] R. Dann and R. Kosloff, Quantum thermo-dynamical construction for driven open quantum systems, Quantum 5, 590 (2021b), 10.22331/q-2021-11-25-590.
  * Esposito  _et al._ [2015] M. Esposito, M. A. Ochoa, and M. Galperin, Nature of heat in strongly coupled open quantum systems, Physical Review B 92, 235440 (2015), 10.1103/PhysRevB.92.235440.
  * Dou  _et al._ [2018] W. Dou, M. A. Ochoa, A. Nitzan, and J. E. Subotnik, Universal approach to quantum thermodynamics in the strong coupling regime, Physical Review B 98, 134306 (2018), 10.1103/PhysRevB.98.134306.
  * Perarnau-Llobet  _et al._ [2018] M. Perarnau-Llobet, H. Wilming, A. Riera, R. Gallego, and J. Eisert, Strong coupling corrections in quantum thermodynamics, Physical Review Letters 120, 120602 (2018), 10.1103/PhysRevLett.120.120602.
  * Strasberg [2019] P. Strasberg, Repeated interactions and quantum stochastic thermodynamics at strong coupling, Physical Review Letters 123, 180604 (2019), 10.1103/PhysRevLett.123.180604.
  * Rivas [2020] Á. Rivas, Strong coupling thermodynamics of open quantum systems, Physical Review Letters 124, 160601 (2020), 10.1103/PhysRevLett.124.160601.
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum computation and quantum information_ (Cambridge University Press, Cambridge, UK, 2010).
  * Note [1] In a private communication, Prof. Shang-Yung Wang presented an approach based on the introduction of gauge potentials such that the local effective Hamiltonian is gauge invariant. Even though this feature is quite neat, it does not preserve the association of the Hamiltonian as the time-translation generator, differing from the approach adopted here.
  * Neves and Brito [2022] L. R. T. Neves and F. Brito, A constraint on local definitions of quantum internal energy, arXiv preprint arXiv:2205.04457 (2022), 10.48550/arXiv.2205.04457.
  * Kondepudi and Prigogine [2015] D. Kondepudi and I. Prigogine,  _Modern thermodynamics_ , 2nd ed. (John Wiley, New York, 2015).


# Symmetry-preserving contact interaction model for heavy-light mesons

F. E. Serna  M. A. Brito  G. Krein

###### Abstract

We use a symmetry-preserving regularization method of ultraviolet divergences
in a vector-vector contact interaction model for low-energy QCD. The contact
interaction is a representation of nonperturbative kernels used Dyson-
Schwinger and Bethe-Salpeter equations. The regularization method is based on
a subtraction scheme that avoids standard steps in the evaluation of divergent
integrals that invariably lead to symmetry violation. Aiming at the study of
heavy-light mesons, we have implemented the method to the pseudoscalar $\pi$
and $K$ mesons. We have solved the Dyson-Schwinger equation for the $u,d$ and
$s$ quark propagators, and obtained the bound-state Bethe-Salpeter amplitudes
in a way that the Ward-Green-Takahashi identities reflecting global symmetries
of the model are satisfied for arbitrary routing of the momenta running in
loop integrals.

###### Keywords:

Regularization of ultraviolet divergences, Ward-Green-Takahashi identities,
Dyson-Schwinger and Bethe-Salpeter equations

###### :

11.10.St, 11.15.Tk, 11.30.-j, 14.40.Df, 14.40.Lb, 24.85.+p

##  1 Introduction

Contact fermionic interactions find widespread applications in hadron physics
with Nambu–Jona-Lasinio (NJL) Nambu:1961tp  type of models. A great deal of
qualitative insight on the phenomenon of hadron mass generation via dynamical
chiral symmetry breaking (D$\chi$SB) and the role of the pions as the
associated (quasi) Goldstone bosons has been gleaned from such models – for
reviews, see Refs. Vogl:1991qt ; Klevansky:1992qe ; Hatsuda:1994pi ;
Bijnens:1995ww . The lack of confinement and nonrenormalizabilty are the major
weaknesses of contact fermionic models. The nonrenormalizability carries with
it the danger of introducing gross violations of global symmetries due to the
regularization procedure; ambiguities arising from momentum shifts in
divergent integrals are the main cause of the problems. Despite well known,
practitioners at large inexplicably ignore these problems. However, a new
perspective has emerged recently with the implementation
GutierrezGuerrero:2010md ; Roberts:2010rn ; Roberts:2011wy  of a confining,
symmetry-preserving treatment of a vector-vector contact interaction as a
representation of the gluon’s two-point Schwinger function used in kernels of
Dyson-Schwinger (DS) equations. By introducing a mechanism that ensures the
absence of quark production thresholds, a feature of a confining theory, and
embedding the interaction in a global-symmetry-preserving, rainbow-ladder
truncation framework of the DS and BS equations Munczek:1994zz ; Bender:1996bb
, the scheme has enhanced the capacity of a contact interaction to describe in
a unified manner a diverse array of nonperturbative QCD phenomena – for a
recent review, see Ref. Cloet:2013jya

In this communication we examine the contact-interaction model within the
perspective of a regularization scheme that allows to separate symmetry-
offending parts in BS amplitudes in a way independent of choices of momentum
routing in divergent integrals, so that the Ward-Green-Takahashi (WGT)
identities reflecting global symmetries of the model are preserved by the
regularization. The scheme is inspired in a method introduced in Refs. Batt-
thesis ; Battistel:1998tj  to manipulate divergent Feynman integrals without
specification of an explicit regulator. The method has been used with the NJL
model in vacuum Battistel:2008fd ; Battistel:2013cja  and at finite
temperature and baryon density Farias:2006cs ; Farias:2005cr ; Farias:2007zz .
We formulate the problem generically with the aim of using it for heavy-light
mesons, and assess its reliability presenting numerical results for the $\pi$
and $K$ mesons.

##  2 DS and BS equations and WGT identities

We start with a brief review of the basic elements of the contact-interaction
scheme introduced in Ref. GutierrezGuerrero:2010md . We consider the
inhomogeneous BSE for a quark and antiquark state of total momentum $P$ (here
and in the following we omit renormalization constants)

$\displaystyle\left[\Gamma_{\cal M}(k;P)\right]_{AB}$ $\displaystyle=$
$\displaystyle{\cal
M}_{AB}+\int_{q}\left[K(k,q;P)\right]_{AC,DB}\left[S(q_{+})\Gamma_{\cal
M}(q;P)S(q_{-})\right]_{CD},$ (1)

where $\int_{q}\equiv\int d^{4}q/(2\pi)^{4}$, ${\cal M}$ gives the Dirac
spinor structure of the state, $K(q,k;P)$ is the fully amputated quark-
antiquark scattering kernel; $A,B,\cdots$ denote collectively color, flavor,
and spinor indices; $q_{\pm}=q\pm\eta_{\pm}P$, with $\eta_{+}+\eta_{-}=1$ and
$q$ is the relative momentum. $S(k)$ is the dressed-quark propagator given by
a DSE; for a given flavor $f$ the general form of this DSE is

$\displaystyle S_{f}(k)^{-1}=i\gamma\cdot
k+m_{f}+\int_{q}\,g^{2}D_{\mu\nu}(k-q)\,\frac{\lambda^{a}}{2}\gamma_{\mu}S_{f}(q)\frac{\lambda^{a}}{2}\Gamma^{f}_{\nu}(q,k),$
(2)

where $m_{f}$ is the current-quark mass. Here we will be interested in the
flavor-nonsinglet axial-vector $\Gamma^{lh}_{5\mu}(k;P)$ and pseudoscalar
$\Gamma^{lh}_{\rm PS}(k;P)$ amplitudes for a quark of flavor $l$ and antiquark
of flavor $h$, for which ${\cal M}_{5\mu}=\gamma_{5}\gamma_{\mu}$ and ${\cal
M}_{PS}=\gamma_{5}$ respectively. Translation invariance requires that
physical observables cannot depend on the choice of the momentum routing in
quark propagators in Eq. (1); that is, physical results must be independent of
$\eta_{\pm}$.

Associated with $\Gamma^{lh}_{5\mu}$ is the Ward-Green-Takahashi (WGT)
identity

$\displaystyle
P_{\mu}\Gamma_{5\mu}^{lh}(k;P)=S_{l}^{-1}(k_{+})i\gamma_{5}+i\gamma_{5}S_{h}^{-1}(k_{-})-\,i\,(m_{l}+m_{h})\,\Gamma^{lh}_{\rm
PS}(k;P)\,,$ (3)

where $\Gamma^{lh}_{PS}(k;P)$ is the pseudoscalar vertex. Pseudoscalar meson
bound states are obtained from the solution of the homogeneous equation for
$\Gamma^{lh}_{\rm PS}(k;P)$

$\displaystyle\left[\Gamma^{lh}_{\rm
PS}(k;P)\right]_{AB}=\int_{q}\left[K(k,q;P)\right]_{AC,DB}\left[S_{l}(q_{+})\Gamma^{lh}_{\rm
PS}(q;P)S_{h}(q_{-})\right]_{CD},$ (4)

where here $A,B,\cdots$ denote color and spinor indices only. The general form
of the of $\Gamma^{lh}_{\rm PS}(k;P)$ is

$\displaystyle\Gamma^{lh}_{\rm PS}(k,P)=\gamma_{5}\bigl{[}iE^{lh}_{\rm
PS}+\gamma\cdot P\,F^{lh}_{\rm PS}+\gamma\cdot k\,G^{lh}_{\rm
PS}+\,\sigma_{\mu\nu}k_{\mu}P_{\nu}\,H^{lh}_{\rm PS}\bigr{]},$ (5)

where $E^{lh}_{\rm PS},F^{lh}_{\rm PS},\cdots$ are functions of $k$, $P$ and
$k\cdot P$. The meson mass, $m_{\rm PS}$, is the eigenvalue $P^{2}=-m^{2}_{\rm
PS}$ that solves Eq. (4).

The contact-interaction scheme introduced in Ref. GutierrezGuerrero:2010md
amounts to make the following replacements in Eq. (2) and the construction of
scattering kernels of all BSEs

$g^{2}D_{\mu\nu}(k-q)\rightarrow\frac{4\pi\alpha_{IR}}{m^{2}_{G}}\delta_{\mu\nu}\
\hskip 14.22636pt{\rm and}\hskip
14.22636pt\Gamma^{a}_{\mu}\rightarrow\frac{\lambda^{a}}{2}\gamma_{\mu}\,,$ (6)

where $m_{G}$ is a gluon mass-scale and $\alpha_{IR}$ is a fitting parameter.
This means that the scattering kernel in all BSEs is

$\left[K(k,q;P)\right]_{AC,DB}=-\frac{4\pi\alpha_{IR}}{m^{2}_{G}}\left(\frac{\lambda^{a}}{2}\gamma_{\mu}\right)_{AC}\left(\frac{\lambda^{a}}{2}\gamma_{\mu}\right)_{DB}.$
(7)

A feature of the momentum independence of the contact interaction is that the
DS and BS equations become nonrenormalizable. This means that mass-scale
parameters introduced with the regularization of divergent integrals cannot be
removed from the calculations and need to be fixed phenomenologically. Another
feature, as mentioned previously, is that regularization of divergent
integrals carries the danger of symmetry violation, in that WGT identities,
like the one in Eq. (3), are not maintained, even when Poincaré-invariant
regularization schemes are employed.

Let us consider the DSE and the homogeneous pseudoscalar BSE, Eqs. (2) and
(4), with the contact interaction. In this case, $S_{f}(k)^{-1}=i\gamma\cdot
k+M_{f}$, with

$\displaystyle
M_{f}=m_{f}+\frac{64\pi\alpha_{IR}}{3m^{2}_{G}}\,\int_{q}\frac{M_{f}}{q^{2}+M^{2}_{f}}\,.$
(8)

In addition, $G^{lh}_{\rm PS}=H^{lh}_{\rm PS}=0$, so the Bethe-Salpeter
amplitude Eq. (5) takes the form

$\displaystyle\Gamma^{lh}_{\rm PS}(P)=\gamma_{5}\left[iE^{lh}_{\rm
PS}(P)+\frac{1}{2M_{lh}}\gamma\cdotp P\,F^{lh}_{\rm PS}(P)\right]~{},$ (9)

where $M_{lh}=M_{l}M_{h}/(M_{l}+M_{h})$. Then, the BSE can be written in the
matrix form

$\displaystyle\left[\begin{array}[]{c}E^{lh}_{\rm PS}(P)\\\\[5.69046pt]
F^{lh}_{\rm PS}(P)\end{array}\right]$ $\displaystyle=$
$\displaystyle\frac{4\pi\alpha_{IR}}{3m^{2}_{G}}\left[\begin{array}[]{cc}{\cal
K}^{EE}_{lh}&{\cal K}^{EF}_{lh}\\\\[5.69046pt] {\cal K}^{FE}_{lh}&{\cal
K}^{FF}_{lh}\end{array}\right]\left[\begin{array}[]{c}E^{lh}_{\rm
PS}(P)\\\\[5.69046pt] F^{lh}_{\rm PS}(P)\end{array}\right],\hskip 14.22636pt$
(16)

where

$\displaystyle{\cal K}^{EE}_{lh}$ $\displaystyle=$
$\displaystyle-\int_{q}\,{\rm
Tr}[\gamma_{5}\,\gamma_{\mu}\,S_{l}(q_{+})\,\gamma_{5}\,S_{h}(q_{-})\,\gamma_{\mu}]~{},$
(17) $\displaystyle{\cal K}^{EF}_{lh}$ $\displaystyle=$
$\displaystyle\frac{i}{2M_{lh}}\int_{q}\,{\rm
Tr}[\gamma_{5}\,\gamma_{\mu}\,S_{l}(q_{+})\,\gamma_{5}\,{\gamma\cdot
P}\,S_{h}(q_{-})\,\gamma_{\mu}]~{},$ (18) $\displaystyle{\cal K}^{FE}_{lh}$
$\displaystyle=$ $\displaystyle\frac{2M^{2}_{lh}}{P^{2}}{\cal K}^{EF}_{lh}$
(19) $\displaystyle{\cal K}^{FF}_{lh}$ $\displaystyle=$
$\displaystyle\frac{1}{P^{2}}\int_{q}{\rm Tr}[\gamma_{5}\,{\gamma\cdot
P}\,\gamma_{\mu}\,S_{l}(q_{+})\gamma_{5}\,{\gamma\cdot
P}\,S_{h}(q_{-})\,\gamma_{\mu}].\hskip 54.06006pt$ (20)

Here, the trace is over Dirac indices. All integrals in Eqs. (8) and (17)-(20)
are ultraviolet divergent; the divergences are quadratic and logarithm. The
vast majority of applications within NJL models ignore the pseudo vector
component $F^{lh}_{\rm PS}(P)$; in doing so, this leads to the random-phase-
approximation (RPA) for the BS equation.

##  3 Symmetry-preserving subtraction scheme

The issue of symmetry violation can be exposed examining the momenta running
in the quark propagators; they are $q_{+}=q+k_{1}$ and $q_{-}=q+k_{2}$, with
$k_{1}=\eta_{+}P$ and $k_{2}=-\eta_{-}P$. The momenta $k_{1}$ and $k_{2}$ are
arbitrary, as while satisfying $\eta_{+}+\eta_{-}=1$, the $\eta_{\pm}$ are
otherwise arbitrary. However, to maintain translation invariance, the
integrals must depend only on the relative momentum $k_{1}-k_{2}$ or,
equivalently, the integrals must not depend on $\eta_{\pm}$, they can depend
only in the combination $\eta_{+}+\eta_{-}=1$. The dependence on $k_{1}-k_{2}$
only is also crucial for maintaining the WGT identities, like the one in Eq.
(3).

The traditional way of handling within NJL-type of models integrals like those
in Eqs. (17)-(20) is as follows: after evaluating the traces, a choice for
$\eta_{+}$ and $\eta_{-}$ is made and Feynman parameters are used to combine
in a single term the product $(q^{2}_{+}+M^{2}_{l})(q^{2}_{-}+M^{2}_{h})$ in
the denominator; then, after a momentum shift to eliminate the angle $q\cdot
P$ in the denominator, the integral over $q$ is performed. In making the
momentum shift, changes in the integration limits are ignored. Invariably,
results depend upon the choices made for $\eta_{\pm}$; in particular, the
value of the leptonic decay constants, which play a key role in the
normalisation of the pion BSE, depend upon the choices made for $\eta_{\pm}$.
Some regularization-independent results, like the Goldberg-Treiman relation at
the quark level and the Gell-Mann–Oakes-Renner relation, can be obtained after
using the gap equation to eliminate the quadratic divergences.

The subtraction scheme of Refs. Batt-thesis ; Battistel:1998tj ;
Battistel:2008fd ; Battistel:2013cja ; Farias:2006cs ; Farias:2005cr ;
Farias:2007zz  is based on the repeated use of the identity

$\displaystyle\frac{1}{q^{2}_{\pm}+M^{2}}$ $\displaystyle=$
$\displaystyle\left(\frac{1}{q^{2}_{\pm}+M^{2}}-\frac{1}{q^{2}+M^{2}}\right)+\frac{1}{q^{2}+M^{2}}$
(21) $\displaystyle=$
$\displaystyle\frac{1}{q^{2}+M^{2}}-\frac{q^{2}_{\pm}-q^{2}}{\left(q^{2}+M^{2}\right)}\frac{1}{\left(q^{2}_{\pm}+M^{2}\right)}.$

Assuming a Poincaré invariant regularization for the integrals in Eqs.
(17)-(20), subtractions are performed for each of the propagators
$S_{l}(q_{+})$ and $S(q_{-})$, the number of which is dictated by the
requirement that a finite integral is obtained - note that while the original
denominator goes as $1/q^{2}$ for $q\rightarrow\infty$, the last term in Eq.
(21) goes as $1/q^{3}$ for $q\rightarrow\infty$. We illustrate the procedure
in detail for the ${\cal K}^{EE}_{lh}(P)$ kernel. Evaluation of the trace in
Eq. (17) leads to

$\displaystyle{\cal K}^{EE}_{lh}$ $\displaystyle=$ $\displaystyle
16\int^{\Lambda}_{q}\frac{q_{+}\cdot
q_{-}+M_{l}M_{h}}{(q^{2}_{+}+M^{2}_{l})(q^{2}_{-}+M^{2}_{h})},$ (22)

where $\Lambda$ indicates the mass scale associated with the regularization.
Next, subtracting thrice each of the denominators, ${\cal K}^{EE}_{lh}$ can be
written as a sum of three kinds of terms: quadratic and logarithmically
divergent integrals that are independent of $\eta_{\pm}$, symmetry violating
terms that are proportional to $\eta_{+}$ and $\eta_{-}$, and a finite
integral also independent of $\eta_{\pm}$. Explicitly, they are given by the
following expressions:

$\displaystyle{\cal K}^{EE}_{\rm PS}$ $\displaystyle=$ $\displaystyle
8\Bigl{\\{}\left[\eta^{2}_{+}\,A_{\mu\nu}(M^{2}_{l})+\eta^{2}_{-}\,A_{\mu\nu}(M^{2}_{h})\right]P_{\mu}P_{\nu}+\,I_{\rm
quad}(M^{2}_{l})+I_{\rm quad}(M^{2}_{h})$ (23)
$\displaystyle-\left[P^{2}+\left(M_{h}-M_{l}\right)^{2}\right]\,\left[I_{\rm
log}(M^{2}_{l})-Z_{0}(M^{2}_{l},M^{2}_{h},P^{2};M^{2}_{l})\right]\Bigr{\\}},$

where $A_{\mu\nu}(M^{2})$, $I_{\rm quad}(M^{2})$ and $I_{\rm log}(M^{2})$ are
the divergent integrals (when $\Lambda\rightarrow\infty$)

$A_{\mu\nu}(M^{2})=\int^{\Lambda}_{q}\frac{4q_{\mu}q_{\nu}-(q^{2}+M^{2})\delta_{\mu\nu}}{(q^{2}+M^{2})^{3}},$
(24) $\displaystyle I_{\rm
quad}(M^{2})=\int^{\Lambda}_{q}\,\frac{1}{q^{2}+M^{2}},\hskip 42.67912ptI_{\rm
log}(M^{2})=\int^{\Lambda}_{q}\,\frac{1}{(q^{2}+M^{2})^{2}}~{},$ (25)

and $Z_{0}(M^{2}_{l},M^{2}_{h},P^{2};M^{2}_{l})$ is a finite integral given by

$\displaystyle
Z_{0}(M^{2}_{l},M^{2}_{h},P^{2};M^{2}_{l})=\frac{1}{(4\pi)^{2}}\int^{1}_{0}\,dz\,\ln\left[\frac{P^{2}z(1-z)-(M^{2}_{l}-M^{2}_{h})z+M^{2}_{l}}{M^{2}_{l}}\right].$
(26)

It is important to note that to arrive at these results, no momentum shifts
were made in divergent integrals. The other amplitudes, ${\cal K}^{EF}_{lh}$,
${\cal K}^{FE}_{lh}$, and ${\cal K}^{FF}_{lh}$ are given by similar integrals
SABK .

Note that Eq. (23) makes it clear that whatever choice made for $\eta_{\pm}$,
unavoidably translation symmetry is broken, unless the regularization scheme
leads to $A_{\mu\nu}(M^{2})=0$. Let us next examine how choices of
$\eta_{\pm}$ lead to violation of the WGT identity in Eq. (3). One needs to
deal with integrals of the form

$\displaystyle\int^{\Lambda}_{q}\frac{q_{\mu}}{q^{2}_{\pm}+M^{2}}\\!$
$\displaystyle=$ $\displaystyle\\!\mp\,\eta_{\pm}P_{\mu}\,I_{\rm
quad}(M^{2})\pm\,\frac{1}{3}\eta^{3}_{\pm}\left[P^{2}A_{\mu\rho}(M^{2})P_{\rho}-P_{\mu}A_{\rho\sigma}(M^{2})P_{\rho}P_{\sigma}\right]$
(27)
$\displaystyle\mp\,\eta_{\pm}B_{\mu\rho}(M^{2})P_{\rho}\mp\frac{1}{3}\eta^{3}_{\pm}C_{\mu\rho\sigma\lambda}(M^{2})P_{\sigma}P_{\rho}P_{\lambda},$

where $B_{\mu\nu}(M^{2})$ and $C_{\mu\nu\rho\sigma}(M^{2})$ are the new
structures

$\displaystyle
B_{\mu\nu}(M^{2})=\int^{\Lambda}_{q}\frac{2q_{\mu}q_{\nu}-(q^{2}+M^{2})\delta_{\mu\nu}}{(q^{2}+M^{2})^{2}},\hskip
28.45274ptC_{\mu\nu\rho\sigma}(M^{2})=\int^{\Lambda}_{q}\frac{c_{\mu\nu\rho\sigma}(q,M^{2})}{(q^{2}+M^{2})^{4}},$
(28)

with

$\displaystyle
c_{\mu\nu\rho\sigma}(q^{2},M^{2})=24q_{\mu}q_{\nu}q_{\rho}q_{\sigma}-\,4(q^{2}+M^{2})(\delta_{\mu\nu}q_{\rho}q_{\sigma}+{\rm
perm.}\;\nu\sigma\rho).$ (29)

Using these integrals in the WGT identity, Eq. (3), one obtains

$\displaystyle
0=(M_{l}-m_{l})+(M_{h}-m_{h})-\,\frac{64\pi\alpha_{IR}}{3m^{2}_{G}}\bigg{\\{}I_{\rm
quad}(M^{2}_{l})+I_{\rm
quad}(M^{2}_{h})+\left[\eta^{2}_{+}\,A_{\mu\nu}(M^{2}_{l})+\eta^{2}_{-}\,A_{\mu\nu}(M^{2}_{h})\right]\bigg{\\}}P_{\mu}P_{\nu},$
(30)

and

$\displaystyle 0=\int^{\Lambda}_{q}\left(\frac{q_{+}\cdot
P}{q^{2}_{+}+M^{2}}-\frac{q_{-}\cdot P}{q^{2}_{-}+M^{2}}\right)\sim{\rm
terms}\;\;{\rm proportional}\;\;{\rm
to}\;\;\eta_{\pm}\left(A_{\mu\nu},B_{\mu\nu},C_{\mu\nu\rho\sigma}\right).$
(31)

We see that for arbitrary momentum routing in the loop integrals, i.e. for
arbitrary values of $\eta_{\pm}$ satisfying $\eta_{+}+\eta_{-}=1$, the
subtraction scheme allows to identify in a systematic way symmetry offending
terms; they are the integrals $A_{\mu\nu}$, $B_{\mu\nu}$ and
$C_{\mu\nu\rho\sigma}$ in the equations above. A consistent regularization
scheme should make the integrals vanish automatically. Otherwise, the
vanishing of the integrals must be imposed; in doing so, the regularization
scheme becomes part of the model. Dimensional regularization and Pauli-Villars
regularization are examples of schemes that lead to $A_{\mu\nu}=0$,
$B_{\mu\nu}=0$, and $C_{\mu\nu\rho\sigma}=0$. Removing the symmetry offending
terms, the kernels ${\cal K}^{EE}_{lh}$, ${\cal K}^{EF}_{lh}$, ${\cal
K}^{FE}_{lh}$, and ${\cal K}^{FF}_{lh}$ become free from ambiguities and
symmetry-preserving.

Let us make contact with the scheme of Ref. GutierrezGuerrero:2010md . For
$M_{l}=M_{h}=M$ and $\eta_{+}=1$ and $\eta_{-}=0$, which is the choice made in
Ref. GutierrezGuerrero:2010md , Eq. (31) becomes

$\displaystyle\int^{\Lambda}_{q}\left(\frac{q_{+}\cdot
P}{q^{2}_{+}+M^{2}}-\frac{q\cdot P}{q^{2}+M^{2}}\right)=0.$ (32)

This is precisely Eq. (15) of Ref. GutierrezGuerrero:2010md . Using Eq. (27)
for this integral, one obtains

$\displaystyle 0=\int^{\Lambda}_{q}\left(\frac{q_{+}\cdot
P}{q^{2}_{+}+M^{2}}-\frac{q\cdot
P}{q^{2}+M^{2}}\right)=P_{\mu}\bigl{[}B_{\mu\nu}(M^{2})+\frac{1}{3}C_{\mu\nu\rho\sigma}(M^{2})P_{\rho}P_{\sigma}+\frac{1}{3}P_{\mu}A_{\rho\nu}(M^{2})P_{\rho}-\frac{4}{3}P^{2}A_{\mu\nu}(M^{2})\bigr{]}P_{\nu},$
(33)

that is

$\displaystyle
B_{\mu\nu}(M^{2})+\frac{1}{3}C_{\mu\nu\rho\sigma}(M^{2})P_{\rho}P_{\sigma}+\,\frac{1}{3}P_{\mu}A_{\rho\nu}(M^{2})P_{\rho}-\frac{4}{3}P^{2}A_{\mu\nu}(M^{2})=0.$
(34)

In the quiral limit, $P=0$, this leads to

$B_{\mu\nu}(M^{2})=0=\delta_{\mu\nu}\int_{\Lambda}\frac{d^{4}q}{(2\pi)^{4}}\frac{\frac{1}{2}q^{2}+M^{2}}{(q^{2}+M^{2})^{2}},$
(35)

which is Eq. (17) of Ref. GutierrezGuerrero:2010md . This result makes it
clear that our scheme, besides being in agreement with Ref.
GutierrezGuerrero:2010md  for the particular choice of $\eta_{\pm}$, it is
also more general, as it is valid for arbitary values $\eta_{\pm}$.

##  4 Numerical results

In order assess the reliability of the subtraction scheme, we apply it to the
calculation of the masses and electroweak decay constants of the $\pi$ and $K$
mesons. The mass of the meson is obtained from Eq. (16); it is an eigenvalue
problem which has a solution at $P^{2}=-m^{2}_{\rm PS}$. We use the canonical
normalization condition for the BS amplitudes, namely

$\displaystyle 1=\frac{\partial}{\partial P^{2}}\Pi^{lh}_{\rm
PS}(Q,P)\bigg{|}_{Q=P},$ (36)

with

$\displaystyle\Pi^{lh}_{\rm
PS}(Q,P)=6\int\frac{d^{4}q}{(2\pi)^{4}}{\rm{Tr}}[\bar{\Gamma}^{hl}_{\rm
PS}(-Q)S_{l}(q+k_{1})\Gamma^{lh}_{\rm PS}(Q)S(q+k_{2})]~{},$ (37)

where $\bar{\Gamma}^{hl}_{\rm PS}(-Q)=\left[C^{-1}\Gamma^{lh}_{\rm
PS}(-Q)C\right]^{T}=\Gamma^{lh}_{\rm PS}(-Q)$. The electroweak decay constant
of the meson, $f_{\rm PS}$, can be expressed in terms of normalized amplitudes
$E$ and $F$ in the form

$\displaystyle f_{\rm PS}=\frac{6}{M_{lh}}\left[E_{\rm PS}(P){\cal
K}^{FE}_{\rm PS}+F_{\rm PS}(P){\cal K}^{FF}_{\rm
PS}\right]\bigg{|}_{P^{2}=-m^{2}_{\rm PS}}.$ (38)

We compare results with Ref. Chen:2012txa , where the scheme of Ref.
GutierrezGuerrero:2010md  was used for the $K$ meson. We use proper-time
regularization for $I^{\Lambda}_{quad}(M_{f})$. In the BS amplitudes where
$I^{\Lambda}_{quad}(M_{f})$ appears, we use the gap equation (8), and for
$I^{\Lambda}_{log}(M_{f})$ use the identity

$I^{\Lambda}_{log}(M_{f})=-\frac{dI^{\Lambda}_{quad}(M_{f})}{dM^{2}_{f}}.$
(39)

On the other hand, for the subtraction mass scale we have used the light mass
$M^{2}_{l}$ and finite integrals like
$Z_{0}(M^{2}_{h},M^{2}_{l},P^{2};M^{2}_{h})$ are performed without imposing a
regulator.

We use the same parameters employed in Ref. Chen:2012txa :
$\alpha_{IR}=0.93\pi$, $m_{G}=0.8$ GeV, $m_{u}=m_{d}=0.007$ GeV and
$m_{s}=0.160$ GeV, $\tau^{2}_{IR}=(0.24\,{\rm GeV})^{2}$ and
$\tau^{2}_{UV}=(0.905\,{\rm GeV})^{2}$. The last two parameters are the
infrared and ultraviolet cutoff parameters of the proper-time regularization.
The solution of the gap equation leads for the constituent quark masses
$M_{u}=M_{d}=0.367$ GeV, and $M_{s}=0.525$ GeV. Table I presents our results
and those of Ref. Chen:2012txa  (indicated with a star).

Table 1: Masses and electroweak decay constants$\pi$ and $K$ mesons. The four first columns show the results obtained using the subtraction scheme and the four last indicated with a star, are the results of Ref. Chen:2012txa . |  Meson |  $E$ |  $F$ |  $m_{\rm PS}$ |  $f_{\rm PS}$ |  $E^{*}$ |  $F^{*}$ |  $m^{*}_{\rm PS}$ |  $f^{*}_{\rm PS}$
---|---|---|---|---|---|---|---|---|---
| $\pi$ | 3.759 | 0.498 | 0.139 | 0.106 | 3.596 | 0.474 | 0.139 | 0.101
| $K$ | 3.984 | 0.632 | 0.494 | 0.115 | 3.864 | 0.591 | 0.493 | 0.107

We see from Table I that the results obtained with the subtraction scheme
compares well with those of Ref. Chen:2012txa : the masses are almost
identical in both approaches, and the decay constants from the subtraction are
a little larger, but by less than 5%. We recall that the differences in the
approaches are that in the subtraction scheme the results are independent of
the choice of $\eta_{\pm}$, while those of Ref. Chen:2012txa  are for
$\eta_{+}=1$ and $\eta_{-}=0$. Another difference is that the finite integral
$Z_{0}(M^{2}_{h},M^{2}_{l},P^{2};M^{2}_{h})$ is integrated without imposing
ultraviolet or infrared cutoffs. There is no difficulty in using an infrared
cutoff in the finite integrals to avoid unphysical quark-antiquark thresholds
in amplitudes.

##  5 Conclusions and Perspectives

We examined the contact-interaction model introduced in Ref.
GutierrezGuerrero:2010md  within the perspective of a regularization scheme
that allows to separate symmetry-offending parts in Bethe-Salpeter amplitudes
in a way independent of choices of momentum routing in divergent integrals. In
doing so, the Ward-Green-Takahashi (WGT) identities reflecting global
symmetries of the model are preserved by the regularization. Symmetry-
offending parts of the amplitudes, the integrals $A_{\mu\nu}$, $B_{\mu\nu}$
and $C_{\mu\nu}$, can be neatly separated. In general, a cutoff regularization
scheme leads to nonzero values for the symmetry-offending integrals, while
Pauli-Villars or dimensional regularization lead to the vanishing of the
symmetry-offending integrals. In a nonrenormalizable model, like the contact-
interaction model discussed here, the vanishing of $A_{\mu\nu}$, $B_{\mu\nu}$
and $C_{\mu\nu}$ must be imposed in an ad hoc manner.

Our aim is to apply the subtraction scheme to heavy-light mesons, where the
heavy quark is much heavier than the light quark, like in the $D$ and $B$
mesons. When the masses are not much different, one can show that in the RPA
approximation to the BS equation, the heavy-light pseudoscalar meson mass
$m^{2}_{\rm PS}$ and corresponding electroweak decay constant $f_{\rm PS}$ can
be expressed in terms of the light-quark condensate $\langle\bar{q}q\rangle$
and the pion decay constant $f_{\pi}$ as

$\displaystyle m^{2}_{\rm PS}$ $\displaystyle\simeq$
$\displaystyle(M_{h}-M_{l})^{2}-m_{h}\frac{M_{l}}{M_{h}}\frac{\langle\bar{q}q\rangle}{f^{2}_{\pi}}$
(40) $\displaystyle f_{\rm PS}$ $\displaystyle\simeq$
$\displaystyle\frac{1}{2}\left(1+\frac{M_{h}}{M_{l}}\right)f_{\pi}$ (41)

These expressions are valid for $M_{h}\leq\Lambda$, where $\Lambda$ is the
regularization mass scale in the divergent integrals $I_{\rm quad}$ and
$I_{\rm log}$. Eq. (41) shows very clearly that $f_{\pi}<f_{K}<f_{D}$, which
is the correct ordering of their experimental values. Using the traditional
cutoff regularization scheme, in which no attention is given to symmetry-
offending terms in amplitudes, one obtains blaschke  $f_{\pi}\simeq f_{K}$ and
$f_{D}<f_{\pi}$. Of course, Eqs. (40) and (41) are not strictly applicable to
this case, but they show the trend one can expect when using the subtraction
scheme. Because of the large diferences between the masses of the quarks
running in the loops, the application of the subtraction method in this case
is nontrivial and needs adaptations – results will be presented elsewhere SABK
.

The work of M.A.B. and F.E.S. were financed by doctoral scholarships by
Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - CAPES and
Conselho Nacional de Desenvolvimento Científico e Tecnológico - CNPq,
respecytively. The work of G.K. was supported in part by Conselho Nacional de
Desenvolvimento Científico e Tecnológico \- CNPq, Grant No. 305894/2009-9, and
Fundação de Amparo à Pesquisa do Estado de São Paulo - FAPESP, Grant No.
2013/01907-0.

## References

  * (1) Y. Nambu and G. Jona-Lasinio,  Phys. Rev. 122, 345 (1961).
  * (2) U. Vogl and W. Weise,  Prog. Part. Nucl. Phys. 27, 195 (1991).
  * (3) S. Klevansky,  Rev. Mod. Phys. 64, 649 (1992).
  * (4) T. Hatsuda and T. Kunihiro,  Phys. Rept. 247, 221 (1994).
  * (5) J. Bijnens,  Phys. Rept. 265, 369 (1996).
  * (6) L. X. Gutiérrez-Guerrero, A. Bashir, I. C. Cloët and C. D. Roberts,  Phys. Rev. C 81, 065202 (2010).
  * (7) H. L. L. Roberts, C. D. Roberts, A. Bashir, L. X. Gutiérrez-Guerrero and P. C. Tandy,  Phys. Rev. C 82, 065202 (2010).
  * (8) H. L. L. Roberts, A. Bashir, L. X. Gutiérrez-Guerrero, C. D. Roberts and D. J. Wilson,  Phys. Rev. C 83, 065206 (2011).
  * (9) H. J. Munczek,  Phys. Rev. D 52, 4736 (1995).
  * (10) A. Bender, C. D. Roberts and L. von Smekal,  Phys. Lett. B 380, 7 (1996).
  * (11) I. C. Cloet and C. D. Roberts, Prog. Part. Nucl. Phys. 77, 1 (2014).
  * (12) O. A Battistel, PhD Thesis 1999, Universidade Federal de Minas Gerais, Brazil (unpublished).
  * (13) O. A. Battistel and M. C. Nemes, Phys. Rev. D 59, 055010 (1999).
  * (14) O. A. Battistel, G. Dallabona and G. Krein, Phys. Rev. D 77, 065025 (2008) .
  * (15) O. A. Battistel and G. Dallabona, Phys. Rev. D 80, 085028 (2009).
  * (16) R. L. S. Farias, G. Dallabona, G. Krein and O. A. Battistel, Phys. Rev. C 77, 065201 (2008).
  * (17) R. L. S. Farias, G. Dallabona, G. Krein and O. A. Battistel, Phys. Rev. C 73, 018201 (2006).
  * (18) R. L. S. Farias, G. Krein, G. Dallabona and O. A. Battistel, Nucl. Phys. A 790, 332 (2007).
  * (19) F. E. Serna , M. A. Brito and G. Krein, in preparation.
  * (20) C. Chen, L. Chang, C. D. Roberts, S. M. Schmidt, S. Wan and D. J. Wilson,  Phys. Rev. C 87, 045207 (2013).
  * (21) D. Blaschke, P. Costa, and Yu. L. Kalinovsky, Phys. Rev. D 85, 034005 (2012).