Units and calendar attributes of time_bnds are dropped by to_netcdf

What happened?

If I define a time_bnds variable and denote it with the bounds attribute to the time variable, sometimes to_netcdf will drop the units and calendar attributes of the time_bnds variable. This seems like a bug to me.

What did you expect to happen?

I expected written variables to have the attributes I defined for them.

Minimal Complete Verifiable Example

# /// script
# requires-python = ">=3.11"
# dependencies = [
#   "xarray[complete]@git+https://github.com/pydata/xarray.git@main",
# ]
# ///
#
# This script automatically imports the development branch of xarray to check for issues.
# Please delete this header if you have _not_ tested this script with `uv run`!

import xarray as xr
xr.show_versions()
# your reproducer code ...
import numpy as np
from datetime import datetime

TIME0 = datetime(1980, 1, 1)
TUNITS = f'days since {TIME0.strftime("%Y-%m-%d")}'
TCAL = 'proleptic_gregorian'

year = 2001
tval = (datetime(year, 1, 1) - TIME0).days
yrdays = (datetime(year+1, 1, 1) - datetime(year, 1, 1)).days
tbvals = np.reshape([tval, tval + yrdays], (1, 2))
time_bnds = xr.DataArray(
    data=tbvals.astype(np.double), dims=['time','nv'],
    attrs={'long_name':'time bounds', 'units':TUNITS, 'calendar':TCAL}
)

nlat = 180
nlon = 360
late = np.linspace( -90,  90, nlat+1)
lone = np.linspace(-180, 180, nlon+1)
lat = 0.5*(late[1:] + late[:-1])
lon = 0.5*(lone[1:] + lone[:-1])

coords = {
    'time':(['time'], np.array([tval]).astype(np.double), {
        'long_name':'time',
        'units':TUNITS,
        'calendar':TCAL,
        'bounds':'time_bnds',
    }),
    'lat':(['lat'], lat.astype(np.single), {
        'long_name':'latitude',
        'units':'degrees_north',
    }),
    'lon':(['lon'], lon.astype(np.single), {
        'long_name':'longitude',
        'units':'degrees_east',
    }),
}

field = xr.DataArray(
    data=np.zeros((1, nlat, nlon)).astype(np.single),
    dims=['time','lat','lon'], coords=coords,
)

ds1 = xr.Dataset(data_vars={'time_bnds':time_bnds})
ds1.to_netcdf('test1.nc4')
ds1.close()

ds2 = xr.Dataset(data_vars={'time_bnds':time_bnds, 'field':field})
ds2.to_netcdf('test2.nc4')
ds2.close()

ds3 = xr.Dataset(data_vars={'something':time_bnds, 'field':field})
ds3.to_netcdf('test3.nc4')
ds3.close()

Steps to reproduce

No response

MVCE confirmation

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• New issue — a search of GitHub Issues suggests this is not a duplicate.
• Recent environment — the issue occurs with the latest version of xarray and its dependencies.

Relevant log output

(newpy) bweir@discover11:/discover/nobackup/bweir> ncdump -h test1.nc4 
netcdf test1 {
dimensions:
	time = 1 ;
	nv = 2 ;
variables:
	double time_bnds(time, nv) ;
		time_bnds:_FillValue = NaN ;
		time_bnds:long_name = "time bounds" ;
		time_bnds:units = "days since 1980-01-01" ;
		time_bnds:calendar = "proleptic_gregorian" ;
}
(newpy) bweir@discover11:/discover/nobackup/bweir> ncdump -h test2.nc4 
netcdf test2 {
dimensions:
	time = 1 ;
	nv = 2 ;
	lat = 180 ;
	lon = 360 ;
variables:
	double time(time) ;
		time:_FillValue = NaN ;
		time:long_name = "time" ;
		time:units = "days since 1980-01-01" ;
		time:calendar = "proleptic_gregorian" ;
		time:bounds = "time_bnds" ;
	double time_bnds(time, nv) ;
		time_bnds:_FillValue = NaN ;
		time_bnds:long_name = "time bounds" ;
	float lat(lat) ;
		lat:_FillValue = NaNf ;
		lat:long_name = "latitude" ;
		lat:units = "degrees_north" ;
	float lon(lon) ;
		lon:_FillValue = NaNf ;
		lon:long_name = "longitude" ;
		lon:units = "degrees_east" ;
	float field(time, lat, lon) ;
		field:_FillValue = NaNf ;
}
(newpy) bweir@discover11:/discover/nobackup/bweir> ncdump -h test3.nc4 
netcdf test3 {
dimensions:
	time = 1 ;
	nv = 2 ;
	lat = 180 ;
	lon = 360 ;
variables:
	double time(time) ;
		time:_FillValue = NaN ;
		time:long_name = "time" ;
		time:units = "days since 1980-01-01" ;
		time:calendar = "proleptic_gregorian" ;
	double something(time, nv) ;
		something:_FillValue = NaN ;
		something:long_name = "time bounds" ;
		something:units = "days since 1980-01-01" ;
		something:calendar = "proleptic_gregorian" ;
	float lat(lat) ;
		lat:_FillValue = NaNf ;
		lat:long_name = "latitude" ;
		lat:units = "degrees_north" ;
	float lon(lon) ;
		lon:_FillValue = NaNf ;
		lon:long_name = "longitude" ;
		lon:units = "degrees_east" ;
	float field(time, lat, lon) ;
		field:_FillValue = NaNf ;
}

Anything else we need to know?

You'll see that test1.nc4 and test3.nc4 have the correct attributes for time_bnds. If you comment out the bounds attribute for the time variable, then time_bnds will have the correct attributes. Clearly there is something going on in the code that's tracing through time to time:bounds to time_bnds:units and time_bnds:calendar and dropping them. Why you would write code to do this horrifies me.

Environment

Start with the provided MVCE and the Dataset.to_netcdf entry point; compare test1.nc4, test2.nc4, and test3.nc4 with ncdump. Trace the bounds-related encoding path and add a regression test showing that units and calendar remain on time_bnds in the test2.nc4 case. Done means the reproduced output preserves both attributes.