Add Support for Symmetric Interior Penalty (SIP) Method for Parabolic Terms in Trixi.jl

I am currently working on implementing the advection-diffusion equation in Trixi.jl. Upon reviewing the source code, I have observed that the semidiscretization for parabolic terms is handled within `solvers/dgsem_tree/dg_1d_parabolic.jl`. Specifically, the function `rhs_parabolic!` utilizes arguments such as `dg::DG` for the solver and `parabolic_scheme`, although it appears that `parabolic_scheme` doesn't play any role in the current implementation. The diffusion term seems to default to the classical **Bassi-Rebay 1 (BR1)** formulation.

While the BR1 formulation is currently supported and the LDG formulation remains on the TO-DO list, I would like to implement the **Symmetric Interior Penalty (SIP)** method for diffusion terms, as described in https://doi.org/10.1007/978-3-642-22980-0.

**Proposed Changes:**
1. Add support for the **Symmetric Interior Penalty (SIP)** method within the `solvers/dgsem_tree/dg_1d_parabolic.jl` module.
2. Update the `parabolic_scheme` argument to recognize and handle this new formulation, allowing users to select between BR1, LDG, and SIP formulations.

Start in solvers/dgsem_tree/dg_1d_parabolic.jl, especially rhs_parabolic!, and trace how parabolic_scheme currently relates to the BR1 formulation. Determine the changes needed for SIP and scheme selection alongside BR1 and the planned LDG option. Done means users can select SIP through parabolic_scheme and the parabolic diffusion implementation supports it.

Written by the indexing model from the issue text.